Lectura, filtrado y visualización de datos vectoriales geoespaciales en R

Luisa Fernanda Carrión Ramírez 24/03/2020

Este cuaderno tiene como objetivo mostrar la utilidad de algunas funciones de R para el manejo de datos geoespaciales. Para ello, se deben instalar dos bibliotecas denominadas sf (simple features) y tidyverse. Se empezará con instalar la biblioteca “sf” con el siguiente comando:

#install.packages("sf")

Posteriormente, se instalará el paquete tidyverse de la misma manera.

#install.packages("tidyverse")

Luego, se deben cargar ambas bibliotecas empezando con sf:

library(sf)

Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

Se aplicará el mismo comando para la biblioteca tidyverse:

library(tidyverse)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.0     v purrr   0.3.3
v tibble  2.1.3     v dplyr   0.8.4
v tidyr   1.0.2     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

Lectura de datos vectoriales

En este apartado, se quiere leer un archivo shapefile previamente descargado en la página DIVA-GIS que contiene la división politica departamental de Colombia. Para ello, se emplea el siguiente comando:

deptos <- read_sf("C:/Users/LUISA CARRION/Documents/Geomática/COL_adm1.shp")

Para saber la información que el shapefile relacionada a los departamentos se procede a escribir el siguiente comando:

head(deptos)

Simple feature collection with 6 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -77.149 ymin: -4.228429 xmax: -69.36835 ymax: 11.10792
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

En la biblioteca sf existe una función que permite saber cuál es el sistema de referencia de coordenadas de los datos vectoriales almacenados. En este caso, se quiere conocer esta información en el objeto deptos aplicando el siguiente comando:

st_crs(deptos)

Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Usando ggplot para visualizar datos geoespaciales

R nos permite gráficar los datos contenidos en un shapefile mediante la función ggplot():

ggplot()+ geom_sf(data = deptos)

Cualquier sistema de referencia de coordenadas puede ser empleado desde que este diseñado para se aplicado en la región geográfica deseada. En la página epsg.io/ se puede obtner información a cerca de que sistemas de referencia de coordenadas pueden usarse en cada país y explica brevemente en que consite.

ggplot() + geom_sf(data = deptos) + coord_sf(crs=st_crs(3978))

En el gráfico anterior se usó el CRS 3978 que es un sistema diseñado para Canadá y algunas zonas de Estados Unidos por lo cual, para trabajar en Colombia se requiere otro CRS. Al buscar en la página anteriormente mencionada, encontramos que un sistema referenciado de coordenadas referenciadas apto para este caso es el 32618. Para hacer el cambio de CRS se emplean los siguientes comandos:

deptos_utm <- st_transform(deptos, crs = st_crs(32618))

deptos_utm

Simple feature collection with 32 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -245935.3 ymin: -469204.3 xmax: 1407491 ymax: 1763314
epsg (SRID):    32618
proj4string:    +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

ggplot() + geom_sf(data = deptos_utm)

Filtrar datos geoespaciales basados en atributos

A lo largo del semestre cada alumno investigará sobre un departamento colombiano de su elección. En mi caso en particular, escogí Casanare por lo cual, debo filtrar este departamento mediante el siguiente comando:

casanare <- deptos %>%  filter(NAME_1 == "Casanare")

Para obtener la gráfica se debe escribir:

ggplot() + geom_sf(data = casanare)

Como se requiere más información para completar el mapa de Casanare, se debe importar el archivo COL_adm2.shp. Luego, se filtra y por último se usa la función ggplot () para graficar los datos:

munic <- read_sf("C:/Users/LUISA CARRION/Documents/Geomática/COL_adm2.shp")
mun_casanare <- munic %>% filter(NAME_1 == "Casanare")
ggplot() + geom_sf(data = mun_casanare)

Si se usa el comando mun_casanare, se puede observar una tabla que nos muestra los municipios pertenecientes a Casanare:

mun_casanare

Simple feature collection with 19 features and 11 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -73.0989 ymin: 4.246699 xmax: -69.84787 ymax: 6.250501
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

A continuación, se quiere ponen una numeración a los municipios del departamento. Se intentó realizar este proceso con dos comandos que arrojaron error:

casanare_points <- st_centroid(mun_casanare)

st_centroid assumes attributes are constant over geometries of xst_centroid does not give correct centroids for longitude/latitude data

casanare_points <- cbind(mun_casanare, st_coordinates(st_centroid(mun_casanare$geometry)))

st_centroid does not give correct centroids for longitude/latitude data

Teniendo esto en cuenta, se optó por escribir un comando más extenso y especícico que se puede observar a continuación:

ggplot(casanare) +
    geom_sf() +
    geom_sf(data = casanare_points, fill = "mistyrose2") + 
    geom_text(data = casanare_points, aes(x=X, y=Y,label = ID_2), size = 2.3) +
    coord_sf(xlim = c(-73.5, -69.5), ylim = c(4, 6.5), expand = FALSE)

Luego, se carga la librería scales

library(scales)

Attaching package: 㤼㸱scales㤼㸲

The following object is masked from 㤼㸱package:purrr㤼㸲:

    discard

The following object is masked from 㤼㸱package:readr㤼㸲:

    col_factor

Con esta librería se puede agregar color al gráfico mediante el siguiente comando:

ggplot(casanare) + 
  geom_sf(data=casanare_points, aes(x=X, y=Y, fill =
                                       ID_2), color = "black", size = 0.50) +
  geom_text(data = casanare_points, aes(x=X, y=Y,label = ID_2), size = 2.5) +
  theme(aspect.ratio=1)+
  scale_fill_distiller(name="ID_2", palette ="YlOrRd", breaks = pretty_breaks(n = 5))+
  labs(title="Otro Mapa de Casanare")

Ignoring unknown aesthetics: x, y

En el gráfico se pueden visualizar distintos colores para cada municipio dependiendo de su número pero, no puede considerarse como un mapa real. De igual manera, es posible descargar este grafico como un pdf o un png con los siguientes comandos:

ggsave("casanare_municipios.pdf")

Saving 7 x 7 in image

ggsave("map_casanare.png", width = 6, height = 6, dpi = "screen")

**Usar leaflet para visualizar datos:

Para los siguientes pasos, es necesario instalar el paquete leaflet:

#install.packages("leaflet")

Luego, se debe rodar esta función:

library(leaflet)

Ahora, se debe convertir simple features a spatial points:

ant_points <- as(casanare_points, 'Spatial')

Usando la opción head() se puede visualizar lo que contiene el objeto casanare_points.

head(casanare_points)

Simple feature collection with 6 features and 13 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -73.0102 ymin: 4.3456 xmax: -69.92316 ymax: 6.250501
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  ID_0 ISO   NAME_0 ID_1   NAME_1 ID_2       NAME_2    TYPE_2    ENGTYPE_2 NL_NAME_2 VARNAME_2         X
1   53 COL Colombia   10 Casanare  388      Aguazul Municipio Municipality      <NA>      <NA> -72.57165
2   53 COL Colombia   10 Casanare  389      Chámeza Municipio Municipality      <NA>      <NA> -72.88951
3   53 COL Colombia   10 Casanare  390 Hato Corozal Municipio Municipality      <NA>      <NA> -71.20900
4   53 COL Colombia   10 Casanare  391    La Salina Municipio Municipality      <NA>      <NA> -72.36621
5   53 COL Colombia   10 Casanare  392         Maní Municipio Municipality      <NA>      <NA> -72.23162
6   53 COL Colombia   10 Casanare  393    Monterrey Municipio Municipality      <NA>      <NA> -72.87900
         Y                       geometry
1 5.099226 MULTIPOLYGON (((-72.5075 4....
2 5.192376 MULTIPOLYGON (((-72.7872 5....
3 6.077862 MULTIPOLYGON (((-72.0873 6....
4 6.135512 MULTIPOLYGON (((-72.4195 6....
5 4.648025 MULTIPOLYGON (((-72.2622 4....
6 4.820158 MULTIPOLYGON (((-72.93552 5...

Para obtener el área de cada municipio se usa el paquete “lwgeom” que debe ser instalado con el siguiente comando:

#install.packages("lwgeom")

Luego, se debe hacer rodar la biblioteca de la siguiente manera:

library(lwgeom)

Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.6.1, PROJ 4.9.3

Ahora, se escribe el siguiente comando para conocer el área en metros cuadrados(m^2) de cada municipio que conforma Casanare:

mun_casanare$area <- st_area(mun_casanare)

Para cambiar la unidad de medida a kilometros cuadrados (km^2) se escribe:

mun_casanare$km2 <- mun_casanare$area/(1000000)

Se puede verificar este cambio escribiendo el siguiente comando:

mun_casanare$km2

Units: [m^2]
 [1]  1433.3676   328.6694  5063.0322   152.3924  3862.1055   720.9689  1046.8378  4742.7239 12435.9210
[10]   784.1395   179.3627   278.6638   495.2043  3238.2907  1105.4737  2441.8742  3037.4892   755.1099
[19]  2520.3069

Ahora, se debe cambiar de simple featuresa spatial polygons aplicando el siguiente comando:

casanare_mun <- as(mun_casanare, 'Spatial')

Usando la función head() corroboramos la información contenida en el objeto.

#head(casanare_mun)

Ahora, se crean algunos vectores:

bins <- c(0, 50, 100, 200, 300, 500, 1000, 2000, Inf)
pal <- colorBin("PuRd", domain = casanare_mun$km2, bins = bins)


labels <- mun_casanare$NAME_2

labels

 [1] "Aguazul"              "Chámeza"              "Hato Corozal"         "La Salina"           
 [5] "Maní"                 "Monterrey"            "Nunchía"              "Orocué"              
 [9] "Paz de Ariporo"       "Pore"                 "Recetor"              "Sácama"              
[13] "Sabanalarga"          "San Luis de Palenque" "Támara"               "Tauramena"           
[17] "Trinidad"             "Villanueva"           "Yopal"               

Se escribe el siguiente comando para generar un mapa:

m <- leaflet(casanare_mun) %>%
  setView(-71.8, 5.1, zoom = 7.8)  %>% addPolygons(  
  fillColor = ~pal(km2),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#2",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels) %>%
  addLegend(pal = pal, values = ~km2, opacity = 0.7, title = NULL,
    position = "bottomright")

Se escribe la letra m en un chunck para observar el mapa:

m

Otra manera para generar un mapa de manera más sintetizada puede ser:

leaflet() %>%
  addProviderTiles(providers$Esri.WorldImagery, options= providerTileOptions(opacity = 0.99)) %>%
  addPolygons(data = casanare_mun, popup= casanare_mun$NAME_2,
    stroke = TRUE, fillOpacity = 0.25, smoothFactor = 0.25)

Ahora, nos enfocaremos en crear una ilustración para Yopal, la capital del departamento:

#library(leaflet)

En esta gráfica podemos ver un mapa normal del territorio:

l <- leaflet() %>%
  addTiles() %>% addMarkers(lng=-72.403, lat=5.327, popup="Yopal")
l  

Sin embargo, al emplear la opción de leaflet providers el gráfico puede tomar diferentes formas y verse un poc más profesionales. A continuación, se pueden ver dos ejemplos:

library(leaflet)
leaflet() %>% setView(lng = -72.4, lat = 5.35, zoom = 13) %>%
  addProviderTiles(providers$OpenTopoMap)

leaflet() %>%setView(lng = -72.4, lat = 5.35, zoom = 9)%>%
  addProviderTiles(providers$Esri.NatGeoWorldMap, options= providerTileOptions(opacity = 1))