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")
::coordinates(meuse) <- c("x", "y")
sp::plot(meuse)
sptitle("Puntos")
Líneas
Un objeto tipo líneas SpatialLines
se conforma de una secuencia de puntos en secuencia:
<- sp::coordinates(meuse)
cc <- sp::SpatialLines(list(sp::Lines(list(sp::Line(cc)), "line1")))
m.sl ::plot(m.sl)
sptitle("Líneas")
Polígonos
data(meuse.riv,package = "sp")
<- list(sp::Polygons(list(sp::Polygon(meuse.riv)),
meuse.lst "meuse.riv"))
<- sp::SpatialPolygons(meuse.lst)
meuse.pol ::plot(meuse.pol, col = "grey")
sptitle("Polígonos")
Pixels
Convertimos el
data(meuse.grid,package = "sp")
::coordinates(meuse.grid) <- c("x", "y")
sp<- as(meuse.grid, "SpatialPixels")
meuse.grid 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")
::plot(meuse.pol, col = "grey", add = TRUE)
sp::plot(meuse, add = TRUE) sp
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))
::plot(meuse.pol, axes = TRUE)
sp::plot(meuse.pol, axes = FALSE)
spaxis(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.
::plot(meuse)
sp::plot(meuse.pol, add = TRUE)
sp::plot(meuse)
sp::SpatialPolygonsRescale(sp::layout.scale.bar(),offset = locator(1),
spscale = 1000, fill = c("transparent", "black"), plot.grid = FALSE)
text(locator(1), "0")
text(locator(1), "1 km")
::SpatialPolygonsRescale(layout.north.arrow(), offset = locator(1),
spscale = 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')`
<- maps::map("world", interior = FALSE, xlim = c(-179, 179), ylim = c(-89, 89), plot = FALSE)
wrld
# `pruneMap` redefine los límites del mapa, El objeto de entrada debe ser creado con la librería `maps`.
<- maptools::pruneMap(wrld, xlim = c(-179, 179))
wrld_p
# 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`
<- sp::CRS("+proj=longlat +ellps=WGS84")
llCRS <- maptools::map2SpatialLines(wrld_p, proj4string = llCRS)
wrld_sp <- sp::CRS("+proj=moll")
prj_new
# Con `spTransform` estamos proyectando el sistema de referencia de `wrld_sp` a `prj_new` y se guarda en `wrld_proj`
# library(rgdal)
<- sp::spTransform(wrld_sp, prj_new)
wrld_proj # Con `gridlines` creamos la cuadrícula de referencia:
<- sp::gridlines(wrld_sp, easts = c(-179, seq(-150,
wrld_grd 150, 50), 179.5), norths = seq(-75, 75, 15), ndiscr = 100)
# Proyectamos la cuadrícula de referencia:
<- sp::spTransform(wrld_grd, prj_new)
wrld_grd_proj
# `gridat` funciona como `gridlines` pero es para puntos específicos
<- sp::gridat(wrld_sp, easts = 0, norths = seq(-75,
at_sp 75, 15), offset = 0.3)
<- sp::spTransform(at_sp, prj_new)
at_proj ::plot(wrld_proj, col = "grey60")
sp::plot(wrld_grd_proj, add = TRUE, lty = 3, col = "grey70")
sptext(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.
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)
= gstat::idw(zinc~1, meuse, meuse.grid) zn.idw
## [inverse distance weighted interpolation]
# Atributos y parámetros
= gray.colors(4, 0.55, 0.95)
grays # grays = palette(hcl.colors(4, "viridis"))
image(zn.idw, col = grays)
::plot(meuse.pol, add = TRUE)
sp::plot(meuse, pch = 1, cex = sqrt(meuse$zinc)/20, add = TRUE)
sp<- c(100, 200, 500, 1000, 2000)
legVals 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
= gstat::idw(log(zinc)~1, meuse,meuse.grid)@data$var1.pred ln.zn.idw
## [inverse distance weighted interpolation]
<- zn.idw
zn names(zn) <- c("direct", "log")
@data$log <- ln.zn.idw
znsummary(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
::spplot(zn[c("direct", "log")]) sp
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`
<- as(meuse.grid, "SpatialPixelsDataFrame")
meuse.grid # Creo una imagen:
<- as.image.SpatialGridDataFrame(meuse.grid["dist"])
im # Se crea un `SpatialLinesDataFrame`
<- maptools::ContourLines2SLDF(contourLines(im))
cl spplot(cl)
4.3.2 Agregando referencias en 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 ) |
<- list("sp.polygons", meuse.pol)
river <- list("SpatialPolygonsRescale", layout.north.arrow(),
north offset = c(178750, 332500), scale = 400)
<- list("SpatialPolygonsRescale", layout.scale.bar(),
scale offset = c(180200, 329800), scale = 1000, fill = c("transparent",
"black"))
<- list("sp.text", c(180200, 329950), "0")
txt1 <- list("sp.text", c(181200, 329950), "1 km")
txt2 <- list("sp.points", meuse, pch = 3, col = "black")
pts <- list(river, north, scale, txt1, txt2,pts)
meuse.layout ::spplot(zn["log"], sp.layout = meuse.layout) sp
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.
<- function(uu)
read_git_shp
{#creamos un par de archivos temporales
<- tempfile()
temp <- tempfile()
temp2 #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
<- list.files(temp2, pattern = ".shp$",full.names=TRUE)
your_SHP_file
= strsplit(your_SHP_file,"/")
ff = unlist(ff)
ff = ff[length(ff)]
ff = strsplit(ff,".shp")
ff = unlist(ff)
ff
= rgdal::readOGR(your_SHP_file,layer = ff)
datos unlink(temp)
unlink(temp2)
return(datos)
}
library(rgdal) #readOGR
library(RColorBrewer)
library(classInt)
library(raster)
### Leo el shape:
<- "https://github.com/vmoprojs/DataLectures/raw/master/SpatialData/shape_Ambato_urbano.zip"
uu <- read_git_shp(uu) poligonos
## 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
<- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/EmergenciasTransito.txt"
uu <- read.csv(uu, sep ="")
datos <- subset(datos,Month == 12) datos
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:
<- cbind(datos$Longitude,datos$Latitude)
coord # Convierto puntos en "Puntos Espaciales" para que R los pueda interpretar:
= SpatialPoints(coord, proj4string=CRS("+proj=longlat +datum=WGS84"))
coord 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
:
<- spTransform(coord, CRS("+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs")) coord
Veamos los resultados gráficamente:
# Grafico el mapa:
::plot(poligonos)
sppoints(coord, col = "blue", pch=13, cex = 0.1)
¿Si uso proyecciones diferentes?
<- cbind(datos$Longitude,datos$Latitude)
coord1 # Convierto puntos en "Puntos Espaciales" para que R los pueda interpretar:
= SpatialPoints(coord, proj4string=CRS("+proj=longlat +datum=WGS84"))
coord1 # Grafico el mapa:
::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) sp