# Primero creamos el mapa de toda la instrumentación

rm(list=ls()) # remove all the variables from the workspace
setwd("C:/Users/Marcos/Documents/TFM/Situacion_pluviometros/")
library(ggmap)
library(rgdal)
library(ggplot2)
library(scales)
library(maptools)

# Convertimos coordenadas de UTM a geográficas:

# Leemos el archivo de texto donde están las coordenadas UTM de los Pluviómetros
coord <- read.delim("C:/Users/Marcos/Documents/TFM/Situacion_pluviometros/pluviometrosVenero.txt", header=TRUE)[,4:5]

# Añadimos coordenadas del TDR, los dos limnímetros y el radar meteorológico (figura 3.1 TFM_Carlos)
coord.todo <- rbind(coord,c(356489,4474915), c(360744,4470212), c(360040,4471782), 
                    c(359088,4474084))
rownames(coord.todo)=make.names(c(rep("Pluviómetros",7), "Radar", "TDR", rep("Limnímetros",2)), unique=TRUE)
coord.todo.spdf <- coord.todo  # ".spdf": Class for spatial attributes that have spatial point locations
colnames(coord.todo.spdf) <- c("x","y")
coordinates(coord.todo.spdf) <- ~x+y

# To retrieve the CRS for a spatial object: proj4string(x)
proj4string(coord.todo.spdf) <- CRS("+proj=utm +zone=30 ellps=WGS84")

# Ahora pasamos de un CRS a otro (pasamos de coordenadas UTM a geográficas)
coord.geog <- as.data.frame(spTransform(coord.todo.spdf, CRS("+proj=longlat +datum=WGS84")))
colnames(coord.geog) <- c("lon", "lat")

# Creamos un data.frame llamado "datos" para poder pintar bien la leyenda del mapa
datos <- coord.todo
datos$lon <- coord.geog$lon
datos$lat <- coord.geog$lat
datos$INSTRUMENTACIÓN <- c("Pluviómetro","Pluviómetro","Pluviómetro","Pluviómetro","Pluviómetro","Pluviómetro","Pluviómetro",
                           "Radar meteorológico","TDR","Limnímetro","Limnímetro")

# Ahora con la función "get_map" consultamos en Google Maps el mapa 
# que queremos crear de nuestra zona de estudio 
mapa <- get_map(location=c(lon=-4.66, lat=40.39), zoom=13, maptype="terrain", 
                source="google")
nombres.pluvio <- c("Arromoro", "Atalaya", "Dehesa", "Gasolinera", "Collado Morales",
                    "Peña Parda", "Trampalones")

# Cercamos el área de la cuenca
area <- readShapePoly("venero_Project.shp")     # Lee el shape con el area de la cuenca
# To change to correct projection:
proj4string(area) <- "+proj=utm +zone=30 +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
area <- spTransform(area, CRS("+proj=longlat +datum=WGS84"))


png("mapa_zona_estudio.png", units="in", width=9, height=8.5, res=300)

# Vamos a pintar los distintos instrumentos de medida sobre el mapa (junto con éste) y añadimos leyenda
ggmap(mapa, extend='device', legend="left", base_layer=ggplot(datos, aes(x=lon, y=lat))) +
  geom_point(data=datos, aes(shape = INSTRUMENTACIÓN, color=INSTRUMENTACIÓN, fill=INSTRUMENTACIÓN), size=7.5) +
  geom_text(data=coord.geog[1:7,], aes(x=coord.geog$lon[1:7], y=coord.geog$lat[1:7], label=nombres.pluvio, hjust=0.2, vjust=1), size=5.2) +
  scale_shape_manual(values=c(21,24,22,23)) +
  geom_polygon(aes(x=long, y=lat, group=id), data=area, colour="red", alpha=1.0, size=1.1)

dev.off()