En esta lección aprenderá a visualizar archivos de geoinformación con diferentes librerias como son sp y ggplot. Fuentes de información para esta práctica son: http://rspatial.r-forge.r-project.org/gallery/ y http://www.kevjohnson.org/making-maps-in-r/

Primero cargaremos y prepararemos los datos a utilizar.

#Librerias
library(sp)
library(rgdal) #leer archivos de geoinformación
## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
## Path to GDAL shared files: C:/Users/DOCENTE11/Documents/R/win-library/3.0/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: C:/Users/DOCENTE11/Documents/R/win-library/3.0/rgdal/proj
library(rgeos) #geoprocesos vectoriales
## Warning: package 'rgeos' was built under R version 3.0.3
## rgeos version: 0.3-6, (SVN revision 450)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Polygon checking: TRUE
library(sp) #visualización
library(ggplot2)#visualización
library(latticeExtra) #complemento visualización
## Warning: package 'latticeExtra' was built under R version 3.0.3
## Loading required package: RColorBrewer
## Loading required package: lattice
## 
## Attaching package: 'latticeExtra'
## 
## The following object is masked from 'package:ggplot2':
## 
##     layer
library(ggmap)#complemento visualización, mapas de google
## Warning: package 'ggmap' was built under R version 3.0.3
#Directorio de trabajo. Aqui deben estar todos los datos almacenados.
setwd("C:/curso_R/")

#Datos a utilizar.
shell.exec("https://drive.google.com/uc?export=download&id=0B0ea6usixQHFbFlBMnBEMXN5Yzg")
shell.exec("https://drive.google.com/uc?export=download&id=0B0ea6usixQHFbmNIWnZoaFk0c3M")#descargar datos

#Cargar archivos de geoinformación.
zona_estudio<-readOGR(".","prov_zona_estudio")
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "prov_zona_estudio"
## with 7 features and 11 fields
## Feature type: wkbPolygon with 2 dimensions
estaciones<-readOGR(".","estaciones_32717") #estaciones metereológicas
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "estaciones_32717"
## with 100 features and 14 fields
## Feature type: wkbPoint with 2 dimensions
rios<-readOGR(".","rio_torrente")#fuente IGM 50k
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "rio_torrente"
## with 2093 features and 5 fields
## Feature type: wkbLineString with 2 dimensions
#Cargar capa raster y procesarla.
elevacion <- readGDAL("elevation.asc")
## elevation.asc has GDAL driver GTiff 
## and has 798 rows and 844 columns
names(elevacion) <- "dem" # asignar nombre a campo elevacion
proj4string(elevacion) 
## [1] "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
str(elevacion@data)
## 'data.frame':    673512 obs. of  1 variable:
##  $ dem: int  0 0 0 0 0 0 0 0 0 0 ...
gridcell <- elevacion@grid@cellsize[1]
elev <- subset(elevacion, dem>0)

#Comprobar sistemas de referencias.
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=longlat +datum=WGS84 +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"
#Reproyectar capa de rios.
utm17 <- "+proj=utm +zone=17 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0" #utm zona 17 sur con datum WGS84
riosutm17 <- spTransform(rios,CRS(utm17))
proj4string(riosutm17)
## [1] "+proj=utm +zone=17 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0"
#Recortar capa de rios a la zona de estudio
rios_zonaestudio <- gIntersection(zona_estudio,riosutm17,  byid = TRUE) # esta función llevará un tiempo en ejecutarse, ten paciencia.
## Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_not_poly,
## "rgeos_intersection"): spgeom1 and spgeom2 have different proj4 strings

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. Será la variable “estaciones” y la columna “ALTURA”

names(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 comenzaremos a incluir elementos al mapa, como son recuadro con coordenadas, leyenda dentro del mapa, leyenda personalizada y cambio de color a los elementos de la leyenda.

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

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

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

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

#Cambio de color de leyenda. Otras alternativas a heat.colors() son rainbow() y terrain.colors().
spplot(estaciones, "ALTURA",key.space=list(x=0.7,y=0.3,corner=c(0,1)),
  scales=list(draw=T), cuts = 3,col.regions=heat.colors(3),   legendEntries = c("low", "intermediate", "high"))# 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 libreria latticeExtra. Deberas instalar esta libreria.

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

l2 = list("sp.lines", rios_zonaestudio,  col="blue")
spplot(zona_estudio, "DPA_PROVIN", col="grey10", fill="grey90", sp.layout = list(l2,l1))

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

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

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

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

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

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

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

Otros elementos del mapa a incorporar son la escala gráfica, la flecha indicando el norte y titulo del mapa.

#Escala gráfica.
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
  offset = c(600000,9700000), scale = 50000, fill=c("transparent","black"))
  text1 = list("sp.text", c(600000,9709000), "0")
  text2 = list("sp.text", c(650000,9709000), "50 km")
#Flecha de norte.
arrow = list("SpatialPolygonsRescale", layout.north.arrow(), 
    offset = c(620000,9715000), scale = 30000)

spplot(estaciones, "ALTURA", sp.layout=list(scale,text1,text2,arrow),scales=list(draw=T),   main = "Estaciones")

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.

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

ggplot

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

#En el caso de polígonos y lineas es necesario aplicar 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")
#Se cambia el nombre del campo long y lat por x e y, para facilitar la nomenclatura de los ejes en el gráfico a crear.
names(zona_estudio.f) <- c("x","y", "order","hole","piece","group","id" )
head(zona_estudio.f)
##          x       y 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=x, y=y, group = group),colour="white", fill="grey50" )

# Superponer capa
ggplot() + 
  geom_polygon(data=zona_estudio.f, aes(x=x, y=y, group = group),colour="white", fill="grey50" )+
  geom_point(data=estaciones@data, aes(x=x, y=y), color="red")

# Transparencia, etiquetas de ejes y título del mapa
ggplot() + 
geom_polygon(data=zona_estudio.f, aes(x=x, y=y, group = group),colour="white", fill="grey50" , 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")

# Transparencia, etiquetas de ejes y título del mapa
ggplot() + 
geom_polygon(data=zona_estudio.f, aes(x=x, y=y, group = group),colour="white", fill="grey50" , 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")

Para guardar un mapa, se procede de la siguiente manera

p <- ggplot() + geom_polygon(data=zona_estudio.f, aes(x=x, y=y, group = group),colour="white", fill="grey50" , 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")

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

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

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 27 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.
wgs84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 
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 = 8, maptype = 'satellite') #Si el zoom es mayor, esto hace que los poligonos queden cortados y la visualización produce un efecto extraño distorcionando los mismos.
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=-3.437407,-79.304519&zoom=8&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 ) +
labs(x="Lat", y="Long", title="Provincias y estaciones") + 
geom_point(data=estaciones@data, aes(x=Long, y=Lat), color="red")

Ha llegado al final de esta lección!