4 Visualización
Tomado de Pebesma and Bivand (2023), capítulo 3 & 4
Librerías usadas:
sf
La siguiente función nos sirve para leer datos comprimidos desde github:
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 = sf::st_read(your_SHP_file,layer = ff)
unlink(temp)
unlink(temp2)
return(datos)
}4.1 Métodos de visualización para objetos sf y sfc
4.1.1 Solo la Geometría: sfc
Las columnas-lista de geometría (objetos de la clase sfc, obtenidos mediante el método st_geometry) solo muestran la geometría:
library(sf)
uu <- "https://github.com/vmoprojs/DataLectures/raw/master/SpatialData/shape_Ambato_urbano.zip"
nc <- read_git_shp(uu)## Reading layer `PARROQUIA_URBANA' from data source
## `/private/var/folders/0p/n_r_hl095sv7nktfp_8n9n_80000gn/T/Rtmp4EHXgX/file2ca3b282cf3/PARROQUIA_URBANA.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 9 features and 5 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 761357 ymin: 9855965 xmax: 770602.4 ymax: 9865474
## Projected CRS: WGS 84 / UTM zone 17S
que se pueden anotar con colores, símbolos, etc., como los gráficos base habituales. Por ejemplo, los puntos se agregan a un polígono mediante:
plot(st_geometry(nc), col = sf.colors(nrow(nc), categorical = TRUE), border = 'grey',
axes = TRUE)
plot(st_geometry(st_centroid(nc)), pch = 3, col = 'red', add = TRUE)## Warning: st_centroid assumes attributes are constant over geometries
y posteriormente se pueden agregar leyendas, títulos, etc. border = NA elimina los bordes del polígono.
Como puede verse, los ejes trazados son sensibles al CRS, y en el caso de las coordenadas de longitud/latitud, los símbolos de grados y la orientación se agregan si axes = TRUE.
4.2 Geometría con atributos: sf
El gráfico predeterminado de un objeto sf es un gráfico múltiple de todos los atributos, hasta un máximo razonable:
con una advertencia cuando no todos los atributos se pueden graficar razonablemente. Se puede aumentar el número máximo de mapas que se trazarán mediante
4.3 Color, lugar y tamaño
En caso de que se seleccione un solo atributo, de forma predeterminada se proporciona una clave de color; para nc esto está a continuación:
nc$Shape_Area <- nc$Shape_Area/1000 # convertir a km
nc$Shape_Leng <- nc$Shape_Leng/1000
plot(nc["Shape_Area"])
La ubicación de la clave de color puede controlarse (1=abajo, 2=izquierda, 3=arriba y 4=derecha):
4.4 Intervalos de clase
Los saltos de color (intervalos de clase) se pueden controlar mediante los argumentos breaks y “nbreaks. nbreaks especifica el número de breaks; breaks es un vector con valores:
## 0% 25% 50% 75% 100%
## 525.4482 3281.1907 3819.9076 4256.6631 5416.2222
Otras formas de marcar las clases es con el argumento style de classInt::classIntervals().
Otros métodos incluyen
pretty(predeterminado), da como resultado saltos de clase redondeados y tiene como efecto secundario quenbreakspuede respetarse sólo de forma aproximada.equalpara dividir el rango de datos ennbreaksclases iguales,quantilepara usar cuantiles como divisiones de clases, y
4.5 Con otros paquetes
4.5.1 ggplot
Contiene una geom especialmente para objetos de características simples. Actualmente, el soporte es bueno para los polígonos; para líneas o puntos.
Se puede colorear los polígonos con aes
library(dplyr)
library(tidyr)
nc2 <- nc %>% dplyr::select(Shape_Area, Shape_Leng, geometry) %>% gather(VAR, SID, -geometry)
ggplot() +
geom_sf(data = nc2, aes(fill = SID)) +
facet_wrap(~VAR, ncol = 1)
4.5.2 mapview
El paquete mapview crea mapas interactivos en páginas html, utilizando el paquete leaflet como caballo de batalla. Se encuentran ejemplos extensos aquí.
Un ejemplo se obtiene mediante
4.6 Transformación
Importamos datos de puntos en coordenadas long-lat
uu <- "https://raw.githubusercontent.com/vmoprojs/DataLectures/master/SpatialData/EmergenciasTransito.txt"
puntos <- read.csv(uu,sep = "\t")
puntos <- subset(puntos,Month == 12)Veamos los elementos de la proyeccion del polígono:
## Coordinate Reference System:
## User input: WGS 84 / UTM zone 17S
## wkt:
## 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.6.1 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
proj4es+ 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 <- data.frame(x=puntos$Longitude,y=puntos$Latitude)
# Convierto puntos en "Puntos Espaciales" para que R los pueda interpretar:
coord = sf::st_as_sf(coord,coords = c("x","y"))
sf::st_crs(coord) <- "+proj=longlat +datum=WGS84"Para transformar la proyección o el datum de un mapa se usa la función spTransform del paquete sp:
# Transformo los datos de tal forma que tengan la misma proyección que el shape:
coord <- sf::st_transform(coord, crs=sf::st_crs(nc))Veamos los resultados gráficamente:
¿Si uso proyecciones diferentes?
coord1 <- data.frame(x=puntos$Longitude,y=puntos$Latitude)
# Convierto puntos en "Puntos Espaciales" para que R los pueda interpretar:
coord1 = sf::st_as_sf(coord1,coords = c("x","y"))
sf::st_crs(coord1) <- "+proj=longlat +datum=WGS84"
# Grafico el mapa:
plot(coord1, col = "blue", pch=13, cex = 0.1)
plot(st_geometry(nc),add = TRUE)
## Warning in plot.sf(nc, add = TRUE): ignoring all but the first attribute