3 Representación de capas

En el resto del documento, basándonos en las operaciones básicas presentadas en los ejemplos anteriores, vamos a usar R como un SIG.

Vídeo: https://youtu.be/0WkOijfG8BQ

3.1 Obtención de los datos

Usaremos los mismos datos que hemos usado con QGIS.

En QGIS, cuando se añaden nuevos campos, se modifican las tablas de atributos de las capas en los archivos de estas. Por este motivo, aunque se hayan obtenido los datos para QGIS, es recomendable usar otra carpeta de trabajo y obtener de nuevo las capas. En el caso del archivo Excel descargado del INE, podemos generar de nuevo la tabla de datos o copiar el archivo de la tabla y eliminar la columna adicional añadida en QGIS.

3.2 Leer capas

En el ejemplo siguiente, comenzamos cargando el paquete sf. Hay muchos paquetes para trabajar con información geográfica. El paquete sp se usaba frecuentemente, sf está ocupando su lugar. Muchos paquetes, inicialmente basados en sp, se han adaptado a sf. Todavía quedan paquetes basados en sp que es posible que nos interese usar: podemos usarlos convirtiendo los datos de la estructura sf a sp.

library(sf)

censo <- st_read("datos/Granada/municipios censo GR.shp")
class(censo)
summary(censo)

hidro <- st_read("datos/Granada/red hidrografica GR.shp")
summary(hidro)

censo_sp <- as(censo, Class = "Spatial")
# funciones de sp.
class(censo_sp)
sp::summary(censo_sp)

censo_sf <- st_as_sf(censo_sp)
class(censo_sf)
summary(censo_sf)

Mediante la función st_read() del paquete sf leemos una capa de la carpeta de trabajo. La función class() nos muestra que hemos leído un data frame. La función summary() nos muestra un resumen de las columnas del data frame. Las columnas son las mismas que muestra QGIS en la tabla de atributos pero, además, tenemos la columna geometry que es la que contiene los datos geográficos de la capa. Para representar solo la geometría de una capa, esta es la columna que debemos utilizar explícitamente.

3.3 Representar capas

La función st_geometry() del paquete sf nos devuelve el componente geometry de una capa¹¹. Representamos la capa mediante la función plot(). Al resultado podemos añadir otras capas mediante la misma función, pero utilizando el parámetro add para indicarlo. El resultado se muestra en la zona inferior derecha de la ventana, en la pestaña Plots. Podemos exportarlo a distintos formatos mediante la opción Export.

plot(st_geometry(censo))
plot(st_geometry(hidro), col = "blue", add = TRUE)

3.4 Cambiar el color y otros elementos de una capa

Para cambiar la representación de una capa debemos representarla de nuevo. La función plot() ofrece muchas posibilidades, algunas de ellas se muestran a continuación: hemos cambiado el color, el grosor de las líneas y añadido ejes a la representación.

plot(st_geometry(censo), col = "lightgreen", lwd = 1.5, axes = TRUE)

3.5 Importar una imagen desde la Web

Vamos a importar una imagen desde OpenStreetMap¹². Hay un paquete R llamado OpenStreetMap que ofrece esta funcionalidad.

library(OpenStreetMap)

censo_cg <- st_transform(censo, 4326)
summary(censo_cg)

# rectángulo más pequeño que incluye la capa
sup_izq <- as.vector(cbind(st_bbox(censo_cg)['ymax'],
                           st_bbox(censo_cg)['xmin']))
inf_der <- as.vector(cbind(st_bbox(censo_cg)['ymin'],
                           st_bbox(censo_cg)['xmax']))

mapa_osm <- openmap(sup_izq, inf_der, type = "bing")
censo_osm <- st_transform(censo_cg, osm())

# restaurar el área de visualización
graphics.off()

plot(mapa_osm)
plot(st_geometry(censo_osm), add = TRUE, lwd = 1.2)

La función de OpenStreetMap que obtiene un mapa de una zona determinada es openmap(). Si miramos su definición, trabaja con la latitud y longitud. Deberemos transformar la capa que queramos representar a un CRS que permita obtener directamente estos datos: por ejemplo, el de código EPSG 4326. La transformación se puede llevar a cabo mediante la función st_transform().

A partir de esta representación de la capa, obtenemos los vértices del rectángulo más pequeño que la incluye. Con dos de esos vértices, mediante la función openmap() obtenemos el mapa de la Web, el tipo de mapa a importar se define mediante el parámetro type. Como vamos a representar el mapa obtenido conjuntamente con la capa de partida, vamos a transformar la capa para que tenga el mismo CRS que OpenStreetMap mediante las funciones st_transform(), que realiza la transformación, y osm(), que devuelve el CRS de OpenStreetMap.

Figura 3.1: Ejemplo con una de las capas y una imagen de OpenStreetMap.

En la figura 3.1, se muestra el resultado obtenido al presentar una de las capas con la capa obtenida de OpenStreetMap mediante el código de los ejemplos.

Es recomendable usar esta función en lugar de acceder directamente al componente geometry del data frame porque el resultado obtenido con ambos métodos de acceso no siempre es idéntico.↩︎
https://www.openstreetmap.org ↩︎