Mapa de referencia en R

Maestría en Ecohidrología
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/

Daniela Ballari (PhD)
dballari@uazuay.edu.ec
Universidad del Azuay
Publicaciones - https://scholar.google.com/citations?user=1VcAm9AAAAAJ

Objetivo de la práctica: Crear un mapa de referencia para publicar. El mapa de referencia documentado aqui fue publicado en: Ballari D, Giraldo R, Campozano L, Samaniego E. Spatial functional data analysis for regionalizing precipitation seasonality and intensity in a sparsely monitored region: Unveiling the spatio-temporal dependencies of precipitation in Ecuador. Int J Climatol. 2018;1-18. https://doi.org/10.1002/joc.5504

La figura es la siguiente. Se encuentra compuesta de 3 figuras.

Cargar librerias y datos

library(rgdal) # Spatial package

## Loading required package: sp

## rgdal: version: 1.3-3, (SVN revision 759)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/Usuario/Documents/R/win-library/3.4/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Usuario/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.3-1

library(sp) 
library(sf)

## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3

library(ggmap) # Visualization

## Loading required package: ggplot2

library(ggsn) # Visualization
library(ggplot2) # Visualization
library(gridExtra) # Panel visualization #grid.arrange(p1, p2, ncol=2)

setwd("C:/R_ecohidrologia/Mapa_referencia") #definir directorio de trabajo 

Preparación de datos

Los datos a utilizar en esta práctica están disponibles Aqui

#Study area
ecuador<- readOGR(".","ecuador")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\R_ecohidrologia\Mapa_referencia", layer: "ecuador"
## with 1 features
## It has 7 fields

zona_estudio.f <- fortify(ecuador) # ggplot visualization

# data used for TRMM correction
prec_raingauge_correction<- readOGR(".","est_corr")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\R_ecohidrologia\Mapa_referencia", layer: "est_corr"
## with 1808 features
## It has 1 fields

#Raingauges for validation. INAMHI 2011 without the ones used for trmm correction
raingauges_validation<- readOGR(".","est_val")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\R_ecohidrologia\Mapa_referencia", layer: "est_val"
## with 145 features
## It has 1 fields

#Altitude of 1000m and DEM for visualization  
cn<- readOGR(".","curva_nivel_1000")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\R_ecohidrologia\Mapa_referencia", layer: "curva_nivel_1000"
## with 1 features
## It has 2 fields

cn.f <- fortify(cn)

#DEM, digital elevation model
load("dem4326.RData") #dem4326.df

Visualización de capas

plot(ecuador)

plot(prec_raingauge_correction)

plot(raingauges_validation)

plot(cn)

Figura a. Reference map

sq_map <- get_map(c(-125,-60,-25,50), zoom = 3, source = "google", maptype = "watercolor",  color = "bw")

## maptype = "watercolor" is only available with source = "stamen".

## resetting to source = "stamen"...

## Map from URL : http://tile.stamen.com/watercolor/3/1/2.jpg

## Map from URL : http://tile.stamen.com/watercolor/3/2/2.jpg

## Map from URL : http://tile.stamen.com/watercolor/3/3/2.jpg

## Map from URL : http://tile.stamen.com/watercolor/3/1/3.jpg

## Map from URL : http://tile.stamen.com/watercolor/3/2/3.jpg

## Map from URL : http://tile.stamen.com/watercolor/3/3/3.jpg

## Map from URL : http://tile.stamen.com/watercolor/3/1/4.jpg

## Map from URL : http://tile.stamen.com/watercolor/3/2/4.jpg

## Map from URL : http://tile.stamen.com/watercolor/3/3/4.jpg

## Map from URL : http://tile.stamen.com/watercolor/3/1/5.jpg

## Map from URL : http://tile.stamen.com/watercolor/3/2/5.jpg

## Map from URL : http://tile.stamen.com/watercolor/3/3/5.jpg

p1<- ggmap(sq_map, height=600) + annotate("rect", xmin = -75, xmax = -82, ymin = -6, ymax = 2, alpha = .5,colour="black", fill="grey70")+
  theme(text = element_text(size=6))+
  annotate("text", x = -121, y = 48, label = "(a)", size=2)+  
  annotate("text", x = -50, y = -13, label = "South ", fontface ="italic", size=2)+
  annotate("text", x = -50, y = -18, label = "America", fontface ="italic", size=2)+
  annotate("text", x = -63, y = 0, label = "Ecuador", fontface =1, size=2)+
  annotate("text", x = -95, y = -20, label = "Pacific Ocean", fontface ="italic", size=2)+
  annotate("text", x = -45, y = 23, label = "Atlantic Ocean", fontface ="italic", size=2)+  
  coord_fixed(ylim = c(-60, 50))

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Figura b. Altitudinal reference map

sbbox <- make_bbox(lon = zona_estudio.f$lon, lat = zona_estudio.f$lat, f = .1)
sq_map <- get_map(location = sbbox, source = "google", maptype = "watercolor",color = "bw")#

## maptype = "watercolor" is only available with source = "stamen".

## resetting to source = "stamen"...

## Map from URL : http://tile.stamen.com/watercolor/7/34/63.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/35/63.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/36/63.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/37/63.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/34/64.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/35/64.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/36/64.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/37/64.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/34/65.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/35/65.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/36/65.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/37/65.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/34/66.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/35/66.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/36/66.jpg

## Map from URL : http://tile.stamen.com/watercolor/7/37/66.jpg

p2<- ggmap(sq_map, height=600) + 
  geom_raster(data= dem4326.df, aes(s1, s2, fill=dem.ecuador.32717)) +
  scale_fill_gradient2(low = 'white', high = 'black')+ 
  theme(legend.title=element_blank(), legend.justification=c(1,0), legend.position=c(1,0), legend.background = element_rect(alpha('white', 0.0)))+
  theme(text = element_text(size=6), legend.key.height=unit(7,"point"),legend.key.width=unit(5,"point"))+
  geom_polygon(data=zona_estudio.f, aes(x=long, y=lat, group = group, alpha = 1),colour="grey70", fill=NA, show.legend = FALSE)+
  coord_cartesian()+
  annotate("text", x = -81.3, y = 1.95, label = "(b)", size=2)+  
  annotate("rect", xmin = -79.15, xmax = -77.85, ymin = -0.85, ymax = -1.55, alpha = .5, fill="white")+
  annotate("text", x = -80, y = -1, label = "Coast", fontface ="italic", size=2)+  
  annotate("text", x = -80, y = -1.3, label = "Region", fontface ="italic", size=2)+  
  annotate("text", x = -78.5, y = -1, label = "Andean", fontface ="italic", size=2)+ 
  annotate("text", x = -78.5, y = -1.3, label = "Region", fontface ="italic", size=2)+ 
  annotate("text", x = -76.6, y = -1, label = "Amazonian", fontface ="italic", size=2)+  
  annotate("text", x = -76.6, y = -1.3, label = "Region", fontface ="italic", size=2)+  
  coord_fixed()

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Figura c. Rain gauges reference map

sbbox <- make_bbox(lon = zona_estudio.f$lon, lat = zona_estudio.f$lat, f = .1)
sq_map <- get_map(location = sbbox, source = "google", maptype = "watercolor",color = "bw")#

## maptype = "watercolor" is only available with source = "stamen".

## resetting to source = "stamen"...

## Map from URL : http://tile.stamen.com/watercolor/7/34/63.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## 'f3ee4c3b4bcba8c81302339815909cd9.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/35/63.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## '3ad6f1241107c629ade5f2670d145776.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/36/63.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## 'a7c98a7e45cfc3519eb880e4ccfc3fd4.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/37/63.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## '80cc7fd564a73da7120fb3193936dc3e.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/34/64.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## '439e33eeef99af16cf925a1bbbf679cd.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/35/64.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## '8671abdf5dc5edfe0796663b4217637f.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/36/64.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## '9275bbce7efe04e53ebe17b71f713b8b.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/37/64.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## '5a4a5cdbc963a7ddf8d496a3e67e2a5d.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/34/65.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## '843734ec1371055353e93f1e9dd5bc2c.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/35/65.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## '47ce10313d45b0c2cfba2aba5f2ac310.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/36/65.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## 'b3a81690be504e9ed174599995347ab6.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/37/65.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## 'a5f8d52930052078042d577e41e7a067.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/34/66.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## '8d9f26ff0664ce8b12a4ae736bb59d43.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/35/66.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## '154b12e0db0cec46cee58af0a5f626f2.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/36/66.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## '92e47afea7fe619c84fba6fe5a293ead.rds', motivo 'No such file or directory'

## Map from URL : http://tile.stamen.com/watercolor/7/37/66.jpg

## Warning in file.remove(index[[url]]): no fue posible abrir el archivo
## 'fb03f6c3e489dfd9534da8521703cc04.rds', motivo 'No such file or directory'

#joing rain gauges data frames
val <- data.frame(raingauges_validation@coords)
names(val)<- c("X", "Y")   
val$Rain_gauges <- "validation rain gauges"
cor <- data.frame(prec_raingauge_correction@coords)
names(cor)<- c("X", "Y")   
cor$Rain_gauges <- "correction rain gauges"
points <- rbind(val, cor)

p3<- ggmap(sq_map,height=600) + 
  geom_polygon(data=cn.f, aes(x=long, y=lat, group = group, alpha = 1),colour="white", fill=NA,show.legend = FALSE)+
  geom_polygon(data=zona_estudio.f, aes(x=long, y=lat, group = group, alpha = 1),colour="grey70", fill=NA, show.legend = FALSE)+
  geom_point(data=points, aes(x = X, y = Y, shape = type, color = Rain_gauges), shape=16, size=0.7)+
  scale_color_manual(values = c("correction rain gauges" = 'black','validation rain gauges' = 'grey60')) +
  theme(text = element_text(size=6), legend.key.height=unit(5,"point"),legend.key.width=unit(1,"point"),
        legend.title=element_blank(), legend.justification=c(1,0), legend.position=c(1,0), legend.background = element_rect(alpha('white', 0.0)))+
  annotate("text", x = -81.3, y = 1.95, label = "(c)", size=2)+  
  annotate("text", x = -79.5, y = 0.6, label = "Coast", fontface ="italic", size=2)+  
  annotate("text", x = -79.5, y = 0.33, label = "Region", fontface ="italic", size=2)+  
  annotate("text", x = -77.5, y = 1.27, label = "Andean", fontface ="italic", size=2)+ 
  annotate("text", x = -77.52, y = 1.0, label = "Region", fontface ="italic", size=2)+ 
  annotate("text", x = -76.6, y = -1.3, label = "Amazonian", fontface ="italic", size=2)+  
  annotate("text", x = -76.6, y = -1.60, label = "Region", fontface ="italic", size=2)+  
  scalebar(zona_estudio.f, dist = 100, dd2km = TRUE, model = 'WGS84', st.size=1.3, location="bottomleft")+
  coord_fixed()

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Unir figuras y guardar como pdf

pdf("Figure1.pdf")
grid.arrange(p1, p2, p3, ncol=3,respect=TRUE)

## Warning: Removed 4 rows containing missing values (geom_point).

dev.off()

## png 
##   2

## Warning: Removed 4 rows containing missing values (geom_point).

Trabajo Calificado 6p

Adapta el código presentado en esta práctica utilizando la libraría sf en lugar de sp. 2p Crea un mapa de referencia para tu zona de estudio usando este código o Qgis. 4p Reporta los pasos realizados y el resultado final.

Has llegado al final de este material!