Entrada y salida de geo-datos

Autor/a

Veronica Andreo

Fecha de publicación

31 de octubre de 2023

library(sf)
library(terra)
library(dplyr)
library(spData)

Introducción

Este capítulo trata sobre la lectura y escritura de datos geográficos. La entrada de datos geográficos es esencial para la geocomputación: las aplicaciones del mundo real son imposibles sin datos. La salida de datos también es vital, ya que permite a otros utilizar valiosos conjuntos de datos nuevos o mejorados resultantes de su trabajo. En conjunto, estos procesos de entrada/salida pueden denominarse E/S de datos (I/O en Inglés).

La E/S de datos geográficos suele realizarse con unas pocas líneas de código al principio y al final de los proyectos. A menudo se pasa por alto como un simple proceso de un solo paso. Sin embargo, los errores cometidos al principio de los proyectos (por ejemplo, utilizar un conjunto de datos obsoleto o defectuoso de alguna manera) pueden dar lugar a grandes problemas más adelante, por lo que vale la pena dedicar un tiempo considerable a identificar qué conjuntos de datos están disponibles, dónde se pueden encontrar y cómo recuperarlos. Por otra parte, existen muchos formatos de archivos geográficos, cada uno de los cuales tiene sus pros y sus contras, como así también diferentes formas de lectura y escritura eficiente.

Recuperación de datos abiertos

Cada vez hay más datos geográficos disponibles en Internet, muchos de los cuales son de acceso y uso gratuitos (siempre que se dé el crédito correspondiente a sus proveedores)1. En cierto modo, ahora hay demasiados datos, en el sentido de que a menudo hay varios sitios donde acceder al mismo conjunto de datos. Algunos conjuntos de datos son de mala calidad. En este contexto, es vital saber dónde buscar. Varios “geoportales” (servicios web que proporcionan conjuntos de datos geoespaciales, como Data.gov) son un buen punto de partida, ya que ofrecen una amplia gama de datos, pero a menudo sólo para lugares concretos (como se ilustra en la página actualizada de Wikipedia sobre el tema).

Algunos geoportales globales superan este problema. El portal GEOSS y el Copernicus Open Access Hub, por ejemplo, contienen muchos conjuntos de datos ráster con cobertura mundial. En el portal SEDAC, gestionado por la Administración Nacional de Aeronáutica y del Espacio (NASA), y en el geoportal INSPIRE de la Unión Europea, con cobertura mundial y regional, se puede acceder a una gran cantidad de conjuntos de datos vectoriales.

La mayoría de los geoportales ofrecen una interfaz gráfica que permite consultar los conjuntos de datos en función de características como la extensión espacial y temporal, siendo el EarthExplorer del Servicio Geológico de Estados Unidos un buen ejemplo. Explorar los conjuntos de datos de forma interactiva en un navegador es una forma eficaz de comprender las capas disponibles. Sin embargo, desde el punto de vista de la reproducibilidad y la eficiencia, es mejor descargar los datos mediante código. Las descargas pueden iniciarse desde la línea de comandos mediante diversas técnicas, principalmente a través de URL y API (véase Sentinel API, por ejemplo). Los archivos alojados en URL estáticas pueden descargarse con download.file(), como se ilustra en el fragmento de código siguiente que accede a PeRL: Permafrost Region Pond and Lake Database de doi.pangaea.de:

download.file(url = "https://hs.pangaea.de/Maps/PeRL/PeRL_permafrost_landscapes.zip",
              destfile = "PeRL_permafrost_landscapes.zip", 
              mode = "wb")
unzip("PeRL_permafrost_landscapes.zip")
canada_perma_land = read_sf("PeRL_permafrost_landscapes/canada_perma_land.shp")

Paquetes de datos geográficos

Se han desarrollado muchos paquetes de R para acceder a datos geográficos, algunos de los cuales se presentan a continuación. Éstos proporcionan interfaces a una o más bibliotecas espaciales o geoportales y tienen como objetivo hacer que el acceso a los datos sea aún más rápido desde la línea de comandos.

Selected R packages for geographic data retrieval.
Package Description
FedData Datasets maintained by the US Federal government, including elevation and land cover.
geodata Download and import imports administrative, elevation, WorldClim data.
osmdata Download and import small OpenStreetMap datasets.
osmextract Download and import large OpenStreetMap datasets.
rnaturalearth Access to Natural Earth vector and raster data.
rnoaa Imports National Oceanic and Atmospheric Administration (NOAA) climate data.

Cabe destacar que la tabla representa sólo un pequeño número de paquetes de datos geográficos disponibles. Por ejemplo, existe un gran número de paquetes R para obtener diversos datos sociodemográficos, como tidycensus y tigris (EE.UU.), cancensus (Canadá), eurostat y giscoR (Unión Europea), o idbr (bases de datos internacionales) – lea Analyzing US Census Data para encontrar algunos ejemplos de cómo analizar dichos datos. Del mismo modo, existen varios paquetes de R que dan acceso a datos espaciales de diversas regiones y países, como bcdata (provincia de Columbia Británica), geobr (Brasil), RCzechia (Chequia) o rgugik (Polonia). Otro paquete destacable es GSODR, que proporciona un resumen global de datos meteorológicos diarios en R (consulte el README del paquete para obtener una visión general de las fuentes de datos meteorológicos).

Cada paquete de datos tiene su propia sintaxis para acceder a los datos. Esta diversidad se demuestra en los siguientes trozos de código, que muestran cómo obtener datos utilizando tres paquetes de la tabla anterior.2. Las fronteras de los países suelen ser útiles y se puede acceder a ellas con la función ne_countries() del paquete rnaturalearth de la siguiente manera:

library(rnaturalearth)
usa = ne_countries(country = "United States of America") # United States borders
class(usa)
# alternative way of accessing the data, with geodata
# geodata::gadm("USA", level = 0, path = tempdir())
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

Por defecto rnaturalearth devuelve objetos de la clase Spatial*. El resultado se puede convertir en un objeto sf con st_as_sf() de la siguiente manera:

usa_sf = st_as_sf(usa)

Un segundo ejemplo descarga una serie de rásters que contienen sumas mensuales globales de precipitación con una resolución espacial de diez minutos (~18,5 km en el ecuador) utilizando el paquete geodata. El resultado es un objeto multicapa de la clase SpatRaster.

library(geodata)
worldclim_prec = worldclim_global("prec", res = 10, path = tempdir())
class(worldclim_prec)

Un tercer ejemplo utiliza el paquete osmdata para buscar parques en la base de datos OpenStreetMap (OSM). Como se ilustra en el fragmento de código siguiente, las consultas comienzan con la función opq() (abreviatura de OpenStreetMap query), cuyo primer argumento es bounding box, o cadena de texto que representa un bounding box (la ciudad de Leeds en este caso). El resultado se pasa a una función para seleccionar los elementos OSM que nos interesan (parques en este caso), representados por pares clave-valor. A continuación, se pasan a la función osmdata_sf() que hace el trabajo de descargar los datos y convertirlos en una lista de objetos sf (ver vignette('osmdata') para más detalles):

library(osmdata)
parks = opq(bbox = "leeds uk") |> 
  add_osm_feature(key = "leisure", value = "park") |> 
  osmdata_sf()

Una limitación del paquete osmdata es que no puede descargar grandes conjuntos de datos OSM (por ejemplo, todos los datos OSM de una gran ciudad). Para superar esta limitación, se desarrolló el paquete osmextract, que se puede utilizar para descargar e importar archivos binarios .pbf que contienen versiones comprimidas de la base de datos OSM para regiones predefinidas.

OpenStreetMap es una vasta base de datos mundial de origen colectivo que crece día a día y cuenta con un amplio ecosistema de herramientas que facilitan el acceso a los datos, desde el servicio web Overpass turbo para desarrollar y probar rápidamente consultas OSM hasta osm2pgsql para importar los datos a una base de datos PostGIS. Aunque la calidad de los conjuntos de datos derivados de OSM varía, la fuente de datos y los ecosistemas OSM más amplios tienen muchas ventajas: proporcionan conjuntos de datos disponibles a nivel mundial, gratuitos y en constante mejora gracias a un ejército de voluntarios. El uso de OSM fomenta la “ciencia ciudadana” y las contribuciones al patrimonio digital (puedes empezar a editar datos que representen una parte del mundo que conozcas bien en www.openstreetmap.org).

A veces, los paquetes vienen con conjuntos de datos incorporados. Se puede acceder a ellos de cuatro maneras: adjuntando el paquete (si el paquete utiliza “carga lenta” como spData), con data(dataset, package = mypackage), haciendo referencia al conjunto de datos con mypackage::dataset, o con system.file(filepath, package = mypackage) para acceder a los archivos de datos sin procesar. El siguiente fragmento de código ilustra las dos últimas opciones utilizando el conjunto de datos world (ya cargado adjuntando su paquete padre con library(spData)):

world2 = spData::world
world3 = read_sf(system.file("shapes/world.gpkg", package = "spData"))

El último ejemplo, system.file("shapes/world.gpkg", package = "spData"), devuelve una ruta al archivo world.gpkg, que se almacena dentro de la carpeta "shapes/" del paquete spData.

Otra forma de obtener información espacial es realizar la geocodificación: transformar la descripción de una ubicación, normalmente una dirección, en sus coordenadas. Esto suele hacerse enviando una consulta a un servicio en línea y obteniendo la ubicación como resultado. Existen muchos servicios de este tipo que difieren en el método de geocodificación utilizado, las limitaciones de uso, los costes o los requisitos de clave API. R dispone de varios paquetes para geocodificación; sin embargo, tidygeocoder parece permitir conectar con el mayor número de servicios de geocodificación con una interfaz consistente. La función principal de tidygeocoder es geocode, que toma un marco de datos con direcciones y añade coordenadas como "lat" y "long". Esta función también permite seleccionar un servicio de geocodificación con el argumento method y tiene muchos parámetros adicionales.

El siguiente ejemplo busca las coordenadas de la placa azul de John Snow situada en un edificio del barrio londinense del Soho.

library(tidygeocoder)
geo_df = data.frame(address = "54 Frith St, London W1D 4SJ, UK")
geo_df = geocode(geo_df, address, method = "osm")
geo_df

El data.frame resultante puede convertirse en un objeto sf con st_as_sf().

geo_sf = st_as_sf(geo_df, coords = c("long", "lat"), crs = "EPSG:4326")

Este paquete también permite realizar el proceso inverso denominado geocodificación inversa utilizado para obtener un conjunto de información (nombre, dirección, etc.) a partir de un par de coordenadas.

Servicios web geográficos

En un esfuerzo por estandarizar las API web para acceder a los datos espaciales, el Open Geospatial Consortium (OGC) ha creado una serie de especificaciones para los servicios web (conocidos colectivamente como OWS, que es la abreviatura de OGC Web Services). Estas especificaciones incluyen el Web Feature Service (WFS), el Web Map Service (WMS), el Web Map Tile Service (WMTS), el Web Coverage Service (WCS) e incluso un Web Processing Service (WPS). Servidores de mapas como PostGIS han adoptado estos protocolos, lo que ha llevado a la estandarización de las consultas. Al igual que otras API web, las API OWS utilizan una “URL base”, un “endpoint” y “argumentos de consulta URL” seguidos de un ? para solicitar datos (véase la viñeta best-practices-api-packages en el paquete httr).

Hay muchas peticiones que se pueden hacer a un servicio OWS. Una de las más fundamentales es getCapabilities, demostrada con las funciones httr GET() y modify_url() a continuación.

El fragmento de código muestra cómo pueden construirse y enviarse consultas API, en este caso para descubrir las capacidades de un servicio gestionado por la Organización de las Naciones Unidas para la Agricultura y la Alimentación (FAO):

library(httr)
base_url = "http://www.fao.org"
endpoint = "/figis/geoserver/wfs"
q = list(request = "GetCapabilities")
res = GET(url = modify_url(base_url, path = endpoint), query = q)
res$url
## [1] "https://www.fao.org/figis/geoserver/wfs?request=GetCapabilities"

El fragmento de código anterior demuestra cómo se pueden construir programáticamente las peticiones API con la función GET(), que toma una URL base y una lista de parámetros de consulta que pueden ampliarse fácilmente. El resultado de la petición se guarda en res, un objeto de la clase response definido en el paquete httr, que es una lista que contiene información de la petición, incluyendo la URL. Como puede verse ejecutando browseURL(res$url), los resultados también pueden leerse directamente en un navegador. Una forma de extraer el contenido de la petición es la siguiente:

txt = content(res, "text")
xml = xml2::read_xml(txt)
xml
#> {xml_document} ...
#> [1] <ows:ServiceIdentification>\n  <ows:Title>GeoServer WFS...
#> [2] <ows:ServiceProvider>\n  <ows:ProviderName>UN-FAO Fishe...
#> ...

Los datos pueden descargarse de los servicios WFS con la petición GetFeature y un typeName específico (como se ilustra en el fragmento de código siguiente).

Los nombres disponibles varían en función del servicio de características web al que se acceda. Se pueden extraer mediante programación utilizando tecnologías web o desplazándose manualmente por el contenido de la salida GetCapabilities en un navegador.

qf = list(request = "GetFeature", typeName = "area:FAO_AREAS")
file = tempfile(fileext = ".gml")
GET(url = base_url, path = endpoint, query = qf, write_disk(file))
fao_areas = read_sf(file)

Nótese el uso de write_disk() para asegurar que los resultados se escriben en disco en lugar de cargarse en memoria, lo que permite importarlos con sf. Este ejemplo muestra cómo obtener acceso de bajo nivel a servicios web utilizando httr, que puede ser útil para entender cómo funcionan los servicios web. Para muchas tareas cotidianas, sin embargo, puede ser más apropiada una interfaz de más alto nivel, y se han desarrollado una serie de paquetes de R, y tutoriales, precisamente con este propósito. El paquete ows4R fue desarrollado para trabajar con servicios OWS.

Formatos de archivo

Los conjuntos de datos geográficos suelen almacenarse como archivos o en bases de datos espaciales. Los formatos de archivo pueden almacenar datos vectoriales o ráster, mientras que las bases de datos espaciales como PostGIS pueden almacenar ambos. Hoy en día, la variedad de formatos de archivo puede parecer desconcertante, pero se ha producido una gran consolidación y estandarización desde los inicios del software SIG en la década de 1960, cuando se creó el primer programa de amplia distribución (SYMAP) para el análisis espacial en la Universidad de Harvard.

GDAL (que debería pronunciarse “goo-dal”, con la doble “o” en referencia a la orientación a objetos), la Biblioteca de Abstracción de Datos Geoespaciales, ha resuelto muchos problemas relacionados con la incompatibilidad entre formatos de archivos geográficos desde su lanzamiento en 2000. GDAL proporciona una interfaz unificada y de alto rendimiento para la lectura y escritura de muchos formatos de datos ráster y vectoriales.3. Muchos programas SIG abiertos y propietarios, incluidos GRASS GIS, ArcGIS y QGIS, utilizan GDAL detrás de sus interfaces gráficas de usuario para realizar el trabajo de ingesta y escupir los datos geográficos en los formatos adecuados.

GDAL proporciona acceso a más de 200 formatos de datos vectoriales y ráster. La siguiente tabla presenta información básica sobre los formatos de archivos espaciales seleccionados y utilizados con frecuencia.

Selected spatial file formats.
Name Extension Info Type Model
ESRI Shapefile .shp (the main file) Popular format consisting of at least three files. No support for: files > 2GB; mixed types; names > 10 chars; cols > 255. Vector Partially open
GeoJSON .geojson Extends the JSON exchange format by including a subset of the simple feature representation; mostly used for storing coordinates in longitude and latitude; it is extended by the TopoJSON format Vector Open
KML .kml XML-based format for spatial visualization, developed for use with Google Earth. Zipped KML file forms the KMZ format. Vector Open
GPX .gpx XML schema created for exchange of GPS data. Vector Open
FlatGeobuf .fgb Single file format allowing for quick reading and writing of vector data. Has streaming capabilities. Vector Open
GeoTIFF .tif/.tiff Popular raster format. A TIFF file containing additional spatial metadata. Raster Open
Arc ASCII .asc Text format where the first six lines represent the raster header, followed by the raster cell values arranged in rows and columns. Raster Open
SQLite/SpatiaLite .sqlite Standalone relational database, SpatiaLite is the spatial extension of SQLite. Vector and raster Open
ESRI FileGDB .gdb Spatial and nonspatial objects created by ArcGIS. Allows: multiple feature classes; topology. Limited support from GDAL. Vector and raster Proprietary
GeoPackage .gpkg Lightweight database container based on SQLite allowing an easy and platform-independent exchange of geodata Vector and (very limited) raster Open

La fundación del Open Geospatial Consortium (OGC) en 1994 supuso un avance importante en la normalización y la apertura de los formatos de archivo. Además de definir el modelo de datos de características simples o simple features, el OGC también coordina el desarrollo de normas abiertas, por ejemplo, las utilizadas en formatos de archivo como KML y GeoPackage. Los formatos de archivo abiertos del tipo respaldado por el OGC tienen varias ventajas sobre los formatos propietarios: los estándares se publican, garantizan la transparencia y abren la posibilidad de que los usuarios sigan desarrollando y ajustando los formatos de archivo a sus necesidades específicas.

ESRI Shapefile es el formato de intercambio de datos vectoriales más popular; sin embargo, no es un formato abierto (aunque su especificación sí lo es). Se desarrolló a principios de la década de 1990 y tiene una serie de limitaciones. En primer lugar, es un formato multiarchivo, que consta de al menos tres archivos. Sólo admite 255 columnas, los nombres de las columnas están restringidos a diez caracteres y el límite de tamaño de los archivos es de 2 GB. Además, ESRI Shapefile no admite todos los tipos de geometría posibles, por ejemplo, no es capaz de distinguir entre un polígono y un multipolígono.4 A pesar de estas limitaciones, hacía tiempo que se echaba en falta una alternativa viable. Mientras tanto, surgió GeoPackage, y parece ser un candidato más que adecuado para sustituir a ESRI Shapefile. Geopackage es un formato de intercambio de información geoespacial y una norma OGC. La norma GeoPackage describe las reglas para almacenar información geoespacial en un pequeño contenedor SQLite. Por lo tanto, GeoPackage es un contenedor de base de datos espacial ligero, que permite almacenar datos vectoriales y ráster, pero también datos no espaciales y extensiones. Aparte de GeoPackage, hay otros formatos de intercambio de datos geoespaciales que merece la pena consultar.

El formato GeoTIFF parece ser el formato de datos ráster más destacado. Permite incrustar información espacial, como CRS, en un archivo TIFF. Al igual que ESRI Shapefile, este formato se desarrolló por primera vez en los años 90, pero como formato abierto. Además, GeoTIFF sigue ampliándose y mejorándose. Una de las recientes incorporaciones más significativas al formato GeoTIFF es su variante denominada COG (Cloud Optimized GeoTIFF). Los objetos ráster guardados como COG pueden alojarse en servidores HTTP, de modo que otras personas puedan leer sólo partes del archivo sin descargarlo entero.

Además, se están desarrollando nuevos formatos de datos espaciales (por ejemplo, GeoParquet y Zarr) y se están mejorando los existentes. Si necesita más información sobre otros formatos, le recomendamos que lea la documentación de GDAL sobre los controladores vector y raster. Además, algunos formatos de datos espaciales pueden almacenar otros modelos (tipos) de datos que no sean vectoriales o ráster. Entre ellos se incluyen los formatos LAS y LAZ para almacenar nubes de puntos lidar, y NetCDF y HDF para almacenar matrices multidimensionales.

Los datos espaciales también suelen almacenarse en formatos de texto tabulares (no espaciales), como archivos CSV u hojas de cálculo Excel. Por ejemplo, esto puede ser conveniente para compartir muestras espaciales con personas que no utilizan herramientas SIG o intercambiar datos con otro software que no acepte formatos de datos espaciales. Sin embargo, este enfoque tiene varios posibles problemas: es bastante complicado para almacenar geometrías más complejas que los POINT y no almacena directamente información sobre los CRS.

Entrada de datos (I)

La ejecución de comandos como sf::read_sf() (la función principal que utilizamos para cargar datos vectoriales) o terra::rast() (la función principal utilizada para cargar datos raster) desencadena silenciosamente una cadena de eventos que lee datos de archivos. Además, existen muchos paquetes de R que contienen una amplia gama de datos geográficos o que proporcionan un acceso sencillo a distintas fuentes de datos. Todos ellos cargan los datos en R o, más exactamente, asignan objetos a su espacio de trabajo, almacenados en RAM accesible desde el .GlobalEnv de la sesión de R.

Vector data

Los datos vectoriales espaciales vienen en una amplia variedad de formatos de archivo. Las representaciones más populares, como los archivos .geojson y .gpkg se pueden importar directamente en R con la función sf read_sf() (o la equivalente st_read()), que utiliza los controladores vectoriales de GDAL entre bastidores. st_drivers() devuelve un marco de datos que contiene name y long_name en las dos primeras columnas, y las características de cada controlador disponible para GDAL (y por lo tanto sf), incluyendo la capacidad de escribir datos y almacenar datos raster en las columnas siguientes, como se ilustra para los formatos de archivo clave en la Tabla (tab:drivers).
Los siguientes comandos muestran los tres primeros controladores de los que informa la instalación de GDAL del ordenador (los resultados pueden variar en función de la versión de GDAL instalada) y un resumen de sus características. Nótese que la mayoría de los drivers pueden escribir datos (51 de 87) mientras que sólo 16 formatos pueden representar eficientemente datos raster además de datos vectoriales (ver ?st_drivers() para más detalles):

sf_drivers = st_drivers()
head(sf_drivers, n = 3)
summary(sf_drivers[-c(1:2)])
Popular drivers/formats for reading/writing vector data.
name long_name write copy is_raster is_vector vsi
ESRI Shapefile ESRI Shapefile TRUE FALSE FALSE TRUE TRUE
GPX GPX TRUE FALSE FALSE TRUE TRUE
KML Keyhole Markup Language (KML) TRUE FALSE FALSE TRUE TRUE
GeoJSON GeoJSON TRUE FALSE FALSE TRUE TRUE
GPKG GeoPackage TRUE TRUE TRUE TRUE TRUE
FlatGeobuf FlatGeobuf TRUE FALSE FALSE TRUE TRUE

El primer argumento de read_sf() es dsn, que debe ser una cadena de texto o un objeto que contenga una única cadena de texto. El contenido de una cadena de texto puede variar entre diferentes controladores. En la mayoría de los casos, como con el ESRI Shapefile (.shp) o el formato GeoPackageindex{GeoPackage} (.gpkg), el dsn sería un nombre de archivo. read_sf() adivina el controlador basado en la extensión del archivo, como se ilustra para un archivo .gpkg a continuación:

f = system.file("shapes/world.gpkg", package = "spData")
world = read_sf(f, quiet = TRUE)

Para algunos controladores, dsn podría proporcionarse como un nombre de carpeta, credenciales de acceso a una base de datos, o una representación de cadena GeoJSON (ver los ejemplos de la página de ayuda read_sf() para más detalles).

Algunos formatos de controladores vectoriales pueden almacenar múltiples capas de datos. Por defecto, read_sf() lee automáticamente la primera capa del archivo especificado en dsn; sin embargo, utilizando el argumento layer se puede especificar cualquier otra capa.

La función read_sf() también permite leer sólo partes del fichero en la RAM con dos mecanismos posibles. El primero está relacionado con el argumento query, que permite especificar qué parte de los datos leer con el texto de consulta OGR SQL. En el ejemplo siguiente se extraen los datos de Tanzania únicamente. Se hace especificando que queremos obtener todas las columnas (SELECT *) de la capa "world" para las que el name_long es igual a "Tanzania":

tanzania = read_sf(f, query = 'SELECT * FROM world WHERE name_long = "Tanzania"')

Si no conoce los nombres de las columnas disponibles, un buen método consiste en leer simplemente una fila de los datos con 'SELECT * FROM world WHERE FID = 1'. El FID representa un ID de característica – la mayoría de las veces, es un número de fila; sin embargo, sus valores dependen del formato de fichero utilizado. Por ejemplo, FID empieza por 0 en ESRI Shapefile, por 1 en algunos otros formatos de archivo, o incluso puede ser arbitrario.

El segundo mecanismo utiliza el argumento wkt_filter. Este argumento espera un texto conocido que representa el área de estudio de la que queremos extraer los datos. Intentémoslo con un pequeño ejemplo: queremos leer los polígonos de nuestro archivo que se cruzan con el buffer de 50.000 metros de las fronteras de Tanzania. Para ello, tenemos que preparar nuestro “filtro” (a) creando el buffer, (b) convirtiendo el objeto buffer sf en un objeto geométrico sfc con st_geometry(), y (c) traduciendo las geometrías a su conocida representación de texto con st_as_text():

tanzania_buf = st_buffer(tanzania, 50000)
tanzania_buf_geom = st_geometry(tanzania_buf)
tanzania_buf_wkt = st_as_text(tanzania_buf_geom)

Ahora, podemos aplicar este “filtro” usando el argumento wkt_filter.

tanzania_neigh = read_sf(f, wkt_filter = tanzania_buf_wkt)

Nuestro resultado contiene Tanzania y todos los países que se encuentran dentro de su buffer de 50 km.

Reading a subset of the vector data using a query (A) and a wkt filter (B).

Naturalmente, algunas opciones son específicas de determinados controladores.5 Por ejemplo, pensemos en coordenadas almacenadas en un formato de hoja de cálculo (.csv). Para leer estos archivos como objetos espaciales, naturalmente tenemos que especificar los nombres de las columnas (X y Y en nuestro ejemplo de abajo) que representan las coordenadas. Podemos hacerlo con la ayuda del parámetro options. Para conocer las posibles opciones, consulte la sección Opciones abiertas de la descripción del controlador GDAL correspondiente. Para el formato de valores separados por comas (csv), visite http://www.gdal.org/drv_csv.html.

cycle_hire_txt = system.file("misc/cycle_hire_xy.csv", package = "spData")
cycle_hire_xy = read_sf(cycle_hire_txt,
  options = c("X_POSSIBLE_NAMES=X", "Y_POSSIBLE_NAMES=Y"))

En lugar de columnas que describan las coordenadas “XY”, una sola columna puede contener también la información geométrica. Los formatos de texto conocido (WKT), binario conocido (WKB) y GeoJSON son ejemplos de ello. Por ejemplo, el archivo world_wkt.csv tiene una columna llamada WKT que representa polígonos de los países del mundo. Una vez más, utilizaremos el parámetro options para indicarlo.

world_txt = system.file("misc/world_wkt.csv", package = "spData")
world_wkt = read_sf(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT")
# the same as
world_wkt2 = st_read(world_txt, options = "GEOM_POSSIBLE_NAMES=WKT", 
                    quiet = TRUE, stringsAsFactors = FALSE, as_tibble = TRUE)
Nota

No todos los formatos de ficheros vectoriales soportados almacenan información sobre su sistema de referencia de coordenadas. En estas situaciones, es posible añadir la información que falta utilizando la función st_set_crs().

Como ejemplo final, mostraremos cómo read_sf() también lee archivos KML. Un archivo KML almacena información geográfica en formato XML, un formato de datos para la creación de páginas web y la transferencia de datos de forma independiente de la aplicación. Aquí, accedemos a un archivo KML desde la web. Este archivo contiene más de una capa. st_layers() lista todas las capas disponibles. Elegimos la primera capa Placemarks y lo decimos con la ayuda del parámetro layer en read_sf().

u = "https://developers.google.com/kml/documentation/KML_Samples.kml"
download.file(u, "KML_Samples.kml")
st_layers("KML_Samples.kml")
kml = read_sf("KML_Samples.kml", layer = "Placemarks")
## Driver: KML 
## Available layers:
##              layer_name  geometry_type features fields crs_name
## 1            Placemarks       3D Point        3      2   WGS 84
## 2      Highlighted Icon       3D Point        1      2   WGS 84
## 3                 Paths 3D Line String        6      2   WGS 84
## 4         Google Campus     3D Polygon        4      2   WGS 84
## 5      Extruded Polygon     3D Polygon        1      2   WGS 84
## 6 Absolute and Relative     3D Polygon        4      2   WGS 84

Todos los ejemplos presentados hasta ahora en esta sección han utilizado el paquete sf para la importación de datos geográficos. Es rápido y flexible, pero puede merecer la pena buscar otros paquetes para formatos de archivo específicos. Un ejemplo es el paquete geojsonsf. Un benchmark sugiere que es unas 10 veces más rápido que el paquete sf para leer .geojson.

Datos raster

Al igual que los datos vectoriales, los datos raster vienen en muchos formatos de archivo, algunos de los cuales soportan archivos multicapa. El comando rast() de terra lee una sola capa cuando se proporciona un archivo con una sola capa.

raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
single_layer = rast(raster_filepath)

También funciona en caso de que quieras leer un fichero multicapa.

multilayer_filepath = system.file("raster/landsat.tif", package = "spDataLarge")
multilayer_rast = rast(multilayer_filepath)

Todos los ejemplos anteriores leen información espacial de archivos almacenados en su disco duro. Sin embargo, GDAL también permite leer datos directamente de recursos en línea, como recursos web HTTP/HTTPS/FTP. Lo único que tenemos que hacer es añadir un prefijo /vsicurl/ antes de la ruta al fichero. Vamos a probarlo conectándonos a la probabilidad de nieve mensual global a 500 m de resolución para el periodo 2000-2012. La probabilidad de nieve para diciembre se almacena como archivo COG (Cloud Optimized GeoTIFF) en . Para leer un archivo en línea, basta con proporcionar su URL junto con el prefijo /vsicurl/.

myurl = "/vsicurl/https://zenodo.org/record/5774954/files/clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0.tif"
snow = rast(myurl)
snow
## class       : SpatRaster 
## dimensions  : 35849, 86400, 1  (nrow, ncol, nlyr)
## resolution  : 0.004166667, 0.004166667  (x, y)
## extent      : -180, 180, -62.00083, 87.37  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0.tif 
## name        : clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0

Debido a que los datos de entrada son COG, en realidad no estamos leyendo este fichero en nuestra RAM, sino creando una conexión con él sin obtener ningún valor. Sus valores serán leídos si aplicamos cualquier operación basada en valores (por ejemplo, crop() o extract()). Esto nos permite también leer sólo una pequeña porción de los datos sin descargar todo el fichero. Por ejemplo, podemos obtener la probabilidad de nieve para diciembre en Reikiavik (70%) especificando sus coordenadas y aplicando la función extract():

rey = data.frame(lon = -21.94, lat = 64.15)
snow_rey = extract(snow, rey)
snow_rey
ID clm_snow.prob_esacci.dec_p.90_500m_s0..0cm_2000..2012_v2.0
1 70

De esta forma, sólo descargamos un único valor en lugar de todo el archivo GeoTIFF de gran tamaño.

El ejemplo anterior sólo muestra un caso sencillo (pero útil), pero hay más por explorar. El prefijo /vsicurl/ también funciona no sólo para los formatos de archivo raster, sino también para los vectoriales. Permite leer vectores directamente desde el almacenamiento en línea con read_sf() simplemente añadiendo el prefijo antes de la URL del archivo vectorial.

Es importante destacar que /vsicurl/ no es el único prefijo proporcionado por GDAL; existen muchos más, como /vsizip/ para leer archivos espaciales de archivos ZIP sin descomprimirlos previamente o /vsis3/ para leer sobre la marcha archivos disponibles en cubos S3 de AWS. Puede obtener más información en https://gdal.org/user/virtual_file_systems.html.

Salida de datos (O)

La escritura de datos geográficos permite convertir de un formato a otro y guardar objetos recién creados. Dependiendo del tipo de datos (vector o raster), clase de objeto (por ejemplo, sf o SpatRaster), y tipo y cantidad de información almacenada (por ejemplo, tamaño del objeto, rango de valores), es importante saber cómo almacenar archivos espaciales de la manera más eficiente. En las dos secciones siguientes se muestra cómo hacerlo.

Vector data

La contrapartida de read_sf() es write_sf(). Permite escribir objetos sf en una amplia gama de formatos de archivos vectoriales geográficos, incluidos los más comunes como .geojson, .shp y .gpkg. Basándose en el nombre del archivo, write_sf() decide automáticamente qué controlador utilizar. La velocidad del proceso de escritura depende también del controlador.

write_sf(obj = world, dsn = "world.gpkg")

Nota: si intentas escribir de nuevo en la misma fuente de datos, la función sobrescribirá el fichero:

write_sf(obj = world, dsn = "world.gpkg")

En lugar de sobrescribir el archivo, podríamos añadir una nueva capa al archivo con append = TRUE, que es soportado por varios formatos espaciales, incluyendo GeoPackage.

write_sf(obj = world, dsn = "world_many_layers.gpkg", append = TRUE)

Alternativamente, puede usar st_write() ya que es equivalente a write_sf(). Sin embargo, tiene diferentes valores por defecto – no sobrescribe archivos (devuelve un error cuando intenta hacerlo) y muestra un breve resumen del formato del archivo escrito y el objeto.

st_write(obj = world, dsn = "world2.gpkg")
## Writing layer `world2' to data source `world2.gpkg' using driver `GPKG'
## Writing 177 features with 10 fields and geometry type Multi Polygon.

El argumento layer_options también se puede utilizar para muchos propósitos diferentes. Uno de ellos es escribir datos espaciales en un fichero de texto. Esto puede hacerse especificando GEOMETRY dentro de layer_options. Puede ser AS_XY para conjuntos de datos de puntos simples (crea dos nuevas columnas para las coordenadas) o AS_WKT para datos espaciales más complejos (se crea una nueva columna que contiene la conocida representación textual de los objetos espaciales).

write_sf(cycle_hire_xy, "cycle_hire_xy.csv", layer_options = "GEOMETRY=AS_XY")
write_sf(world_wkt, "world_wkt.csv", layer_options = "GEOMETRY=AS_WKT")

Datos Raster

La función writeRaster() guarda objetos SpatRaster en archivos de disco. La función espera información sobre el tipo de datos de salida y el formato de archivo, pero también acepta opciones de GDAL específicas para el formato de archivo seleccionado (véase ?writeRaster para más detalles).

El paquete terra ofrece siete tipos de datos a la hora de guardar una trama: INT1U, INT2S, INT2U, INT4S, INT4U, FLT4S, y FLT8S.6 El tipo de datos determina la representación en bits del objeto ráster escrito en el disco. El tipo de datos a utilizar depende del rango de valores del objeto raster. Cuantos más valores pueda representar un tipo de datos, más grande será el archivo en disco. Los enteros sin signo (INT1U, INT2U, INT4U) son adecuados para datos categóricos, mientras que los números flotantes (FLT4S y FLT8S) suelen representar datos continuos. writeRaster() utiliza FLT4S por defecto. Aunque esto funciona en la mayoría de los casos, el tamaño del fichero de salida será innecesariamente grande si guarda datos binarios o categóricos. Por lo tanto, recomendamos utilizar el tipo de datos que necesite menos espacio de almacenamiento, pero que sea capaz de representar todos los valores (compruebe el rango de valores con la función summary()).

Data types supported by the terra package.
Data type Minimum value Maximum value
INT1U 0 255
INT2S -32,767 32,767
INT2U 0 65,534
INT4S -2,147,483,647 2,147,483,647
INT4U 0 4,294,967,296
FLT4S -3.4e+38 3.4e+38
FLT8S -1.7e+308 1.7e+308

Por defecto, el formato del archivo de salida se deriva del nombre del archivo. Al nombrar un archivo *.tif se creará un archivo GeoTIFF, como se muestra a continuación:

writeRaster(single_layer, filename = "my_raster.tif", datatype = "INT2U")

Algunos formatos de archivo ráster tienen opciones adicionales, que pueden establecerse proporcionando parámetros GDAL al argumento options de writeRaster(). Los archivos GeoTIFF se escriben en terra, por defecto, con la compresión LZW gdal = c("COMPRESS=LZW"). Para cambiar o desactivar la compresión, necesitamos modificar este argumento.

writeRaster(x = single_layer, filename = "my_raster.tif",
            gdal = c("COMPRESS=NONE"), overwrite = TRUE)

Además, podemos guardar nuestro objeto raster como COG (Cloud Optimized GeoTIFF) con las opciones filetype = "COG".

writeRaster(x = single_layer, filename = "my_raster.tif",
            filetype = "COG", overwrite = TRUE)

Resultados visuales

R soporta muchos formatos diferentes de gráficos estáticos e interactivos. El método más general para guardar un gráfico estático es abrir un dispositivo gráfico, crear un gráfico y cerrarlo, por ejemplo:

png(filename = "lifeExp.png", width = 500, height = 350)
plot( world["lifeExp"])
dev.off()

Otros dispositivos gráficos disponibles son pdf(), bmp(), jpeg() y tiff(). Puede especificar varias propiedades del gráfico de salida, como la anchura, la altura y la resolución.

Además, varios paquetes gráficos proporcionan sus propias funciones para guardar una salida gráfica. Por ejemplo, el paquete tmap tiene la función tmap_save(). Puede guardar un objeto tmap en diferentes formatos gráficos o en un archivo HTML especificando el nombre del objeto y una ruta a un nuevo archivo.

library(tmap)
tmap_obj = tm_shape(world) + tm_polygons(col = "lifeExp")
tmap_save(tmap_obj, filename = "lifeExp_tmap.png")

Por otra parte, puede guardar los mapas interactivos creados en el paquete mapview como un archivo HTML o una imagen utilizando la función mapshot():

library(mapview)
mapview_obj = mapview(world, zcol = "lifeExp", legend = TRUE)
mapshot(mapview_obj, file = "my_interactive_map.html")

Notas

  1. Por ejemplo, visite https://freegisdata.rtwilson.com/ para consultar una lista de sitios web con conjuntos de datos geográficos de libre acceso↩︎

  2. En https://rspatialdata.github.io/ se pueden encontrar más ejemplos de descarga de datos utilizando paquetes R dedicados↩︎

  3. GDAL también contiene un conjunto de funciones de utilidad que permiten la creación de mosaicos ráster, el remuestreo, el recorte y la reproyección, etc.↩︎

  4. Para obtener más información sobre las limitaciones de ESRI Shapefile y posibles formatos de archivo alternativos, visite http://switchfromshapefile.org/.↩︎

  5. Puede encontrar una lista de los formatos vectoriales y opciones compatibles en http://gdal.org/ogr_formats.html.↩︎

  6. No se recomienda usar INT4U ya que R no soporta enteros sin signo de 32 bits.↩︎