Información espacial en formato vector en R
Existen varios paquetes en R que permiten manejar datos espaciales, ya sea en formato de vector o raster. Algunos de ellos incluyen: sp, rgdal, rgeos, sf, stars, raster, terra. Sin embargo, en este curso nos enfocaremos en utilizar sf para el manejo de información vectorial y stars para el manejo de rasters.
Introducción a sf
sf es un paquete creado para trabajar de manera sencilla con información vectorial.
Features
El objeto básico en sf es un feature, es decir un rasgo o una característica. Este objeto usualmente contiene una geometría que describe la localización en el espacio de dicho rasgo y que puede contener más atributos que describen otras propiedades del mismo.
Existen varios tipos de features de acuerdo a su geometría, los básicos y más comunes son:
- Punto (Point).
- Línea (Linestring).
- Polígono (Polygons).
- Multipunto (Multipoint).
- Multilínea (Multilinestring).
- Mulitpolígono (Multipolygon).
- Colección de geometrías (Geometrycollection).
Ahora, la geometría de cada rasgo va a estar asociado a un sistema de coordenadas de referencia (CRS). Por lo cual, dicho CRS indicará la proyección de los datos y el datum.
Funciones en sf
Todas las funciones en sf estarán precedidas por st_ que significa spatial type.
Estructura de datos sf
Los features van a estar organizados en un formato de data.frame, es decir, de una tabla o cuadro de datos. Sin embargo, como las geometrías usualmente no contienen un solo valor, se guardarán en una lista.
Para ver un ejemplos carguemos un archivo shapefile.
library(sf)## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUEnc <- st_read(system.file("shape/nc.shp", package="sf"))## Reading layer `nc' from data source
## `D:\JonathanVSV\Documents\R\win-library\4.2\sf\shape\nc.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27nc## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27
## First 10 features:
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
## 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1
## 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
## 3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188 5
## 4 0.070 2.968 1831 1831 Currituck 37053 37053 27 508 1
## 5 0.153 2.206 1832 1832 Northampton 37131 37131 66 1421 9
## 6 0.097 1.670 1833 1833 Hertford 37091 37091 46 1452 7
## 7 0.062 1.547 1834 1834 Camden 37029 37029 15 286 0
## 8 0.091 1.284 1835 1835 Gates 37073 37073 37 420 0
## 9 0.118 1.421 1836 1836 Warren 37185 37185 93 968 4
## 10 0.124 1.428 1837 1837 Stokes 37169 37169 85 1612 1
## NWBIR74 BIR79 SID79 NWBIR79 geometry
## 1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3...
## 2 10 542 3 12 MULTIPOLYGON (((-81.23989 3...
## 3 208 3616 6 260 MULTIPOLYGON (((-80.45634 3...
## 4 123 830 2 145 MULTIPOLYGON (((-76.00897 3...
## 5 1066 1606 3 1197 MULTIPOLYGON (((-77.21767 3...
## 6 954 1838 5 1237 MULTIPOLYGON (((-76.74506 3...
## 7 115 350 2 139 MULTIPOLYGON (((-76.00897 3...
## 8 254 594 2 371 MULTIPOLYGON (((-76.56251 3...
## 9 748 1190 2 844 MULTIPOLYGON (((-78.30876 3...
## 10 160 2038 5 176 MULTIPOLYGON (((-80.02567 3...Aquí podemos ver la tabla de atributos del archivo, así como el tipo de geometría (multipolígono) y las geometrías de cada entrada. A partir de esta información podemos ver que se trata de una tabla con 100 rasgos con 14 atributos + 1 columna con las geometrías de cada rasgo. Cada renglón de la tabla constituye un simple feature (sf), es decir un rasgo con ciertos atributos y su geometría. La geometría de cada rasgo constituye un simple feature geometry (sfg), mientras que la columna con todas las geometrías de los rasgos, una simple feature geometry list-column (sfc).
Conversión de datos tabulados a sf
Un tipo de conversión bastante común es que contemos con una tabla de datos con coordenadas geográficas y ciertos atributos extra que deseamos convertir a sf. Para realizar dicha operación se puede usar la función st_as_sf.
Primero se creará un conjunto de datos de coordenadas y atributos.
tabla <- data.frame(lon = c(-101.2737, -101.0627, -101.15),
lat = c(19.99155, 19.96613, 19.98),
Atributo1 = c("A", "B", "C"))
tabla## lon lat Atributo1
## 1 -101.2737 19.99155 A
## 2 -101.0627 19.96613 B
## 3 -101.1500 19.98000 CDespués se convierte a sf. Para ello, hay que indicar el CRS de la información, para lo cual se puede consultar los códigos EPSG en https://epsg.io/.
vec1 <- st_as_sf(tabla,
# Nombre de las columnas con coordenadas x y
coords = c( "lon", "lat"),
# CRS en código EPSG
crs = 4326,
# ¿Quitar columnas con información espacial?
remove = T)
vec1## Simple feature collection with 3 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -101.2737 ymin: 19.96613 xmax: -101.0627 ymax: 19.99155
## Geodetic CRS: WGS 84
## Atributo1 geometry
## 1 A POINT (-101.2737 19.99155)
## 2 B POINT (-101.0627 19.96613)
## 3 C POINT (-101.15 19.98)Ahora supongamos que queremos hacer un polígono a partir de ese vector.
poly1 <- vec1 %>%
# Si se hacen varios polígonos de acuerdo a un atributo
# dplyr::group_by() %>%
dplyr::summarise() %>%
st_cast("POLYGON")
poly1## Simple feature collection with 1 feature and 0 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -101.2737 ymin: 19.96613 xmax: -101.0627 ymax: 19.99155
## Geodetic CRS: WGS 84
## geometry
## 1 POLYGON ((-101.15 19.98, -1...Lectura y escritura de datos
Para leer y escribir datos desde archivos externos se utilizan las funciones st_read y st_write. Por ejemplo, carguemos un archivo que viene en el paquete sf:
filename <- system.file("shape/nc.shp", package="sf")
nc <- st_read(filename)## Reading layer `nc' from data source
## `D:\JonathanVSV\Documents\R\win-library\4.2\sf\shape\nc.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27nc## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27
## First 10 features:
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
## 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1
## 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
## 3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188 5
## 4 0.070 2.968 1831 1831 Currituck 37053 37053 27 508 1
## 5 0.153 2.206 1832 1832 Northampton 37131 37131 66 1421 9
## 6 0.097 1.670 1833 1833 Hertford 37091 37091 46 1452 7
## 7 0.062 1.547 1834 1834 Camden 37029 37029 15 286 0
## 8 0.091 1.284 1835 1835 Gates 37073 37073 37 420 0
## 9 0.118 1.421 1836 1836 Warren 37185 37185 93 968 4
## 10 0.124 1.428 1837 1837 Stokes 37169 37169 85 1612 1
## NWBIR74 BIR79 SID79 NWBIR79 geometry
## 1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3...
## 2 10 542 3 12 MULTIPOLYGON (((-81.23989 3...
## 3 208 3616 6 260 MULTIPOLYGON (((-80.45634 3...
## 4 123 830 2 145 MULTIPOLYGON (((-76.00897 3...
## 5 1066 1606 3 1197 MULTIPOLYGON (((-77.21767 3...
## 6 954 1838 5 1237 MULTIPOLYGON (((-76.74506 3...
## 7 115 350 2 139 MULTIPOLYGON (((-76.00897 3...
## 8 254 594 2 371 MULTIPOLYGON (((-76.56251 3...
## 9 748 1190 2 844 MULTIPOLYGON (((-78.30876 3...
## 10 160 2038 5 176 MULTIPOLYGON (((-80.02567 3...Para escribir un vector al disco:
st_write(nc,
"nc_exp.shp",
# sobreescribir si ya existe
append = F)## Deleting layer `nc_exp' using driver `ESRI Shapefile'
## Writing layer `nc_exp' to data source `nc_exp.shp' using driver `ESRI Shapefile'
## Writing 100 features with 14 fields and geometry type Multi Polygon.Las funciones st_read y st_write contienen más argumentos para definir
Visualización de datos
Por último, el visualizar los datos nos puede dar una muy buena idea de los productos intermedios en un flujo de trabajo o verificar que la información que importamos o exportamos es la correcta. Para ver el objecto sf con todos sus atributos:
plot(nc)## Warning: plotting the first 10 out of 14 attributes; use max.plot = 14 to plot
## all
O para ver algún atributo en particular.
plot(nc["AREA"])
Operaciones espaciales básicas
Descargar información de trabajo
Primero vamos a cargar información desde mi github.
## [1] "D:\\Descargas\\SIGmaterial\\16mun.zip"
## [2] "D:\\Descargas\\SIGmaterial\\MX.zip"
## [3] "D:\\Descargas\\SIGmaterial\\Mich.zip"
## [4] "D:\\Descargas\\SIGmaterial\\SerieVI_Mich.zip"## [,1]
## [1,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/16mun.cpg"
## [2,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/16mun.dbf"
## [3,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/16mun.prj"
## [4,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/16mun.shp"
## [5,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/16mun.shx"
## [,2]
## [1,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/Mexico.shx"
## [2,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/Mexico.cpg"
## [3,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/Mexico.dbf"
## [4,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/Mexico.prj"
## [5,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/Mexico.shp"
## [,3]
## [1,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/Mich.cpg"
## [2,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/Mich.dbf"
## [3,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/Mich.prj"
## [4,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/Mich.shp"
## [5,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/Mich.shx"
## [,4]
## [1,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/SerieVI_Mich.cpg"
## [2,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/SerieVI_Mich.dbf"
## [3,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/SerieVI_Mich.prj"
## [4,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/SerieVI_Mich.shp"
## [5,] "D:/Drive/Jonathan_trabaggio/Doctorado/R/SIG/Curso/Data/SerieVI_Mich.shx"Veamos qué descargamos.
# Ver qué archivos en la carpeta indicada cumplen con el patrón de ".shp"
files <- list.files(paste0(getwd(), "/Data"), ".shp")
files## [1] "16mun.shp" "Mexico.shp" "Mich.shp" "roi.shp"
## [5] "SerieVI_Mich.shp"# Leer los archivos
mun <- st_read(paste0(getwd(), "/Data/","16mun.shp"))## Reading layer `16mun' from data source
## `D:\Drive\Jonathan_trabaggio\Doctorado\R\SIG\Curso\Data\16mun.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 113 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -50.89613 ymin: 1983660 xmax: 388695.4 ymax: 2258033
## Projected CRS: WGS 84 / UTM zone 14Nmx <- st_read(paste0(getwd(), "/Data/","Mexico.shp"))## Reading layer `Mexico' from data source
## `D:\Drive\Jonathan_trabaggio\Doctorado\R\SIG\Curso\Data\Mexico.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -118.4042 ymin: 14.55055 xmax: -86.73862 ymax: 32.71846
## Geodetic CRS: WGS 84mich <- st_read(paste0(getwd(), "/Data/","Mich.shp"))## Reading layer `Mich' from data source
## `D:\Drive\Jonathan_trabaggio\Doctorado\R\SIG\Curso\Data\Mich.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -103.7455 ymin: 17.92189 xmax: -100.057 ymax: 20.40305
## Geodetic CRS: WGS 84veg <- st_read(paste0(getwd(), "/Data/","SerieVI_Mich.shp"))## Reading layer `SerieVI_Mich' from data source
## `D:\Drive\Jonathan_trabaggio\Doctorado\R\SIG\Curso\Data\SerieVI_Mich.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 6529 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -103.7381 ymin: 17.91491 xmax: -100.063 ymax: 20.39456
## Geodetic CRS: WGS 84mun## Simple feature collection with 113 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -50.89613 ymin: 1983660 xmax: 388695.4 ymax: 2258033
## Projected CRS: WGS 84 / UTM zone 14N
## First 10 features:
## CVEGEO CVE_ENT CVE_MUN NOMGEO geometry
## 1 16001 16 001 Acuitzio MULTIPOLYGON (((256032.1 21...
## 2 16002 16 002 Aguililla MULTIPOLYGON (((113031.8 21...
## 3 16003 16 003 Álvaro Obregón MULTIPOLYGON (((290455.6 22...
## 4 16004 16 004 Angamacutiro MULTIPOLYGON (((222887 2240...
## 5 16005 16 005 Angangueo MULTIPOLYGON (((365943.8 21...
## 6 16006 16 006 Apatzingán MULTIPOLYGON (((160165.4 21...
## 7 16007 16 007 Aporo MULTIPOLYGON (((355615.2 21...
## 8 16008 16 008 Aquila MULTIPOLYGON (((20860.2 207...
## 9 16009 16 009 Ario MULTIPOLYGON (((221485.2 21...
## 10 16010 16 010 Arteaga MULTIPOLYGON (((173820.4 20...mx## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -118.4042 ymin: 14.55055 xmax: -86.73862 ymax: 32.71846
## Geodetic CRS: WGS 84
## COUNTRY geometry
## 1 Mexico MULTIPOLYGON (((-97.77687 2...mich## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -103.7455 ymin: 17.92189 xmax: -100.057 ymax: 20.40305
## Geodetic CRS: WGS 84
## CODIGO ESTADO geometry
## 1 MX16 Michoacán POLYGON ((-103.4796 18.9672...veg## Simple feature collection with 6529 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -103.7381 ymin: 17.91491 xmax: -100.063 ymax: 20.39456
## Geodetic CRS: WGS 84
## First 10 features:
## TIP_INFO CLAVE TIP_ECOV
## 1 OTRO RASGO AH NO APLICABLE
## 2 OTRO RASGO AH NO APLICABLE
## 3 OTRO RASGO AH NO APLICABLE
## 4 OTRO RASGO AH NO APLICABLE
## 5 OTRO RASGO AH NO APLICABLE
## 6 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BA BOSQUE DE CONÍFERAS
## 7 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BP BOSQUE DE CONÍFERAS
## 8 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BP BOSQUE DE CONÍFERAS
## 9 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BPQ BOSQUE DE CONÍFERAS
## 10 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BPQ BOSQUE DE CONÍFERAS
## TIP_VEG DESVEG FASE_VS OTROS CAL_POS
## 1 NO APLICABLE NO APLICABLE NO APLICABLE URBANO CONSTRUIDO Aproximada
## 2 NO APLICABLE NO APLICABLE NO APLICABLE URBANO CONSTRUIDO Aproximada
## 3 NO APLICABLE NO APLICABLE NO APLICABLE URBANO CONSTRUIDO Aproximada
## 4 NO APLICABLE NO APLICABLE NO APLICABLE URBANO CONSTRUIDO Aproximada
## 5 NO APLICABLE NO APLICABLE NO APLICABLE URBANO CONSTRUIDO Aproximada
## 6 BOSQUE DE OYAMEL PRIMARIA NINGUNO NO APLICABLE Aproximada
## 7 BOSQUE DE PINO PRIMARIA NINGUNO NO APLICABLE Aproximada
## 8 BOSQUE DE PINO PRIMARIA NINGUNO NO APLICABLE Aproximada
## 9 BOSQUE DE PINO-ENCINO PRIMARIA NINGUNO NO APLICABLE Aproximada
## 10 BOSQUE DE PINO-ENCINO PRIMARIA NINGUNO NO APLICABLE Aproximada
## geometry
## 1 MULTIPOLYGON (((-101.3781 1...
## 2 MULTIPOLYGON (((-101.3348 1...
## 3 MULTIPOLYGON (((-101.3585 1...
## 4 MULTIPOLYGON (((-101.3107 1...
## 5 MULTIPOLYGON (((-101.4364 1...
## 6 MULTIPOLYGON (((-101.42 19....
## 7 MULTIPOLYGON (((-101.2158 1...
## 8 MULTIPOLYGON (((-100.7489 1...
## 9 MULTIPOLYGON (((-101.3228 1...
## 10 MULTIPOLYGON (((-101.3273 1...Ver los gráficos de la información.
plot(mun)
plot(mx)
plot(mich)
# Como está pesada esta capa mostrar sólo un atributo
plot(veg["TIP_ECOV"])
Operaciones basadas en geometrías
En sf existe un conjunto de funciones diseñadas para realizar operaciones con las geometrías de los objetos sf. Aquí revisaremos algunas de las más utilizadas.
Comprobación de geometrías válidas
Comprobar que contiene geometrías válidas
st_is_valid(mx)## [1] TRUEEn caso de que haya geometrías no válidas se puede forzar la eliminación de geometrías no válidas mediante st_make_valid.
Reproyección
Para reproyectar una capa a otro CRS simplemente se puede indicar el CRS objetivo o utilizar otra capa para extraer el CRS objetivo mediante st_crs. En este ejemplo, primero veamos el CRS de las capas mun y las demás.
st_crs(mun)## Coordinate Reference System:
## User input: WGS 84 / UTM zone 14N
## wkt:
## PROJCRS["WGS 84 / UTM zone 14N",
## 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 14N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-99,
## 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",0,
## 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",32614]]st_crs(mx)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]st_crs(mich)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]st_crs(veg)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]Vemos que mun está en otro CRS. Para transformar la capa al mismo CRS que las demás capas.
mun <- st_transform(mun, st_crs(mx))
st_crs(mun)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]Vemos que ya se realizó la transformación deseada.
Cálculo de área
Para calcular el área:
# Normalmente el valor de área está en m^2
st_area(mx)## 1.961269e+12 [m^2]st_area(mich)## 59654936304 [m^2]Recuerden que estamos trabajando con datos en coordenadas geográficas (no proyectados), por lo que sería mejor idea primero proyectarlos y después calcular su área para obtener medidas cómo el área o distancia. Esto se podría hacer mediante:
# Normalmente el valor de área está en m^2
mx |>
st_transform(32614) |>
st_area()## 1.98135e+12 [m^2]Más adelante explicaremos a qué se refiere el operador |> pero por el momento lo podemos interpretar como una función que permite encadenar procesos.
Intersección
Revisar si existe intersección entre dos sf.
st_intersects(mx, mich)## Sparse geometry binary predicate list of length 1, where the predicate
## was `intersects'
## 1: 1Retorna un valor por cada feature que contiene cada objeto. Si la salida es (empty) implica que no se intersectan, mientras que cuando la salida contiene números, estos indican para cada feature del objeto 1 con qué features del objeto2 se intersectan.
Por el contrario, si se desea extraer la intersección entre dos objetos se usa st_intersection
int_mx_mich <- st_intersection(mx, mich)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries¿Qué creen que será el resultado?
plot(int_mx_mich)
int_mx_mich## Simple feature collection with 1 feature and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -103.7455 ymin: 17.92189 xmax: -100.057 ymax: 20.40305
## Geodetic CRS: WGS 84
## COUNTRY CODIGO ESTADO geometry
## 1 Mexico MX16 Michoacán POLYGON ((-103.5775 18.8816...Noten que la tabla de atributos contiene los atributos del objeto 1 (mx) y del objeto 2 (mich).
Diferencia
Para calcular la diferencia:
mx_sin_mich <- st_difference(mx, mich)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesplot(mx_sin_mich["COUNTRY"])
Unión
Para unir varias geometrías en una sola.
mx_con_mich <- st_union(mx_sin_mich, mich)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesplot(mx_con_mich["COUNTRY"])
Unión de polígonos o varios sf en una misma colección
Para unir varios polígonos o rasgos espaciales dentro de una misma colección se utiliza la notación de tidyverse para unir tablas, es decir, bind_rows.
mx_con_mich_sinunion <- dplyr::bind_rows(mx_sin_mich, mich)
mx_con_mich_sinunion## Simple feature collection with 2 features and 3 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -118.4042 ymin: 14.55055 xmax: -86.73862 ymax: 32.71846
## Geodetic CRS: WGS 84
## COUNTRY CODIGO ESTADO geometry
## 1 Mexico MX16 Michoacán MULTIPOLYGON (((-103.5508 1...
## 2 <NA> MX16 Michoacán POLYGON ((-103.4796 18.9672...plot(mx_con_mich_sinunion["COUNTRY"])
Regresa Michoacán al polígono completo de México, aunque noten que quedan algunos polígonos basura.
Centroides
Se puede calcular el centroide de un polígono.
mich_centroid <- st_centroid(mich)## Warning in st_centroid.sf(mich): st_centroid assumes attributes are constant
## over geometries of xplot(mich_centroid)
Buffers
Además, se puede calcular buffers.
# Normalmente el valor de buffer está en m
mich_centr_buff <- st_buffer(mich_centroid,
10000)
plot(mich_centr_buff)
Join espacial
Muestreo aleatorio
Esta función es bastante utilizada para muestrear una superficie, lo cual es muy utilizado para verificar mapas clasificados u productos similares. Para realizar este muestreo se utiliza la función st_sample.
Operaciones basadas en atributos
En sf existe un conjunto de funciones diseñadas para realizar operaciones con basados en los atributos de los objetos sf. Aquí revisaremos algunas de las más utilizadas.
Tidyverse-esque sintaxis
En R existe un paquete llamado tidyverse que contiene un conjunto de paquetes muy útiles para gráficas (ggplot2), limpiar y arreglar bases de datos (tidyr y dplyr), los cuales hacen uso del operador %>% o pipe para concatenar varias funciones de manera más resumida. Esta sintaxis será de mucha ayuda para concatenar un flujo de trabajo o para simplificar el código, lo cual muchas veces facilita su lectura. Por ello, para la mayoría de las operaciones basadas en atributos usaremos el operador %>% para realizar las operaciones. Actualmente, en R > 4.0.0, dicho operador fue incluido en R base como |>.
Selección por rasgos
Primero revisemos parte de los datos, así como los niveles de un atributo. Para consultar los niveles de un atributo lo podemos hacer con la sintaxis tradicional de R o a la manera tidyverse.
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.1 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x readr::parse_date() masks curl::parse_date()head(veg)## Simple feature collection with 6 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -101.4439 ymin: 19.40706 xmax: -101.3107 ymax: 19.54956
## Geodetic CRS: WGS 84
## TIP_INFO CLAVE TIP_ECOV TIP_VEG
## 1 OTRO RASGO AH NO APLICABLE NO APLICABLE
## 2 OTRO RASGO AH NO APLICABLE NO APLICABLE
## 3 OTRO RASGO AH NO APLICABLE NO APLICABLE
## 4 OTRO RASGO AH NO APLICABLE NO APLICABLE
## 5 OTRO RASGO AH NO APLICABLE NO APLICABLE
## 6 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BA BOSQUE DE CONÍFERAS BOSQUE DE OYAMEL
## DESVEG FASE_VS OTROS CAL_POS
## 1 NO APLICABLE NO APLICABLE URBANO CONSTRUIDO Aproximada
## 2 NO APLICABLE NO APLICABLE URBANO CONSTRUIDO Aproximada
## 3 NO APLICABLE NO APLICABLE URBANO CONSTRUIDO Aproximada
## 4 NO APLICABLE NO APLICABLE URBANO CONSTRUIDO Aproximada
## 5 NO APLICABLE NO APLICABLE URBANO CONSTRUIDO Aproximada
## 6 PRIMARIA NINGUNO NO APLICABLE Aproximada
## geometry
## 1 MULTIPOLYGON (((-101.3781 1...
## 2 MULTIPOLYGON (((-101.3348 1...
## 3 MULTIPOLYGON (((-101.3585 1...
## 4 MULTIPOLYGON (((-101.3107 1...
## 5 MULTIPOLYGON (((-101.4364 1...
## 6 MULTIPOLYGON (((-101.42 19....# Ver niveles únicos del atributo TIP_ECOV
unique(veg$TIP_ECOV)## [1] "NO APLICABLE" "BOSQUE DE CONÍFERAS"
## [3] "BOSQUE DE ENCINO" "VEGETACIÓN INDUCIDA"
## [5] "SELVA CADUCIFOLIA" "PASTIZAL"
## [7] "VEGETACIÓN HIDRÓFILA" "MATORRAL XERÓFILO"
## [9] "SELVA SUBCADUCIFOLIA" "BOSQUE MESÓFILO DE MONTAÑA"
## [11] "SELVA ESPINOSA"# Lo mismo pero en forma tidyverse
veg |>
select(TIP_ECOV) |>
distinct(TIP_ECOV)## TIP_ECOV
## 1 NO APLICABLE
## 2 BOSQUE DE CONÍFERAS
## 3 BOSQUE DE ENCINO
## 4 VEGETACIÓN INDUCIDA
## 5 SELVA CADUCIFOLIA
## 6 PASTIZAL
## 7 VEGETACIÓN HIDRÓFILA
## 8 MATORRAL XERÓFILO
## 9 SELVA SUBCADUCIFOLIA
## 10 BOSQUE MESÓFILO DE MONTAÑA
## 11 SELVA ESPINOSASeleccionemos entonces todos los polígonos de selva caducifolia.
# Lo mismo pero en forma tidyverse
sbc <- veg |>
filter(TIP_ECOV == "SELVA CADUCIFOLIA")
sbc## Simple feature collection with 1286 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -103.7093 ymin: 17.99039 xmax: -100.2101 ymax: 20.37962
## Geodetic CRS: WGS 84
## First 10 features:
## TIP_INFO CLAVE TIP_ECOV
## 1 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 2 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSA/SBC SELVA CADUCIFOLIA
## 3 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 4 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 5 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 6 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 7 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 8 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 9 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 10 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## TIP_VEG DESVEG FASE_VS OTROS CAL_POS
## 1 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 2 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBÓREA NO APLICABLE Aproximada
## 3 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 4 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 5 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 6 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 7 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 8 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 9 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 10 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## geometry
## 1 MULTIPOLYGON (((-101.3513 1...
## 2 MULTIPOLYGON (((-101.4022 1...
## 3 MULTIPOLYGON (((-102.8018 1...
## 4 MULTIPOLYGON (((-102.8374 1...
## 5 MULTIPOLYGON (((-102.7461 1...
## 6 MULTIPOLYGON (((-102.7968 1...
## 7 MULTIPOLYGON (((-102.685 18...
## 8 MULTIPOLYGON (((-102.7767 1...
## 9 MULTIPOLYGON (((-102.7778 1...
## 10 MULTIPOLYGON (((-102.7534 1...# Para hacer lo mismo en sintaxis tradicional
veg[veg$TIP_ECOV == "SELVA CADUCIFOLIA",]## Simple feature collection with 1286 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -103.7093 ymin: 17.99039 xmax: -100.2101 ymax: 20.37962
## Geodetic CRS: WGS 84
## First 10 features:
## TIP_INFO CLAVE TIP_ECOV
## 66 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 67 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSA/SBC SELVA CADUCIFOLIA
## 192 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 193 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 194 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 262 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 263 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 264 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 265 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 266 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## TIP_VEG DESVEG FASE_VS OTROS CAL_POS
## 66 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 67 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBÓREA NO APLICABLE Aproximada
## 192 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 193 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 194 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 262 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 263 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 264 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 265 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 266 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## geometry
## 66 MULTIPOLYGON (((-101.3513 1...
## 67 MULTIPOLYGON (((-101.4022 1...
## 192 MULTIPOLYGON (((-102.8018 1...
## 193 MULTIPOLYGON (((-102.8374 1...
## 194 MULTIPOLYGON (((-102.7461 1...
## 262 MULTIPOLYGON (((-102.7968 1...
## 263 MULTIPOLYGON (((-102.685 18...
## 264 MULTIPOLYGON (((-102.7767 1...
## 265 MULTIPOLYGON (((-102.7778 1...
## 266 MULTIPOLYGON (((-102.7534 1...A partir de ahora utilizaremos únicamente la notación “tidy”. También podemos seleccionar entradas a partir de uno o más niveles de algún atributo. El operador %in% indica que se busquen las entradas donde su atributo TIP_ECOV sea alguna de las entradas indicadas en el vector siguiente, es decir: SELVA CADUCIFOLIA o BOSQUE DE CONÍFERAS.
# Lo mismo pero en forma tidyverse
sbc_bc <- veg |>
filter(TIP_ECOV %in% c("SELVA CADUCIFOLIA", "BOSQUE DE CONÍFERAS"))
sbc_bc## Simple feature collection with 2346 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -103.7093 ymin: 17.99039 xmax: -100.1266 ymax: 20.37962
## Geodetic CRS: WGS 84
## First 10 features:
## TIP_INFO CLAVE TIP_ECOV
## 1 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BA BOSQUE DE CONÍFERAS
## 2 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BP BOSQUE DE CONÍFERAS
## 3 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BP BOSQUE DE CONÍFERAS
## 4 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BPQ BOSQUE DE CONÍFERAS
## 5 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BPQ BOSQUE DE CONÍFERAS
## 6 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BPQ BOSQUE DE CONÍFERAS
## 7 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BPQ BOSQUE DE CONÍFERAS
## 8 ECOLÓGICA-FLORÍSTICA-FISONÓMICA BPQ BOSQUE DE CONÍFERAS
## 9 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/BA BOSQUE DE CONÍFERAS
## 10 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSA/BP BOSQUE DE CONÍFERAS
## TIP_VEG DESVEG FASE_VS OTROS CAL_POS
## 1 BOSQUE DE OYAMEL PRIMARIA NINGUNO NO APLICABLE Aproximada
## 2 BOSQUE DE PINO PRIMARIA NINGUNO NO APLICABLE Aproximada
## 3 BOSQUE DE PINO PRIMARIA NINGUNO NO APLICABLE Aproximada
## 4 BOSQUE DE PINO-ENCINO PRIMARIA NINGUNO NO APLICABLE Aproximada
## 5 BOSQUE DE PINO-ENCINO PRIMARIA NINGUNO NO APLICABLE Aproximada
## 6 BOSQUE DE PINO-ENCINO PRIMARIA NINGUNO NO APLICABLE Aproximada
## 7 BOSQUE DE PINO-ENCINO PRIMARIA NINGUNO NO APLICABLE Aproximada
## 8 BOSQUE DE PINO-ENCINO PRIMARIA NINGUNO NO APLICABLE Aproximada
## 9 BOSQUE DE OYAMEL SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 10 BOSQUE DE PINO SECUNDARIA ARBÓREA NO APLICABLE Aproximada
## geometry
## 1 MULTIPOLYGON (((-101.42 19....
## 2 MULTIPOLYGON (((-101.2158 1...
## 3 MULTIPOLYGON (((-100.7489 1...
## 4 MULTIPOLYGON (((-101.3228 1...
## 5 MULTIPOLYGON (((-101.3273 1...
## 6 MULTIPOLYGON (((-101.2634 1...
## 7 MULTIPOLYGON (((-101.4312 1...
## 8 MULTIPOLYGON (((-101.2528 1...
## 9 MULTIPOLYGON (((-101.4514 1...
## 10 MULTIPOLYGON (((-101.3691 1...Calcular nuevos atributos
Para agregar una nueva columna o un nuevo atributo se utiliza la función mutate seguido del nombre del nuevo atributo y de la definición del nuevo atributo de acuerdo a una función o un único valor.
sbc <- sbc |>
# as.numeric se utiliza para convertir los valores a numéricos y quitarle la propiedad de units (ver sección cálculo de área). Además, al dividir entre 10000 se pasa de m2 a ha.
mutate(area_ha = as.numeric(st_area(sbc)/10000))
sbc## Simple feature collection with 1286 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -103.7093 ymin: 17.99039 xmax: -100.2101 ymax: 20.37962
## Geodetic CRS: WGS 84
## First 10 features:
## TIP_INFO CLAVE TIP_ECOV
## 1 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 2 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSA/SBC SELVA CADUCIFOLIA
## 3 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 4 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 5 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 6 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 7 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 8 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 9 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 10 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## TIP_VEG DESVEG FASE_VS OTROS CAL_POS
## 1 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 2 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBÓREA NO APLICABLE Aproximada
## 3 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 4 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 5 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 6 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 7 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 8 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 9 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 10 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## geometry area_ha
## 1 MULTIPOLYGON (((-101.3513 1... 57.60691
## 2 MULTIPOLYGON (((-101.4022 1... 156.68668
## 3 MULTIPOLYGON (((-102.8018 1... 125.93153
## 4 MULTIPOLYGON (((-102.8374 1... 557.72104
## 5 MULTIPOLYGON (((-102.7461 1... 2433.63242
## 6 MULTIPOLYGON (((-102.7968 1... 122.14498
## 7 MULTIPOLYGON (((-102.685 18... 2680.84067
## 8 MULTIPOLYGON (((-102.7767 1... 31.28422
## 9 MULTIPOLYGON (((-102.7778 1... 331.41608
## 10 MULTIPOLYGON (((-102.7534 1... 102.77817Selección por rasgos2
Ahora vamos a filtrar usando variables numéricas. En este caso, para quedarnos con entradas con un área entre 100 y 5000 ha.
sbc_filt <- sbc |>
filter(area_ha <= 5000 & area_ha >= 100)
sbc_filt## Simple feature collection with 855 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -103.7093 ymin: 17.9907 xmax: -100.2101 ymax: 20.37962
## Geodetic CRS: WGS 84
## First 10 features:
## TIP_INFO CLAVE TIP_ECOV
## 1 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSA/SBC SELVA CADUCIFOLIA
## 2 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 3 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 4 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 5 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 6 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 7 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 8 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 9 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 10 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## TIP_VEG DESVEG FASE_VS OTROS CAL_POS
## 1 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBÓREA NO APLICABLE Aproximada
## 2 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 3 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 4 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 5 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 6 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 7 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 8 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 9 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 10 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## geometry area_ha
## 1 MULTIPOLYGON (((-101.4022 1... 156.6867
## 2 MULTIPOLYGON (((-102.8018 1... 125.9315
## 3 MULTIPOLYGON (((-102.8374 1... 557.7210
## 4 MULTIPOLYGON (((-102.7461 1... 2433.6324
## 5 MULTIPOLYGON (((-102.7968 1... 122.1450
## 6 MULTIPOLYGON (((-102.685 18... 2680.8407
## 7 MULTIPOLYGON (((-102.7778 1... 331.4161
## 8 MULTIPOLYGON (((-102.7534 1... 102.7782
## 9 MULTIPOLYGON (((-102.7617 1... 144.2841
## 10 MULTIPOLYGON (((-102.7179 1... 269.7079Finalmente, veamos una de las ventajas de usar la sintaxis “tidy”, al concatenar varios procesos para obtener el mismo resultado.
sbc_filt2 <- veg |>
filter(TIP_ECOV == "SELVA CADUCIFOLIA") |>
mutate(area_ha = as.numeric(st_area(sbc)/10000)) |>
filter(area_ha <= 5000 & area_ha >= 100)
sbc_filt2## Simple feature collection with 855 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -103.7093 ymin: 17.9907 xmax: -100.2101 ymax: 20.37962
## Geodetic CRS: WGS 84
## First 10 features:
## TIP_INFO CLAVE TIP_ECOV
## 1 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSA/SBC SELVA CADUCIFOLIA
## 2 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 3 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 4 ECOLÓGICA-FLORÍSTICA-FISONÓMICA SBC SELVA CADUCIFOLIA
## 5 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 6 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 7 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 8 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 9 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## 10 ECOLÓGICA-FLORÍSTICA-FISONÓMICA VSa/SBC SELVA CADUCIFOLIA
## TIP_VEG DESVEG FASE_VS OTROS CAL_POS
## 1 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBÓREA NO APLICABLE Aproximada
## 2 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 3 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 4 SELVA BAJA CADUCIFOLIA PRIMARIA NINGUNO NO APLICABLE Aproximada
## 5 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 6 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 7 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 8 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 9 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## 10 SELVA BAJA CADUCIFOLIA SECUNDARIA ARBUSTIVA NO APLICABLE Aproximada
## geometry area_ha
## 1 MULTIPOLYGON (((-101.4022 1... 156.6867
## 2 MULTIPOLYGON (((-102.8018 1... 125.9315
## 3 MULTIPOLYGON (((-102.8374 1... 557.7210
## 4 MULTIPOLYGON (((-102.7461 1... 2433.6324
## 5 MULTIPOLYGON (((-102.7968 1... 122.1450
## 6 MULTIPOLYGON (((-102.685 18... 2680.8407
## 7 MULTIPOLYGON (((-102.7778 1... 331.4161
## 8 MULTIPOLYGON (((-102.7534 1... 102.7782
## 9 MULTIPOLYGON (((-102.7617 1... 144.2841
## 10 MULTIPOLYGON (((-102.7179 1... 269.7079Fusionar por atributos (disolver)
mun_diss <- mun |>
# Agrupar por CVE_ENT, lo cual puede ser útil cuando se tienen datos más complejos
mutate(area_ha = as.numeric(st_area(mun)/10000)) |>
group_by(CVE_ENT) |>
# Resumir, se puede aprovechar para resumir un atributo
summarise(area_ha = sum(area_ha)) |>
st_cast("POLYGON")
mun_diss## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -103.7381 ymin: 17.91491 xmax: -100.063 ymax: 20.39456
## Geodetic CRS: WGS 84
## # A tibble: 1 x 3
## CVE_ENT area_ha geometry
## <chr> <dbl> <POLYGON [°]>
## 1 16 5877654. ((-101.4021 20.03894, -101.4021 20.03992, -101.4021 20.04014~plot(mun_diss)
Noten que los únicos atributos son la variable por la que fueron agrupados los datos y el atributo que se resumió.
Unir tablas por atributos
Otra función bastante utilizada es unir tablas no espaciales a objetos espaciales mediante sus atributos. Así que cargaremos unos datos del censo agropecuario del 2007 desde mi github.
cen_agrop <- read.csv("https://github.com/JonathanVSV/CursoSIG/raw/main/material/cag_2007_04.csv")Vemos que hizo parcialmente la unión. ¿Quién adivina por qué?
mun |>
left_join(cen_agrop,
by = c("NOMGEO" = "entidad_y_municipio"))## Simple feature collection with 113 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -103.7381 ymin: 17.91491 xmax: -100.063 ymax: 20.39456
## Geodetic CRS: WGS 84
## First 10 features:
## CVEGEO CVE_ENT CVE_MUN NOMGEO superficie_total_a
## 1 16001 16 001 Acuitzio 13430.637
## 2 16002 16 002 Aguililla 119237.074
## 3 16003 16 003 Álvaro Obregón NA
## 4 16004 16 004 Angamacutiro 15308.272
## 5 16005 16 005 Angangueo 2476.861
## 6 16006 16 006 Apatzingán NA
## 7 16007 16 007 Aporo 3617.149
## 8 16008 16 008 Aquila 70122.757
## 9 16009 16 009 Ario 48324.613
## 10 16010 16 010 Arteaga 236101.807
## regimen_de_tenencia_de_la_tierra_ejidal
## 1 1349.7480
## 2 41364.4160
## 3 NA
## 4 11273.3090
## 5 2075.7966
## 6 NA
## 7 886.5822
## 8 5915.8852
## 9 22239.7259
## 10 53052.6941
## regimen_de_tenencia_de_la_tierra_comunal
## 1 19.2265
## 2 0.0000
## 3 NA
## 4 4.5000
## 5 0.0000
## 6 NA
## 7 0.0000
## 8 27862.1605
## 9 0.0000
## 10 5028.4367
## regimen_de_tenencia_de_la_tierra_privada
## 1 12060.5622
## 2 77858.5831
## 3 NA
## 4 4026.4626
## 5 400.0645
## 6 NA
## 7 2721.5673
## 8 36344.7113
## 9 26078.8871
## 10 176369.9331
## regimen_de_tenencia_de_la_tierra_de_colonia
## 1 0
## 2 0
## 3 NA
## 4 0
## 5 0
## 6 NA
## 7 0
## 8 0
## 9 0
## 10 0
## regimen_de_tenencia_de_la_tierra_publica unidad_de_medida
## 1 1.100 hectareas
## 2 14.075 hectareas
## 3 NA <NA>
## 4 4.000 hectareas
## 5 1.000 hectareas
## 6 NA <NA>
## 7 9.000 hectareas
## 8 0.000 hectareas
## 9 6.000 hectareas
## 10 1650.744 hectareas
## geometry
## 1 MULTIPOLYGON (((-101.3251 1...
## 2 MULTIPOLYGON (((-102.6761 1...
## 3 MULTIPOLYGON (((-101.0023 1...
## 4 MULTIPOLYGON (((-101.6524 2...
## 5 MULTIPOLYGON (((-100.2788 1...
## 6 MULTIPOLYGON (((-102.2316 1...
## 7 MULTIPOLYGON (((-100.3776 1...
## 8 MULTIPOLYGON (((-103.54 18....
## 9 MULTIPOLYGON (((-101.6504 1...
## 10 MULTIPOLYGON (((-102.0921 1...Al parecer el problema está en los caracteres con acentos. Vamos a tener que hacer dos cosas para arreglar los datos. 1) Substituir los acentos por letras sin acentos en mun y substituir la primera letra del nombre del municipio a su equivalente en mayúsculas en cen_agrop.
mun <- mun |>
# Aquí pasamos la fórmula al estilo purrr, otro paquete de R
# La lógica es indicar una función ~ aplicada a lo que viene antes, indicado por el .x
# Sustituimos todas las letras por su equivalente sin acento
mutate(across(NOMGEO, ~ str_replace_all(.x, c("Á" = "A", "á" = "a", "É" = "E", "é" = "e", "Í" = "I", "i" = "i", "Ó" = "O", "ó" = "o", "Ú" = "U", "u" = "u"))))
cen_agrop <- cen_agrop |>
# Vamos a tener que usar Regular Expressions para arreglar otro problema con los datos: algunos no empiezan con mayúsculas. Vamos a hacer esto para cambiar la primera letra de la primera palabra "^\\w{1}" a mayúscula
mutate(across(entidad_y_municipio, ~ str_replace_all(.x, "^\\w{1}", toupper)))Ahora ya podemos unir los datos de manera correcta.
mun <- mun |>
left_join(cen_agrop,
by = c("NOMGEO" = "entidad_y_municipio")) Una vez realizado estos pasos, ya obtenemos nuestra capa sf con los atributos de la tabla sin información espacial. Este ejemplo, sirve para demostrar algunas de las capacidades de R para automatizar procesos y limpiar datos. Además, es un ejemplo de la vida real, en el que las bases de datos muchas veces contienen errores o no hay una compatibilidad al 100 % entre distintas bases de datos.
Visualizacion avanzada
Para generar visualizaciones más avanzadas podemos utilizar otros paquetes como ggplot2 o tmap. Aquí veremos un ejemplo de cada uno.
ggplot2
En ggplot se puede determinar directamente el color y relleno de cada capa vectorial.
ggplot() +
# Graficar cada capa por separado con sus características de visualización
geom_sf(data = mich, fill = "gray70", colour = "red") +
geom_sf(data = mun, fill = "transparent", colour = "gray40") +
geom_sf(data = sbc, fill = "orange2")
Sin embargo, también se puede utilizar un atributo de cada capa para determinar el color o relleno de los polígonos
ggplot() +
# Graficar cada capa por separado con sus características de visualización
geom_sf(data = mich, fill = "gray70", colour = "red") +
geom_sf(data = mun, fill = "transparent", colour = "gray40") +
geom_sf(data = sbc_bc, aes(fill = TIP_ECOV)) +
# Para sobreescribir valores de relleno por default
scale_fill_manual(values = c("forestgreen", "orange2"))
tmap
La opción de tmap require de utilizar algún atributo de la información para determinar el color de relleno de cada polígono de acuerdo a los valores de ese atributo. Es similar a ggplot utilizando la opción de aes. Podemos ver algunas de las paletas que ya vienen pre hechas tanto para ggplot como para tmap con tmaptools::palette_explorer()
library(tmap)
# Cargar shapr
tm_shape(mich) +
# Elegir forma de visualización
tm_polygons() +
tm_shape(mun) +
tm_polygons() +
tm_shape(sbc_bc) +
tm_fill(col = "TIP_ECOV", palette = "-RdYlGn") 
</div>
</div> </div>