library(sf)Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUElibrary(ggplot2)5 Importación de datos espaciales
Los datos geoespaciales en formato vectorial suelen almacenarse en formato shapefile. Dado que la estructura de puntos, líneas y polígonos es diferente, cada shapefile individual sólo puede contener un tipo de vector (todos los puntos, todas las líneas o todos los polígonos). No encontrarás una mezcla de objetos de puntos, líneas y polígonos en un mismo shapefile.
5.1 Leer un shapefile
Usaremos el paquete sf para trabajar con datos vectoriales en R.
Los shapefiles que vamos a importar son: - Puntos: Paraderos en vía Cuenca - Molleturo. - Lineas: Vía princial Cuenca - Molleturo y caminos pedestres hacia las lagunas en el Parque Nacional El Cajas. - Polígonos: Lagunas del Parque Nacional El Cajas.
Para leer el archivo shapefile en R usaremos el comando st_read().
puntos <- st_read("../data/shp/paraderos_cajas.shp")Reading layer `paraderos_cajas' from data source
`C:\Users\alex_ergostats\Documents\geo_stats_2024_nb\data\shp\paraderos_cajas.shp'
using driver `ESRI Shapefile'
Simple feature collection with 2 features and 2 fields
Geometry type: MULTIPOINT
Dimension: XY
Bounding box: xmin: -79.27525 ymin: -2.79308 xmax: -79.22255 ymax: -2.783234
Geodetic CRS: WGS 84lineas <- st_read("../data/shp/vias_cajas.shp")Reading layer `vias_cajas' from data source
`C:\Users\alex_ergostats\Documents\geo_stats_2024_nb\data\shp\vias_cajas.shp'
using driver `ESRI Shapefile'
Simple feature collection with 191 features and 26 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: -79.33892 ymin: -2.896237 xmax: -79.20061 ymax: -2.769916
Geodetic CRS: WGS 84poligonos <- st_read("../data/shp/lagunas_cajas.shp")Reading layer `lagunas_cajas' from data source
`C:\Users\alex_ergostats\Documents\geo_stats_2024_nb\data\shp\lagunas_cajas.shp'
using driver `ESRI Shapefile'
Simple feature collection with 220 features and 8 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -79.32098 ymin: -2.874182 xmax: -79.23114 ymax: -2.798271
Geodetic CRS: WGS 845.1.1 Metadatos y atributos del shapefile
Cuando importamos shapefile a R, la función st_read() guarda información acerca de los datos. Nos interesan especialmente los metadatos geoespaciales, que describen el formato, el SRC, la extensión y otros componentes de los datos vectoriales, así como los atributos que describen las propiedades asociadas a cada objeto vectorial individual.
5.1.2 Metadados espaciales
Los metadatos clave de todos los shapefiles incluyen:
- Tipo de objeto: la clase del objeto importado.
- Sistema de referencia de coordenadas (SRC): la proyección de los datos.
- Extensión: la extensión espacial (es decir, el área geográfica que cubre el shapefile) del shapefile.
Podemos ver los metadatos del shapefile simplemente llamando al objeto en la consola:
poligonosSimple feature collection with 220 features and 8 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -79.32098 ymin: -2.874182 xmax: -79.23114 ymax: -2.798271
Geodetic CRS: WGS 84
First 10 features:
full_id osm_id osm_type natural place water type
1 r11367469 11367469 relation water <NA> <NA> multipolygon
2 w76301259 76301259 way water <NA> <NA> <NA>
3 w296765494 296765494 way water <NA> lake <NA>
4 w297308408 297308408 way water <NA> lake <NA>
5 w297361636 297361636 way water <NA> <NA> <NA>
6 w297388787 297388787 way <NA> islet <NA> <NA>
7 w297389551 297389551 way water <NA> <NA> <NA>
8 w297389552 297389552 way water <NA> <NA> <NA>
9 w297389553 297389553 way water <NA> <NA> <NA>
10 w297391586 297391586 way water <NA> <NA> <NA>
name geometry
1 Luspa POLYGON ((-79.26245 -2.8013...
2 Hunanchi POLYGON ((-79.23757 -2.8721...
3 S/N POLYGON ((-79.23825 -2.8432...
4 Patos colorados POLYGON ((-79.2336 -2.83552...
5 Togllacocha POLYGON ((-79.25033 -2.7988...
6 Isla Luspa POLYGON ((-79.2597 -2.80711...
7 Canutillos POLYGON ((-79.25842 -2.8164...
8 Mangacocha POLYGON ((-79.25549 -2.8254...
9 Mangacocha POLYGON ((-79.25316 -2.8299...
10 Chachacocha POLYGON ((-79.25634 -2.8639...Nuestro objeto poligonos son polígonos de la clase POLYGON, en el sistema de coordenadas geográficas SRC WGS 84.
El SRC juega un papel crucial en la interpretación de los valores de extensión de un objeto espacial, pues define las unidades en las que se expresan. Esta información resulta fundamental para comprender la cobertura geográfica real de un shapefile u objeto espacial R, ya que estos representan la extensión espacial del objeto en el mundo real.
5.1.3 Atributos de los datos espaciales
Cada objeto de un shapefile tiene uno o más atributos asociados. Los atributos de un shapefile son similares a los campos o columnas de una hoja de cálculo. Cada fila de la hoja de cálculo tiene un conjunto de columnas asociadas que describen el elemento de la fila. En el caso de un shapefile, cada fila representa un objeto espacial - por ejemplo, una carretera, representada como una línea en un shapefile de líneas, tendrá asociada una fila de atributos. Estos atributos pueden incluir diferentes tipos de información que describen los objetos almacenados en un shapefile. Así, nuestra carretera, puede tener un nombre, longitud, número de carriles, límite de velocidad, tipo de carretera y otros atributos almacenados con ella.
5.2 Graficar un shapefile
A continuación, vamos a visualizar los datos en nuestro objeto sf utilizando el paquete ggplot2.
Vamos a personalizar nuestro trazado de límites estableciendo el tamaño, el color y el relleno de nuestro trazado. Al trazar objetos sf con ggplot2, es necesario utilizar el sistema de coordenadas coord_sf().
ggplot() +
geom_sf(data = poligonos) +
ggtitle("Parque Nacional El Cajas") +
coord_sf()
Para agregar colores a los polígonos usamos los argumentos fill para el relleno y color para el borde.
ggplot() +
geom_sf(data = poligonos, fill = "steelblue", color = "steelblue") +
ggtitle("Parque Nacional El Cajas") +
coord_sf()
5.2.1 Visualizar multiples shapefiles
Usando la librería ggplot2 esta tarea es bastante sencilla, simplemente tenemos que agregar las nuevas capas usando el símbolo +.
ggplot() +
geom_sf(data = poligonos, fill = "steelblue", color = "steelblue") +
geom_sf(data = lineas) +
geom_sf(data = puntos) +
ggtitle("Parque Nacional El Cajas") +
coord_sf()
Ahora que tenemos un gráfico básico, podemos personalizarlo cambiando la simbología. Como queremos mostrar los diferentes tipos de caminos en este sitio, también haremos que el color de los datos lines se establezca de acuerdo con el atributo highway. Esto también creará una leyenda para nuestro gráfico. Y finalmente, como también queremos incluir la ubicación de los paraderos en la leyenda, haremos que el relleno de los datos puntos se establezca de acuerdo con el atributo nombre.
ggplot() +
geom_sf(data = poligonos, fill = "steelblue", color = "steelblue") +
geom_sf(data = lineas, aes(color = highway)) +
geom_sf(data = puntos, aes(fill = nombre)) +
ggtitle("Parque Nacional El Cajas") +
coord_sf()Warning in vp$just: encuentros parciales de 'just' to 'justification'
A continuación, vamos a construir una leyenda personalizada utilizando la simbología (los colores y símbolos) que hemos utilizado para crear el gráfico anterior. Utilizaremos las funciones de ggplot2 para especificar la simbología y cambiar los nombres de las leyendas. Añadiremos color para las líneas y fill para el relleno de los puntos. También ajustaremos los títulos de las leyendas pasando un nombre a las respectivas paletas de color y relleno.
col_lineas <- c("orange", "brown")
col_puntos <- c("black", "black")
ggplot() +
geom_sf(data = poligonos, fill = "steelblue", color = "steelblue") +
geom_sf(data = lineas, aes(color = highway)) +
geom_sf(data = puntos, aes(fill = nombre), color = "black") +
scale_color_manual(values = col_lineas, name = "Tipos de camino") +
scale_fill_manual(values = col_puntos, name = "Paraderos") +
ggtitle("Parque Nacional El Cajas") +
coord_sf()Warning in vp$just: encuentros parciales de 'just' to 'justification'
5.2.2 Mapas interactivos
En R podemos crear mapas interactivos usando el paquete mapview.
library(mapview)The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
(status 2 uses the sf package in place of rgdal)mapview(poligonos)5.3 Leer un raster
Usaremos el paquete terra para trabajar con datos raster en R.
El raster que vamos a importar es:
- Un modelo de elevación digital (DEM) del volcán Chimborazo.
Para leer el archivo raster en R usaremos el comando rast().
library(terra)terra 1.7.39dem_chimborazo <- rast("../data/raster/DEM_Chimborazo.tif") 5.3.1 Metadatos del raster
En esta sección trabajaremos con una serie de archivos GeoTIFF. El formato GeoTIFF contiene un conjunto de etiquetas incrustadas con metadatos sobre los datos ráster. Podemos utilizar la función describe() para obtener información acerca del SRC, tamaño del raster, tamaño del pixel, entre otros.
describe("../data/raster/DEM_Chimborazo.tif") [1] "Driver: GTiff/GeoTIFF"
[2] "Files: ../data/raster/DEM_Chimborazo.tif"
[3] " ../data/raster/DEM_Chimborazo.tif.aux.xml"
[4] "Size is 294, 272"
[5] "Coordinate System is:"
[6] "GEOGCRS[\"WGS 84\","
[7] " ENSEMBLE[\"World Geodetic System 1984 ensemble\","
[8] " MEMBER[\"World Geodetic System 1984 (Transit)\"],"
[9] " MEMBER[\"World Geodetic System 1984 (G730)\"],"
[10] " MEMBER[\"World Geodetic System 1984 (G873)\"],"
[11] " MEMBER[\"World Geodetic System 1984 (G1150)\"],"
[12] " MEMBER[\"World Geodetic System 1984 (G1674)\"],"
[13] " MEMBER[\"World Geodetic System 1984 (G1762)\"],"
[14] " MEMBER[\"World Geodetic System 1984 (G2139)\"],"
[15] " ELLIPSOID[\"WGS 84\",6378137,298.257223563,"
[16] " LENGTHUNIT[\"metre\",1]],"
[17] " ENSEMBLEACCURACY[2.0]],"
[18] " PRIMEM[\"Greenwich\",0,"
[19] " ANGLEUNIT[\"degree\",0.0174532925199433]],"
[20] " CS[ellipsoidal,2],"
[21] " AXIS[\"geodetic latitude (Lat)\",north,"
[22] " ORDER[1],"
[23] " ANGLEUNIT[\"degree\",0.0174532925199433]],"
[24] " AXIS[\"geodetic longitude (Lon)\",east,"
[25] " ORDER[2],"
[26] " ANGLEUNIT[\"degree\",0.0174532925199433]],"
[27] " USAGE["
[28] " SCOPE[\"Horizontal component of 3D system.\"],"
[29] " AREA[\"World.\"],"
[30] " BBOX[-90,-180,90,180]],"
[31] " ID[\"EPSG\",4326]]"
[32] "Data axis to CRS axis mapping: 2,1"
[33] "Origin = (-78.857916666999998,-1.426527778000000)"
[34] "Pixel Size = (0.000277777778912,-0.000277777775735)"
[35] "Metadata:"
[36] " AREA_OR_POINT=Area"
[37] "Image Structure Metadata:"
[38] " INTERLEAVE=BAND"
[39] "Corner Coordinates:"
[40] "Upper Left ( -78.8579167, -1.4265278) ( 78d51'28.50\"W, 1d25'35.50\"S)"
[41] "Lower Left ( -78.8579167, -1.5020833) ( 78d51'28.50\"W, 1d30' 7.50\"S)"
[42] "Upper Right ( -78.7762500, -1.4265278) ( 78d46'34.50\"W, 1d25'35.50\"S)"
[43] "Lower Right ( -78.7762500, -1.5020833) ( 78d46'34.50\"W, 1d30' 7.50\"S)"
[44] "Center ( -78.8170833, -1.4643056) ( 78d49' 1.50\"W, 1d27'51.50\"S)"
[45] "Band 1 Block=294x6 Type=Float32, ColorInterp=Gray"
[46] " Min=4198.327 Max=6286.772 "
[47] " Minimum=4198.327, Maximum=6286.772, Mean=4884.034, StdDev=416.768"
[48] " NoData Value=-32768"
[49] " Metadata:"
[50] " STATISTICS_MAXIMUM=6286.7724609375"
[51] " STATISTICS_MEAN=4884.0342173796"
[52] " STATISTICS_MINIMUM=4198.3271484375"
[53] " STATISTICS_STDDEV=416.76784964225"
[54] " STATISTICS_VALID_PERCENT=100" Al igual que los metadatos del shapefile podemos obtener información del objeto raste simplemente llamando al objeto en la consola:
dem_chimborazoclass : SpatRaster
dimensions : 272, 294, 1 (nrow, ncol, nlyr)
resolution : 0.0002777778, 0.0002777778 (x, y)
extent : -78.85792, -78.77625, -1.502083, -1.426528 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : DEM_Chimborazo.tif
name : DEM_Chimborazo
min value : 4198.327
max value : 6286.772 Nuestro objeto dem_chimborazo es un raster de la clase SpatRaster, en el sistema de coordenadas geográficas SRC WGS 84 con unas dimensiones de dimensions : 272, 294, 1 (nrow, ncol, nlyr). Donde nlyr se refiere al número de bandas que contiene el raster.
La información anterior incluye un informe de valores mínimos y máximos, pero no otras estadísticas de rangos de datos. De forma similar a otras estructuras de datos de R, como vectores y dataframes, las estadísticas descriptivas de los datos ráster pueden obtenerse de la siguiente forma:
summary(values(dem_chimborazo)) DEM_Chimborazo
Min. :4198
1st Qu.:4587
Median :4775
Mean :4884
3rd Qu.:5089
Max. :6287 5.4 Graficar un raster
Para visualizar estos datos en R utilizando ggplot2, necesitamos convertirlos en un data.frame. El paquete terra tiene una función integrada para la conversión a un data.frame ploteable.
dem_chimborazo_df <- as.data.frame(dem_chimborazo, xy = TRUE)
ggplot() +
geom_raster(data = dem_chimborazo_df,
aes(x = x, y = y, fill = DEM_Chimborazo)) +
scale_fill_viridis_c() +
coord_quickmap()Warning in vp$just: encuentros parciales de 'just' to 'justification'