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

Daniela Ballari (PhD)
daniela.ballari@ucuenca.edu.ec
http://www.researchgate.net/profile/Daniela_Ballari/publications

Johanna Orellana Alvear (MSc.)
johanna.orellana@ucuenca.edu.ec

Curso completo en: http://rpubs.com/Johanna_Orellana_Alvear/MHidro_indice_2018



En esta sesión revisaremos algunas opciones para visualizar archivos de geoinformación con diferentes librerías como ggplot y ggmap. Fuentes de informaci?n utilizados en este material son: http://rspatial.r-forge.r-project.org/gallery/ y http://www.kevjohnson.org/making-maps-in-r/

Temario

A- Librerías y datos
B- Visualización con ggplot
C- Utiliza mapas de Google

A- Librerías y datos

Carguemos las librerías y datos para generar los mapas.

library(sp)#visualización y gestión archivos vectoriales
library(rgdal) #leer archivos de geoinformación
library(rgeos) #geoprocesos vectoriales
library(ggplot2)#visualización
library(latticeExtra) #complemento visualización
library(ggmap)#complemento visualización, mapas de google
library(ggsn) #escala gráfica y norte para ggplot
#Directorio de trabajo. Aqui deben estar todos los datos almacenados.
setwd("~/Documents/R_Workspace/Maestria_Hidrologia_2018/Dani_clase_adicional/datos/SRS")

#Cargar archivos de geoinformación vectoriales.
zona_estudio<-readOGR(".","prov_zona_estudio")
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "prov_zona_estudio"
## with 7 features
## It has 11 fields
## Integer64 fields read as strings:  DPA_VALOR id
estaciones<-readOGR(".","estaciones_32717") #estaciones meteorológicas
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "estaciones_32717"
## with 40 features
## It has 14 fields
## Integer64 fields read as strings:  ZONA_HIDRO ALTURA PROVINCIA
rios<-readOGR(".","rio_32717_zona_estudio")#fuente IGM 50k
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "rio_32717_zona_estudio"
## with 524 features
## It has 16 fields
## Integer64 fields read as strings:  DPA_VALOR id
#Cargar archivos de geoinformaci?n raster.
elevacion <- readGDAL("elevacion.tif")
## elevacion.tif has GDAL driver GTiff 
## and has 798 rows and 844 columns
names(elevacion)
## [1] "band1"
names(elevacion) <- "dem" # asignar nombre "dem" a campo elevacion
names(elevacion)
## [1] "dem"
#Comprobar que los sistemas de referencias coincidan.
proj4string(zona_estudio)
## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(estaciones)
## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(rios)
## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(elevacion)
## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(zona_estudio)==proj4string(estaciones)
## [1] TRUE
proj4string(zona_estudio)==proj4string(rios)
## [1] TRUE
proj4string(zona_estudio)==proj4string(elevacion)
## [1] TRUE

B- Visualización con ggplot

La libreria ggplot ofrece herramientas muy interesantes para visualización.

#En el caso de polígonos y líneas es necesario aplicar primero la función `fortify` para generar un listado de coordenadas, esto es un `data.frame` util para plotear.
zona_estudio.f <- fortify(zona_estudio, region="id")
head(zona_estudio.f)
##       long     lat order  hole piece id group
## 1 770262.3 9716934     1 FALSE     1  1   1.1
## 2 770300.8 9716975     2 FALSE     1  1   1.1
## 3 770358.2 9717035     3 FALSE     1  1   1.1
## 4 770416.0 9717096     4 FALSE     1  1   1.1
## 5 770422.6 9717101     5 FALSE     1  1   1.1
## 6 770483.7 9717143     6 FALSE     1  1   1.1
#Mapa básico
ggplot() + 
  geom_polygon(data=zona_estudio.f, aes(x=long, y=lat, group = group),colour="white", fill="grey70" )

# Etiquetas de ejes y título del mapa
ggplot() + 
  geom_polygon(data=zona_estudio.f, aes(x=long, y=lat, group = group),colour="white", fill="grey70" ) +
  labs(x="coord x", y="coord y", title="Provincias y estaciones")

# Superponer capa
ggplot() + 
  geom_polygon(data=zona_estudio.f, aes(x=long, y=lat, group = group),colour="white", fill="grey70" ) +
  labs(x="coord x", y="coord y", title="Provincias y estaciones") +
  geom_point(data=estaciones@data, aes(x=x, y=y), color="red")

# Transparencia y mapa proporcional en x e y
ggplot() + 
geom_polygon(data=zona_estudio.f, aes(x=long, y=lat, group = group),colour="white", fill="grey70" , alpha=0.5) +
  labs(x="coord x", y="coord y", title="Provincias y estaciones")+
  geom_point(data=estaciones@data, aes(x=x, y=y), color="red") +
  coord_equal()

#Escala gráfica y norte, libreria ggsn. Más información http://oswaldosantos.github.io/ggsn/
map <- ggplot(data=zona_estudio.f, aes(x=long, y=lat, group=group)) +
  geom_polygon(colour="white", fill="grey70" , alpha=0.5) +
  labs(x="coord x", y="coord y", title="Provincias y estaciones") +
  coord_equal()
  
map + 
  north(zona_estudio.f, symbol=10)+
  scalebar(zona_estudio.f, dist = 100,st.size=4) 

map + 
  blank()+
  north(zona_estudio.f, symbol=10)+
  scalebar(zona_estudio.f, dist = 100,st.size=4)

Para guardar un mapa se procede de la siguiente manera. Observa que en el directorio de trabajo se ha creado un nuevo archivo png.

p <- map + 
  north(zona_estudio.f, symbol=10)+
  scalebar(zona_estudio.f, dist = 100,st.size=4) 

ggsave(p, file = "map1.png", width = 6, height = 4.5, type = "cairo-png")

C- Mapas de Google

La libreria ggplot también puede combinarse con la libreria ggmap para obtener mapas base de google.Las siguientes funciones son de ggmap.

map1 <- get_map(location = c(lon = -79, lat = -3), zoom = 11, maptype = 'terrain')
## Source : https://maps.googleapis.com/maps/api/staticmap?center=-3,-79&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN
ggmap(map1)

map2 <- get_map(location = c(lon = -79, lat = -3), zoom = 9, maptype = 'satellite')
## Source : https://maps.googleapis.com/maps/api/staticmap?center=-3,-79&zoom=9&size=640x640&scale=2&maptype=satellite&language=en-EN
#Otros mapas (maptype) dispnibles: "terrain", "terrain-background", "satellite","roadmap", "hybrid", "toner", "watercolor", "terrain=labels", "terrain#-lines", "toner-2010", "toner-2011", "toner-background", "toner-hybrid", "toner-labels", "toner-lines", "toner-lite
ggmap(map2)

#La variable estaciones tiene dos campos llamados Long y Lat y los utilizaremos para los ejes x e y.
ggmap(map2)+ geom_point(data=estaciones@data, aes(x=Long, y=Lat, size=ALTURA), color="red")
## Warning: Using size for a discrete variable is not advised.
## Warning: Removed 9 rows containing missing values (geom_point).

#Dado que las capas de google están en wgs84, debemos transformar la capa de zona de estudio ese sistema de referencia y luego aplicar la función fortify.
zona_estudiowgs84 <- spTransform(zona_estudio,CRS("+init=epsg:4326"))#transformación
zona_estudiowgs84.f <- fortify(zona_estudiowgs84, region="id")
head(zona_estudiowgs84.f)
##        long       lat order  hole piece id group
## 1 -78.56954 -2.558652     1 FALSE     1  1   1.1
## 2 -78.56920 -2.558284     2 FALSE     1  1   1.1
## 3 -78.56868 -2.557736     3 FALSE     1  1   1.1
## 4 -78.56817 -2.557184     4 FALSE     1  1   1.1
## 5 -78.56811 -2.557143     5 FALSE     1  1   1.1
## 6 -78.56756 -2.556762     6 FALSE     1  1   1.1
map3 <- get_map(location = c(lon = mean(zona_estudiowgs84.f$long), lat = mean(zona_estudiowgs84.f$lat)), zoom = 7, maptype = 'satellite') #Si el zoom es mayor, esto hace que los polígonos queden cortados y la visualización produce un efecto de distorción de los mismos.
## Source : https://maps.googleapis.com/maps/api/staticmap?center=-3.437407,-79.304519&zoom=7&size=640x640&scale=2&maptype=satellite&language=en-EN
ggmap(map3)+
geom_polygon(data=zona_estudiowgs84.f, aes(x=long, y=lat, group = group), colour="white",fill="grey50", alpha=0.5 ) +
geom_point(data=estaciones@data, aes(x=Long, y=Lat), color="red")  + 
labs(x="Lat", y="Long", title="Provincias y estaciones")