El código fuente de este documento está disponible en https://github.com/pf0953-programaciongeoespacialr-2020/leccion-09-r-datos-vectoriales-operaciones-espaciales.

Recursos de interés

Preparativos

Paquetes y datos para ejemplos:

# Paquete para manejo de datos vectoriales
library(sf)

# Paquete de Tidyverse para manipulación de datos
library(dplyr)

# Paquete con conjuntos de datos geoespaciales
library(spData)

# Paquete para mapas en la Web
library(leaflet)

Adicionalmente, se utilizan algunas capas geoespaciales de la Infraestructura Nacional de Datos Espaciales de Costa Rica (SNIT):

# URL base del servicio WFS IGN 1:5mil
url_base_wfs_ign_5mil <- "http://geos.snitcr.go.cr/be/IGN_5/wfs?"

# URL base del servicio WFS IGN 1:200mil
url_base_wfs_ign_200mil <- "http://geos.snitcr.go.cr/be/IGN_200/wfs?"

# URL base del servicio WFS del Sinac
url_base_wfs_sinac <- "http://geos1pne.sirefor.go.cr/wfs?"

# URL de las solicitudes de las capas
solicitud_provincias_wfs <- 
  "request=GetFeature&service=WFS&version=2.0.0&typeName=IGN_5:limiteprovincial_5k&outputFormat=application/json"

solicitud_cantones_wfs <- 
  "request=GetFeature&service=WFS&version=2.0.0&typeName=IGN_5:limitecantonal_5k&outputFormat=application/json"

solicitud_aerodromos_wfs <- 
  "request=GetFeature&service=WFS&version=1.0.0&typename=IGN_200:aerodromos_200k&outputFormat=application/json"

solicitud_redvial_wfs <- 
  "request=GetFeature&service=WFS&version=1.0.0&typename=IGN_200:redvial_200k&outputFormat=application/json"

solicitud_asp_wfs <- 
  "request=GetFeature&service=WFS&version=1.0.0&typename=PNE:areas_silvestres_protegidas&outputFormat=application/json"

# Recuperación y simplificación de las capas
# Provincias de Costa Rica
cr_provincias <-
  st_read(paste0(url_base_wfs_ign_5mil, solicitud_provincias_wfs)) %>%
  st_simplify(dTolerance = 1000)

# Cantones de Costa Rica
cr_cantones <-
  st_read(paste0(url_base_wfs_ign_5mil, solicitud_cantones_wfs)) %>%
  st_simplify(dTolerance = 1000)

# Áreas silvestres protegidas (ASP) de Costa Rica
cr_asp <-
  st_read(paste0(url_base_wfs_sinac, solicitud_asp_wfs)) %>%
  st_simplify(dTolerance = 1000)

# Red vial de Costa Rica
cr_redvial <-
  st_read(paste0(url_base_wfs_ign_200mil, solicitud_redvial_wfs)) %>%
  st_simplify(dTolerance = 1000)

# Aeródromos de Costa Rica
cr_aerodromos <-
  st_read(paste0(url_base_wfs_ign_200mil, solicitud_aerodromos_wfs))

También se incorpora una capa con registros de presencia de especies de Reptilia (clase de los reptiles), en formato CSV, con datos agrupados por la Infraestructura Mundial de Información en Biodiversidad (GBIF):

# Reptiles de Costa Rica
cr_reptilia <- 
  st_read(
    "https://raw.githubusercontent.com/pf0953-programaciongeoespacialr-2020/datos/master/biodiversidad/registros-presencia/cr/cr-reptilia.csv", 
    options=c(
      "X_POSSIBLE_NAMES=decimalLongitude",
      "Y_POSSIBLE_NAMES=decimalLatitude"
    )
  ) %>%
  st_set_crs(4326) %>% 
  st_transform(5367)

Visualización de las capas (1):

# Visualización de las capas (1)
plot(cr_provincias$geometry, axes=TRUE, graticule=TRUE, reset=FALSE)
plot(cr_redvial$geometry, col="grey", lwd=0.2, add=TRUE)
plot(cr_aerodromos$geometry, col="orange", pch=16, add=TRUE)

Visualización de las capas (2):

# Visualización de las capas (2)
plot(cr_asp$geometry, axes=TRUE, graticule=TRUE, reset=FALSE)
plot(cr_reptilia$geometry, col="green", pch=16, cex=0.1, add=TRUE)

Introducción

Esta lección brinda una visión general de las operaciones espaciales en datos vectoriales del paquete sf. Estas operaciones incluyen creación de subconjuntos espaciales (spatial subsetting), agregación espacial (spatial aggregation) y cruce de datos espaciales (spatial joining).

Creación de subconjuntos

Es el proceso de seleccionar objetos espaciales basados en su relación con otros objetos espaciales. Esta relación puede ser de intersección, contención, “roce” (touch), entre otras.

En el siguiente ejemplo, se seleccionan los puntos contenidos en un polígono:

# Aeródromos en Limón

# Geometría (polígonos) de la provincia de Limón
limon <-
  cr_provincias %>%
  filter(provincia == "Limón")

# Geometrías (puntos) de aeródromos contenidos en la geometría de Limón
aerodromos_limon <- cr_aerodromos[limon, , op = st_within]

# Mapeo
plot(limon$geometry, axes=TRUE, graticule=TRUE, reset=FALSE)
plot(aerodromos_limon$geometry, col="orange", pch=16, add=TRUE)

En el siguiente ejemplo, se seleccionan los polígonos contenidos en otro polígono:

# ASP contenidas en Limón

# Geometría (polígonos) de la provincia de Limón
limon <-
  cr_provincias %>%
  filter(provincia == "Limón")

# Geometrías (polígonos) de ASP contenidos en la geometría de Limón
asp_limon <- cr_asp[limon, , op = st_within]

# Mapeo
plot(limon$geometry, axes=TRUE, graticule=TRUE, reset=FALSE)
plot(asp_limon$geometry, col="green", add=TRUE)

# Nombres de las ASP
asp_limon %>%
  st_drop_geometry() %>%
  distinct(asp_limon$nombre_asp)
##         asp_limon$nombre_asp
## 1               Hitoy Cerere
## 2                  Limoncito
## 3 Acuiferos Guacimo y Pococi
## 4                 Tortuguero
## 5                 Rio Banano
## 6       Cuenca Rio Siquirres

La función st_within() es una de las que se conoce como un predicado geométrico binario. Estas funciones relacionan, de diferentes maneras, dos conjuntos de geometrías.

Nótese la diferencia al usar la función st_intersects() en lugar de st_within():

# ASP que se intersecan con la provincia de Limón

# Geometrías (polígonos) de ASP que se intersecan en la geometría de Limón
asp_limon <- cr_asp[limon, , op = st_intersects]

# Mapeo
plot(limon$geometry, axes=TRUE, graticule=TRUE, reset=FALSE)
plot(asp_limon$geometry, col="green", add=TRUE)

# Nombres de las ASP
asp_limon %>%
  st_drop_geometry() %>%
  distinct(asp_limon$nombre_asp)
##            asp_limon$nombre_asp
## 1                  Hitoy Cerere
## 2              Nacional Cariari
## 3              Braulio Carrillo
## 4                    Tortuguero
## 5                       Cahuita
## 6                      Chirripo
## 7      Internacional La Amistad
## 8                      Barbilla
## 9                   Rio Pacuare
## 10 Cordillera Volcanica Central
## 11                    Rio Macho
## 12           Barra del Colorado
## 13           Gandoca Manzanillo
## 14          Corredor Fronterizo
## 15                    Limoncito
## 16   Acuiferos Guacimo y Pococi
## 17                   Rio Banano
## 18         Cuenca Rio Siquirres

Sugerencia: pruebe el ejemplo anterior con las funciones st_overlaps(), st_disjoint() y st_touches().

Cruce de datos

El cruce “no espacial” de dos conjuntos de datos se basa en uno o varios campos (llamdos llaves) que tienen en común los cojuntos que van a cruzarse. Los cruces espaciales se basan en un principio similar, pero en lugar de campos comunes, la relación entre los conjuntos se realiza a través de una operación topológica, a veces llamada spatial overlay.

El siguiente ejemplo cruza los registros (puntos) de presencia de una especie con el conjunto de datos de cantones (polígonos), para obtener el nombre del cantón de cada registro.

# Registros de Bothrops asper (serpiente terciopelo) en Costa Rica
cr_bothrops_asper <-
  cr_reptilia %>%
  filter(species == "Bothrops asper")

# Cantones con registros de Bothrops asper
cr_cantones_bothrops_asper <- cr_cantones[cr_bothrops_asper, ]

# Cantidad de cantones con registros de Bothrops asper
nrow(cr_cantones_bothrops_asper)
## [1] 33
# Se "cruzan" los datos con la tabla de cantones, para obtener el nombre del cantón
cr_bothrops_asper_joined <- st_join(cr_bothrops_asper, cr_cantones["canton"])

# Despliegue de los datos cruzados
cr_bothrops_asper_joined[1:30, 
                         c("species", "canton", "stateProvince",
                           "decimalLongitude", "decimalLatitude"), 
                         drop=TRUE]
##           species        canton stateProvince decimalLongitude decimalLatitude
## 1  Bothrops asper        Quepos    Puntarenas        -84.16191        9.431868
## 2  Bothrops asper     Sarapiquí       Heredia        -84.00697       10.430623
## 3  Bothrops asper     Sarapiquí       Heredia        -84.00697       10.430623
## 4  Bothrops asper       Golfito    Puntarenas        -83.33660        8.403588
## 5  Bothrops asper     Siquirres         Limón        -83.54633       10.053453
## 6  Bothrops asper     Siquirres         Limón        -83.54633       10.053453
## 7  Bothrops asper        Pococí         Limón        -83.51974       10.571449
## 8  Bothrops asper     Siquirres         Limón        -83.54633       10.053453
## 9  Bothrops asper     Sarapiquí       Heredia        -84.00797       10.431169
## 10 Bothrops asper       Golfito    Puntarenas        -83.07377        8.423071
## 11 Bothrops asper       Golfito    Puntarenas        -83.42522        8.519863
## 12 Bothrops asper       Golfito    Puntarenas        -83.20215        8.700622
## 13 Bothrops asper     Turrialba       Cartago        -83.46435        9.972058
## 14 Bothrops asper    San Carlos      Alajuela        -84.76592       10.480142
## 15 Bothrops asper       Golfito    Puntarenas        -83.20246        8.698931
## 16 Bothrops asper          <NA>         Limón        -82.62003        9.624067
## 17 Bothrops asper          <NA>         Limón        -82.62002        9.623880
## 18 Bothrops asper     Sarapiquí       Heredia        -84.10306       10.361914
## 19 Bothrops asper           Osa    Puntarenas        -83.73276        8.624270
## 20 Bothrops asper    San Carlos      Alajuela        -84.76552       10.480279
## 21 Bothrops asper       Golfito    Puntarenas        -83.33791        8.409462
## 22 Bothrops asper       Golfito    Puntarenas        -83.33791        8.409462
## 23 Bothrops asper       Golfito    Puntarenas        -83.33791        8.409462
## 24 Bothrops asper       Golfito    Puntarenas        -83.33791        8.409462
## 25 Bothrops asper Pérez Zeledón      San José        -83.85162        9.289689
## 26 Bothrops asper       Golfito    Puntarenas        -83.33791        8.409462
## 27 Bothrops asper        Quepos    Puntarenas        -84.15538        9.409287
## 28 Bothrops asper         Upala      Alajuela        -85.01374       10.720918
## 29 Bothrops asper       Golfito    Puntarenas        -83.33659        8.411975
## 30 Bothrops asper       Golfito    Puntarenas        -83.33854        8.404359

La función st_join() realiza por defecto un lef_join, pero puede realizar cruces de otros tipos también. Por ejemplo, con el parámetro left=FALSE, puede realizarse un inner_join. Tambié por defecto, la operación topológica que se aplica es st_intersects().

Agregación

De manera similar al caso de la agregación de atributos, la agregación espacial es una forma de “condensar” o “resumir” datos. Los datos agregados muestran estadísticas de una variable (ej. promedio, suma) en relación con una variable de agrupación.

El siguiente ejemplo muestra el promedio de altitud de los sitios altos de Nueva Zelanda para cada región del país.

# Promedio de altitud de sitios altos para cada región de NZ
nz_avheight <- aggregate(x = nz_height, by = nz, FUN = mean)

# Mapa
plot(nz_avheight["elevation"], axes=TRUE, graticule=TRUE)

Un resultado similar puede lograrse con funciones tidy:

nz_avheight2 <- nz %>%
  st_join(nz_height) %>%
  group_by(Name) %>%
  summarize(elevation = mean(elevation, na.rm = TRUE))
## `summarise()` ungrouping output (override with `.groups` argument)
plot(nz_avheight2["elevation"], axes=TRUE, graticule=TRUE)

Ejercicios

  1. Despliegue la lista de especies de vipéridos (familia Viperidae) reportados en la provincia de Guanacaste. Despliegue también el mapa con las localidades.
# Localidades de Viperidae en Guanacaste
# Geometría (polígono) de la provincia de Puntarenas
guanacaste <-
  cr_provincias %>%
  filter(provincia == "Guanacaste")
# Registros de Viperidae en Costa Rica
cr_viperidae <-
  cr_reptilia %>%
  filter(family == "Viperidae")
# Geometrías (puntos) de Viperidae en Guanacaste
guanacaste_viperidae <- cr_viperidae[guanacaste, , op = st_within]
# Mapeo
plot(guanacaste$geometry, axes=TRUE, graticule=TRUE, reset=FALSE)
plot(guanacaste_viperidae$geometry, col="brown", pch=4, add=TRUE)

# Nombres de las especies
guanacaste_viperidae  %>%
  st_drop_geometry() %>%
  distinct(guanacaste_viperidae$species)
##    guanacaste_viperidae$species
## 1         Bothriechis lateralis
## 2        Bothriechis schlegelii
## 3        Porthidium ophryomegas
## 4                Crotalus simus
## 5      Agkistrodon howardgloydi
## 6            Porthidium nasutum
## 7                Bothrops asper
## 8          Atropoides mexicanus
## 9                              
## 10       Agkistrodon bilineatus
## 11            Crotalus durissus
## 12           Atropoides picadoi
  1. Presente el resultado anterior en un mapa de Leaflet.
guanacaste_viperidae_wgs84 <- 
  guanacaste_viperidae %>%
  st_transform(4326)

m <- leaflet() %>% 
  addProviderTiles(providers$Esri.WorldImagery, group = "Imágenes de ESRI") %>% 
  addProviderTiles(providers$Stamen.TonerLite, group = "Stamen Toner Lite") %>%
  addProviderTiles(providers$OpenStreetMap.Mapnik, group = "OpenStreetMap") %>%
  addCircleMarkers(data = guanacaste_viperidae_wgs84,
                   stroke = F,
                   radius = 4,
                   fillColor = 'red',
                   fillOpacity = 1,
                   popup = paste(guanacaste_viperidae_wgs84$genus, 
                                 guanacaste_viperidae_wgs84$species, 
                                 sep = '<br/>')
                   ) %>%
  addLayersControl(baseGroups = c("OpenStreetMap", "Stamen Toner Lite", "Imágenes de ESRI")) %>%
  addMiniMap(
    toggleDisplay = TRUE,
    tiles = providers$Stamen.TonerLite
  )
# Despliegue del mapa
m

  1. Repita el ejercicio de agregar el nombre del cantón a los registros de Bothrops asper, de manera que se incluyan solamente los registros en los que hay conincidencia (match).

  2. Repita el ejercicio de agregar los datos de los sitios altos de Nueva Zelanda, pero de acuerdon con el valor máximo.