4 Visualización

Tomado de R. S. Bivand et al. (2013), capítulo 3 & 4

Librerías usadas:

  • sp
  • maps
  • maptools
  • raster
  • rgdal
  • gstat

4.1 Introducción

El paquete sp tiene el método plot que usa el sistema gráfico tradicional de R (plot, image, lines, points, etc.) así como el método ssplot que usa un sistema parecido al de lattice (Sarkar (2008)).

4.2 El sistema tradicional

4.2.1 Puntos, líneas, polígonos y cuadrículas.

Puntos

Notemos que meuse es un objeto data.frame creado al usar coordinates y lo convierte en un SpatialPointsDataFrame.

# library(sp)
data(meuse,package = "sp")
sp::coordinates(meuse) <- c("x", "y")
sp::plot(meuse)
title("Puntos")

Líneas

Un objeto tipo líneas SpatialLines se conforma de una secuencia de puntos en secuencia:

cc <- sp::coordinates(meuse)
m.sl <- sp::SpatialLines(list(sp::Lines(list(sp::Line(cc)), "line1")))
sp::plot(m.sl)
title("Líneas")

Polígonos

data(meuse.riv,package = "sp")
meuse.lst <- list(sp::Polygons(list(sp::Polygon(meuse.riv)),
"meuse.riv"))
meuse.pol <- sp::SpatialPolygons(meuse.lst)
sp::plot(meuse.pol, col = "grey")
title("Polígonos")

Pixels

Convertimos el

data(meuse.grid,package = "sp")
sp::coordinates(meuse.grid) <- c("x", "y")
meuse.grid <- as(meuse.grid, "SpatialPixels")
image(meuse.grid, col = "grey")
title("Cuadrícula")

Elementos

Podemos agregar elementos al mapa usando el parámetro add = TRUE:

image(meuse.grid, col = "lightgrey")
sp::plot(meuse.pol, col = "grey", add = TRUE)
sp::plot(meuse, add = TRUE)

Hasta ahora solo hemos mostrado la geometría (topología, formas) que trabaja bien cuando no se tiene un CRS o es NA, ahora veamos los tributos medidos en esos objetos.

4.2.2 Ejes & Diseño

El argumento lógico axes se puede configurar para los ejes:

layout(matrix(c(1, 2), 1, 2))
sp::plot(meuse.pol, axes = TRUE)
sp::plot(meuse.pol, axes = FALSE)
axis(1, at = c(178000 + 0:2 * 2000), cex.axis = 0.7)
axis(2, at = c(326000 + 0:3 * 4000), cex.axis = 0.7)
box() # completa la caja del segundo grafico

La modificación de los márgenes mar en el comando par, por ejemplo en par(mar=c (3,3,2,1)) optimiza aún más el uso del espacio cuando se dibujan los ejes, dejando (poco) espacio para un título.

sp::plot(meuse)
sp::plot(meuse.pol, add = TRUE)
sp::plot(meuse)
sp::SpatialPolygonsRescale(sp::layout.scale.bar(),offset = locator(1),
scale = 1000, fill = c("transparent", "black"), plot.grid = FALSE)
text(locator(1), "0")
text(locator(1), "1 km")
sp::SpatialPolygonsRescale(layout.north.arrow(), offset = locator(1),
scale = 400, plot.grid = FALSE)

4.2.3 Grados en ejes y cuadrícula de referencia

Datos no proyectados suelen estar en latitud y longtitud expresadas en grados, con grados negativos al oeste del primer meridiano y al sur del Ecuador.

Se puede usar degAxis para controlar las marcas en los ejes.

Si se necesita agregar una cuadrícula de referencia, la función gridLines puede usarse para generar un objeto SpatialLines.

Las líneas de cuadrícula pueden ser cuadrículas de latitud/longitud, y estas no son líneas rectas. Esto se logra generando una cuadrícula para datos no proyectados, proyectándolos y trazándolos sobre el mapa.

# library(maptools)
# library(maps)

# La función `map` dibuja lineas y poligonos de una base de datos (principalmente de USA), mas info en `help(package='maps')`
wrld <- maps::map("world", interior = FALSE, xlim = c(-179, 179), ylim = c(-89, 89), plot = FALSE)

# `pruneMap` redefine los límites del mapa, El objeto de entrada debe ser creado con la librería `maps`.
wrld_p <- maptools::pruneMap(wrld, xlim = c(-179, 179))

# Notemos que `wrld_sp` no tiene un CRS (hasta ahora es solo una lista), con la función `map2SpatialLines` se transforma en un objeto `SpatialLines`
llCRS <- sp::CRS("+proj=longlat +ellps=WGS84")
wrld_sp <- maptools::map2SpatialLines(wrld_p, proj4string = llCRS)
prj_new <- sp::CRS("+proj=moll")

# Con `spTransform` estamos proyectando el sistema de referencia de `wrld_sp` a `prj_new` y se guarda en `wrld_proj`
# library(rgdal)
wrld_proj <- sp::spTransform(wrld_sp, prj_new)
# Con `gridlines` creamos la cuadrícula de referencia:
wrld_grd <- sp::gridlines(wrld_sp, easts = c(-179, seq(-150,
150, 50), 179.5), norths = seq(-75, 75, 15), ndiscr = 100)
# Proyectamos la cuadrícula de referencia:
wrld_grd_proj <- sp::spTransform(wrld_grd, prj_new)

# `gridat` funciona como `gridlines` pero es para puntos específicos
at_sp <- sp::gridat(wrld_sp, easts = 0, norths = seq(-75,
75, 15), offset = 0.3)
at_proj <- sp::spTransform(at_sp, prj_new)
sp::plot(wrld_proj, col = "grey60")
sp::plot(wrld_grd_proj, add = TRUE, lty = 3, col = "grey70")
text(sp::coordinates(at_proj), pos = at_proj$pos, offset = at_proj$offset,
labels = parse(text = as.character(at_proj$labels)),
cex = 0.6)

4.2.4 Atributos y legenda

La escala del mapa es la relación entre la longitud de una unidad en el mapa y una unidad en el mundo real.

Si queremos mostrar características o atributos de los objetos, necesitamos usar el tipo, tamaño o color de los símbolos, líneas o polígonos.

Parámetros útiles usados en plot & image. * Después de usar image.
Clase Parámetro Significado Más ayuda
SpatialLinesDataFrame col Color ?lines
lwd Color ?lines
lty Color ?lines
SpatialPolygonsDataFrame border Color del borde ?polygon
density Densidad de líneas de relleno ?polygon
angle Ángulo de líneas de relleno ?polygon
lty tipo de línea ?polygon
pbg Color del polígono ?polygon
SpatialPointsDataFrame pch Símbolo ?points
col Color ?points
bg Color de relleno ?points
cex Tamaño del símbolo ?points
SpatialPixelsDataFrame* & SpatialGridDataFrame zlim Límites del atributo ?image.default
col Colores ?image.default
breaks Puntos de separación ?image.default
# `idw` Interpolación ponderada por la inversa de la distancia, devuelve un objeto `SpatialPixelsDataFrame`
# inverse distance weightings
require(gstat)
zn.idw = gstat::idw(zinc~1, meuse, meuse.grid)
## [inverse distance weighted interpolation]
# Atributos y parámetros
grays = gray.colors(4, 0.55, 0.95)
# grays = palette(hcl.colors(4, "viridis"))
image(zn.idw, col = grays)
sp::plot(meuse.pol, add = TRUE)
sp::plot(meuse, pch = 1, cex = sqrt(meuse$zinc)/20, add = TRUE)
legVals <- c(100, 200, 500, 1000, 2000)
legend("left", legend = legVals, pch = 1, pt.cex = sqrt(legVals)/20,
bty = "n", title = "Observado")
legend("topleft", legend = c("100-200", "200-400", "400-800",
"800-1800"), fill = grays, bty = "n", title = "Interpolación")

4.3 ssplot

Las gráficas con ssplot son un poco más difíciles de manejar inicialmente porque la anotación de la gráfica, la adición de información como leyendas, líneas, texto, etc., se maneja de manera diferente y debe pensarse primero.

La ventaja que ofrecen es que muchos mapas se pueden componer en (conjuntos de) gráficos únicos, fácil y eficientemente.

Ejemplo

ln.zn.idw = gstat::idw(log(zinc)~1, meuse,meuse.grid)@data$var1.pred
## [inverse distance weighted interpolation]
zn <- zn.idw
names(zn) <- c("direct", "log")
zn@data$log <- ln.zn.idw
summary(zn@data)
##      direct            log       
##  Min.   : 128.4   Min.   :4.791  
##  1st Qu.: 293.2   1st Qu.:5.484  
##  Median : 371.4   Median :5.694  
##  Mean   : 423.2   Mean   :5.777  
##  3rd Qu.: 499.8   3rd Qu.:6.041  
##  Max.   :1805.8   Max.   :7.482
sp::spplot(zn[c("direct", "log")])

4.3.1 Líneas, polígonos y cruadrículas

El primer argumento de ssplot es un objeto espacial de puntos, líneas, polígonos o cuadrícula.

El segundo atributo indica las columnas (o número de columnas) a graficar.

El resto de argumentos controlan: color, símbolos, leyendas, tamaño, etc.

# library(maptools)
# Importo el data frame
rm(meuse.grid)
data(meuse.grid)
# Con coordinates se crea un objeto `SpatialPointsDataFrame`
coordinates(meuse.grid) <- c("x", "y")
# Lo convierto en un `SpatialPixelsDataFrame`
meuse.grid <- as(meuse.grid, "SpatialPixelsDataFrame")
# Creo una imagen:
im <- as.image.SpatialGridDataFrame(meuse.grid["dist"])
# Se crea un `SpatialLinesDataFrame`
cl <- maptools::ContourLines2SLDF(contourLines(im))
spplot(cl)

4.3.2 Agregando referencias en la salida gráfica

Funciones para la salida gráfica
Función Clase Argumentos
sp.points SpatialPoints pch, cex, col
sp.polygons SpatialPolygons lty, lwd, col
sp.lines SpatialLines lty, lwd, col
sp.text text (ver panel.text)
river <- list("sp.polygons", meuse.pol)
north <- list("SpatialPolygonsRescale", layout.north.arrow(),
offset = c(178750, 332500), scale = 400)
scale <- list("SpatialPolygonsRescale", layout.scale.bar(),
offset = c(180200, 329800), scale = 1000, fill = c("transparent",
"black"))
txt1 <- list("sp.text", c(180200, 329950), "0")
txt2 <- list("sp.text", c(181200, 329950), "1 km")
pts <- list("sp.points", meuse, pch = 3, col = "black")
meuse.layout <- list(river, north, scale, txt1, txt2,pts)
sp::spplot(zn["log"], sp.layout = meuse.layout)

4.4 Proyecciones & Transformación

R nos proporciona el CRS de sus datos en formato proj4. Analicemos por partes el proj4. Contiene todos los elementos CRS individuales que R u otro GIS pueden necesitar. Cada elemento se especifica con un signo +, similar a cómo un archivo .csv está delimitado o dividido por un ,. Después de cada +, verá que se define el elemento CRS. Por ejemplo, proyección (proj =) y datum (datum =).

4.4.1 proj4 en UTM

La proyección UTM se especifica de la siguiente manera:

+ proj = utm + zone = 18 + datum = WGS84 + unidades = m + no_defs + ellps = WGS84 + towgs84 = 0,0,0

  • proj = utm: la proyección es UTM, UTM tiene varias zonas

  • zone = 18: la zona es 18

  • datum = WGS84: el datum WGS84 (el datum se refiere a la referencia 0,0 para el sistema de coordenadas utilizado en la proyección)

  • unidades = m: las unidades de las coordenadas están en METROS

  • ellps = WGS84: el elipsoide (cómo se calcula la redondez de la Tierra) para los datos es WGS84

Toma en cuenta que zone es exclusiva de la proyección UTM. No todos los CRS tendrán una zona.

read_git_shp <- function(uu)
{
  #creamos un par de archivos temporales
  temp <- tempfile()
  temp2 <- tempfile()
  #decargamos el zip folder y lo guardamos en 'temp' 
  
  download.file(uu,temp)
  #descomprimir en 'temp' y guardarlo en 'temp2'
  unzip(zipfile = temp, exdir = temp2)
  #encontramos los archivos SHP
  #el $ al final de ".shp$" asegura que no encontremos archivos del tipo .shp.xml 
  your_SHP_file <- list.files(temp2, pattern = ".shp$",full.names=TRUE)
  
  ff = strsplit(your_SHP_file,"/")
  ff = unlist(ff)
  ff = ff[length(ff)]
  ff = strsplit(ff,".shp")
  ff = unlist(ff)
  
  datos = rgdal::readOGR(your_SHP_file,layer = ff)
  unlink(temp)
  unlink(temp2)
  return(datos)
}
library(rgdal) #readOGR
library(RColorBrewer)
library(classInt)
library(raster)


### Leo el shape:
uu <- "https://github.com/vmoprojs/DataLectures/raw/master/SpatialData/shape_Ambato_urbano.zip"
poligonos <- read_git_shp(uu)
## OGR data source with driver: ESRI Shapefile 
## Source: "/private/var/folders/0p/n_r_hl095sv7nktfp_8n9n_80000gn/T/RtmpZgpGxU/file1556340d0f298/PARROQUIA_URBANA.shp", layer: "PARROQUIA_URBANA"
## with 9 features
## It has 5 fields
uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/EmergenciasTransito.txt"
datos <- read.csv(uu, sep ="")
datos <- subset(datos,Month == 12)

Veamos los elementos de proj.4

crs(poligonos)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 17S",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 17S",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-81,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",10000000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     ID["EPSG",32717]]

4.4.2 proj4 geográfico (lat/long)

La proyección lat/lon se especifica de la siguiente manera:

+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

  •   **proj = longlat**: los datos están en un sistema de coordenadas geográficas (latitud y longitud)
  •   **datum = WGS84**: el datum WGS84 (el datum se refiere a la referencia 0,0 para el sistema de coordenadas utilizado en la proyección)
  •   **ellps = WGS84**: el elipsoide (cómo se calcula la redondez de la tierra) es WGS84

Ten en cuenta que no hay unidades especificadas arriba. Esto se debe a que este sistema de referencia de coordenadas geográficas se encuentra en latitud y longitud, que con mayor frecuencia se registra en grados decimales.

¡Pilas!: la última parte de cada cadena proj4 es + towgs84 = 0,0,0. Este es un factor de conversión que se utiliza si se requiere una conversión de datum.

# Coord: tiene los datos de longitud y latitud:
coord <- cbind(datos$Longitude,datos$Latitude)
# Convierto puntos en "Puntos Espaciales" para que R los pueda interpretar:
coord = SpatialPoints(coord, proj4string=CRS("+proj=longlat +datum=WGS84"))
crs(coord)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

4.4.3 Transformación

Para transformar la proyección o el datum de un mapa se usa la función spTransform del paquete sp:

coord <- spTransform(coord, CRS("+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs"))

Veamos los resultados gráficamente:

# Grafico el mapa:
sp::plot(poligonos)
points(coord, col = "blue", pch=13, cex = 0.1)

¿Si uso proyecciones diferentes?

coord1 <- cbind(datos$Longitude,datos$Latitude)
# Convierto puntos en "Puntos Espaciales" para que R los pueda interpretar:
coord1 = SpatialPoints(coord, proj4string=CRS("+proj=longlat +datum=WGS84"))
# Grafico el mapa:
sp::plot(coord, col = "blue", pch=13, cex = 0.1)
sp::plot(poligonos,add = TRUE)

sp::plot(coord1, col = "red", pch=13, cex = 0.1)
sp::plot(poligonos,add = TRUE)