Maestría en Ecohidrologí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

Curso completo en: http://rpubs.com/daniballari/analisis_espacial_indice



En este material aprenderá a visualizar archivos de geoinformación con diferentes librerias como son sp, 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 sp
C- Visualización con ggplot
D- Utiliza mapas de Google

A- Librerías y datos

Comienza cargando las librerías y preparando los datos a utilizar.

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("C:/R_ecohidrologia/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
estaciones<-readOGR(".","estaciones_32717") #estaciones metereológicas
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "estaciones_32717"
## with 40 features
## It has 14 fields
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
#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 sp

Primero comenzaremos utilizando la visualización de la libreria sp con la función spplot. Para ello, debemos elegir una variable y una columna a representar y que en este caso será la variable estaciones y la columna ALTURA

names(estaciones)#nombres de los campos de estaciones
##  [1] "C1"         "CODIGO"     "NOMBRE_DE"  "TIPO"       "ZONA_HIDRO"
##  [6] "Lat"        "Long"       "ALTURA"     "PROVINCIA"  "INSTITUCIO"
## [11] "FECHA_DE_I" "FECHA_DE_L" "x"          "y"
#Tres alternativas para indicar la columna a representar.

#spplot(estaciones, "ALTURA")
#spplot(estaciones["ALTURA"])
spplot(estaciones[8])

Ahora incluiremos elementos al mapa, como son grilla con coordenadas, leyenda dentro del mapa, leyenda personalizada y cambio de color en la leyenda.

#Cuadricula con coordenadas.
spplot(estaciones, "ALTURA", scales=list(draw=T))

#Leyenda dentro del mapa
spplot(estaciones, "ALTURA",scales=list(draw=T),
       key.space=list(x=0.6,y=0.3,corner=c(0,1)))

#Leyenda personalizada.
spplot(estaciones, "ALTURA",scales=list(draw=T),
       key.space=list(x=0.6,y=0.3,corner=c(0,1)), 
       cuts = 3,
       legendEntries = c("low", "intermediate", "high"))

#Cambio de color de leyenda
spplot(estaciones, "ALTURA",scales=list(draw=T),
       key.space=list(x=0.6,y=0.3,corner=c(0,1)),
       cuts = 3,
       legendEntries = c("low", "intermediate", "high"),
       col.regions=c('blue', 'green', 'red'))

#Cambio de color de leyenda. Otras alternativas a heat.colors() son rainbow() y terrain.colors().
spplot(estaciones, "ALTURA",scales=list(draw=T),
       key.space=list(x=0.6,y=0.3,corner=c(0,1)),
       cuts = 3,
       legendEntries = c("low", "intermediate", "high"),
       col.regions=heat.colors(3)) # col.region = rainbow() or terrain.colors()

Hasta el momento hemos representado una sola capa. Para superponer mas de una capa, se presentan varias altarnativas. Una de ellas es con el argumento sp.layout y otra con la librería latticeExtra.

#Alternativa sp.layout
l1 = list("sp.points", estaciones, pch = 19, col="red")
spplot(zona_estudio, "DPA_PROVIN", sp.layout = list(l1))

#Distintos colores en R pueden ser observados aqui: http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf

l2 = list("sp.lines", rios,  col="cadetblue4")
spplot(zona_estudio, "DPA_PROVIN", col="grey20", fill="beige", sp.layout = list(l2,l1), colorkey=FALSE)

#Visualizar capa raster.
spplot(elevacion, "dem", col.regions=rev(gray(seq(0,1,.01))))

spplot(elevacion, "dem", col.regions=terrain.colors(100))

spplot(elevacion, "dem", col.regions=terrain.colors(100), sp.layout = list(l2,l1))

#Alternativa latticeExtra, funcion as.layer().
spplot(estaciones, "ALTURA",scales=list(draw=T),
       key.space=list(x=0.6,y=0.3,corner=c(0,1)),
       cuts = 3,
       legendEntries = c("low", "intermediate", "high"),
       col.regions=heat.colors(3)) +
       as.layer(spplot(zona_estudio[1], fill="grey90",  col="grey50"), under = TRUE)

#Para dibujar las líneas de los ríos, debemos seguir recurriendo al argumento sp.layout ya que spplot por si solo no permite dibujar lineas.
spplot(estaciones, "ALTURA",scales=list(draw=T),
       key.space=list(x=0.6,y=0.3,corner=c(0,1)),
       cuts = 3,
       legendEntries = c("low", "intermediate", "high"),
       col.regions=heat.colors(3)) +
       as.layer(spplot(zona_estudio[1], fill="grey90",  col="grey50", sp.layout = list(l2)), under = TRUE)

spplot(estaciones, "ALTURA",scales=list(draw=T),
       key.space=list(x=0.6,y=0.3,corner=c(0,1)),
       cuts = 3,
       legendEntries = c("low", "intermediate", "high"),
       col.regions=heat.colors(3)) +
       as.layer(spplot(zona_estudio[1], fill="transparent",  col="white", lwd=3, sp.layout = list(l2)), under = FALSE) + as.layer(spplot(elevacion, "dem", col.regions=terrain.colors(100)) , under = TRUE)

spplot(elevacion, "dem", col.regions=terrain.colors(100)) + as.layer(spplot(zona_estudio[1], fill="transparent", col="white", sp.layout = list(l1,l2)), under = FALSE)

También es posible dibujar varios paneles o gráficos al mismo tiempo, a partir de varias columnas contenida en el data.frame de la variable. Aqui los valores deben ser numéricos y estar representados en una misma escala. Es decir en los mapas de ejemplo no observaremos diferencia entre los puntos ya que la escala es muy dispar.

spplot(estaciones, c("Long", "Lat","ALTURA", "PROVINCIA"))

l3 = list("sp.polygons", zona_estudio,  col="grey10", fill="grey90")
spplot(estaciones, c("Long", "Lat","ALTURA", "PROVINCIA"),
  key.space = "right", 
    sp.layout=list(l3), 
    main = "Estaciones")

C- Visualización con ggplot

La libreria ggplot también 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 group id
## 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")

D- Utiliza 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')
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=-3,-79&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
ggmap(map1)

map2 <- get_map(location = c(lon = -79, lat = -3), zoom = 9, maptype = 'satellite')
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=-3,-79&zoom=9&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=false
#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: 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 group id
## 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 poligonos queden cortados y la visualización produce un efecto de distorción de los mismos.
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=-3.437407,-79.304519&zoom=7&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=false
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")

Ha llegado al final de este material!