Visualización de Datos Georreferenciados

R Markdown

Sobre R Markdown, puede ser consultada Aqui

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(latticeExtra) #complemento visualización
library(ggmap)#complemento visualización, mapas de google
library(ggsn) #escala gráfica y norte para ggplot
library(raster)  # grids, rasters
library(rasterVis)  # raster visualisation
library(maptools)
library(dismo)

#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")

Visualización del Guatemala

g = gmap("Guatemala")  # elegir cualquier país

## Loading required namespace: XML

plot(g, inter=TRUE)

Localización específica

gs = gmap('Antigua Guatemala, Guatemala', type='terrain')
plot(gs, inter=TRUE)

Ha llegado al final de este material!

Visualización de Datos Georreferenciados

Olinda, R. A.

Octubre, 2016

Universidad de San Carlos de Guatemala

Centro Universitario del Norte

USAC

R Markdown

Temario

A- Librerías y datos

B- Visualización con sp