Lectura, filtrado y visualización de datos vectoriales geoespaciales (Meta)

Introducción

El trabajo que se llevará a cabo se realizó con base en la información adquirida en la clase de geomática básica, en la Universidad Nacional de Colombia en 2020-1. Es necesario instalar los paquetes tidyverse y sf, además cargar sus librerías.

install.packages("tidyverse", "sf", "CRAN")

install.packages("tidyverse", "sf", "CRAN")

Error in install.packages : Updating loaded packages

2. Lectura de datos vectoriales

Ahora se puede leer un archivo que representa a los departamentos de Colombia, descargado de DIVA-GIS (consultado en el buscador), click en Datos espaciales gratuitos y luego datos a nivel de país. Se selecciona el país con el que se va a trabajar, en este caso Colombia, y luego los temas en los que se tiene interés. Se descargó el tema: Áreas administrativas, se debe descomprimir la carpeta e indicar su ubicación para que pueda ser leída.

library("tidyverse", "sf", "read", "CRAN")

Para saber cuál es el contenido del objeto deptos, denominado anteriormente.

deptos <-  read_sf("D:/Documentos/Geomática/Meta/COL_adm/COL_adm1.shp")

La biblioteca sf ofrece diversas funciones tales como saber cuál es el sistema de referencia de coordenadas de los datos vectoriales almacenados en el objeto deptos, esto por medio de la función st_crs:

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
geographic CRS: WGS 84

En el cuadro se puede apreciar que el sistema que se maneja es WGS84. ### 3. Usar ggplot para visualizar los datos geoespaciales Para trazar los datos se usa la función ggplot. Se podrá apreciar el mapa de Colombia y las divisiones de los departamentos que lo conforman, además de algunas coordenadas para su ubicación.

st_crs(deptos)

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]]

Es posible utilizar cualquier sistema de referencia para trazar los datos. Sin embargo, no es correcto usar un sistema de referencia de coordenadas que se haya definido explícitamentepara otro país o región. Se puede encontrar información dobre los códigos EPSG en este link: https://epsg.io/

ggplot() + geom_sf(data = deptos) 

Busque las propiedades del CRS con el código EPSG 32618. Corresponde a UTM 18 N. En caso de que necesitemos usar dicho CRS, es necesario convertir el objeto espacial de EPSG4326 a EPSG: 32618.

#  El CRS 3978 se usa en Canadá

Warning message:
In readChar(file, size, TRUE) : truncating string with embedded nuls

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

deptos_utm

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
projected CRS:  WGS 84 / UTM zone 18N

4. Filtrar datos geoespaciales basados en atributos

En esta ocasión se visualizará el departamento del Meta, para ello se debemn filtrar los datos.

ggplot() + geom_sf(data = deptos_utm)

Se trazará el nuevo objeto:

meta <-  deptos %>%   filter(NAME_1 == "Meta")

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls

Se puede repetir el procedimiento anterior para cargar los datos de los municipios de los departamentos de Colombia y así poder filtrar los del Meta. }

ggplot() + geom_sf(data = meta) 

munic <-  read_sf("D:/Documentos/Geomática/Meta/COL_adm/COL_adm2.shp")
mun_meta <- munic %>% filter(NAME_1 == "Meta")
ggplot() + geom_sf(data = mun_meta) 

mun_meta

Simple feature collection with 29 features and 11 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -74.9206 ymin: 1.570401 xmax: -71.0749 ymax: 4.860201
geographic CRS: WGS 84

meta_points <- cbind(mun_meta, st_coordinates(st_centroid(mun_meta$geometry)))

Se puede producir una mejor representación utilizando el siguiente fragmento (vale aclarar que las coordenadas deben ser ajustadas para cada departamento con el que se desee trabajar, se puede guiar de los mapas anteriores):

meta_points <- cbind(mun_meta, st_coordinates(st_centroid(mun_meta$geometry)))

st_centroid does not give correct centroids for longitude/latitude data

ggplot(meta) +
    geom_sf() +
    geom_sf(data = meta_points, fill = "antiquewhite") + 
    geom_text(data = meta_points, aes(x=X, y=Y,label = ID_2), size = 2) +
    coord_sf(xlim = c(-75, -70.6), ylim = c(1.5, 5), expand = FALSE)

Warning in readChar(file, size, TRUE) :
  truncating string with embedded nuls
Warning in readChar(file, size, TRUE) :
  truncating string with embedded nuls
Warning in readChar(file, size, TRUE) :
  truncating string with embedded nuls
Warning in readChar(file, size, TRUE) :
  truncating string with embedded nuls
Warning in readChar(file, size, TRUE) :
  truncating string with embedded nuls
Warning in readChar(file, size, TRUE) :
  truncating string with embedded nuls

library(scales)

Attaching package: 㤼㸱scales㤼㸲

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

    discard

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

    col_factor

There were 12 warnings (use warnings() to see them)

Sin embargo, esta visualización no es un mapa real. De todos modos, la salida se puede guardar como un pdf o como un png. Lo anterior por medio de la función ggsave.

ggplot(meta) + 
  geom_sf(data=meta_points, aes(x=X, y=Y, fill =
                                       ID_2), color = "black", size = 0.25) +
  geom_text(data = meta_points, aes(x=X, y=Y,label = ID_2), size = 2) +
  theme(aspect.ratio=1)+
  scale_fill_distiller(name="ID_2", palette = "YlGn", breaks = pretty_breaks(n = 5))+
  labs(title="Otro mapa de Meta")

Ignoring unknown aesthetics: x, y

ggsave("antioquia_municipios.pdf")

Saving 7 x 7 in image

5. Usando el folleto para visualizar datos

Se debe instalar el paquete “leaflet”

#install.packages("leaflet")

library(leaflet)

Para usar la biblioteca, necesitamos convertir de características simples a puntos espaciales.

library(leaflet)

package 㤼㸱leaflet㤼㸲 was built under R version 4.0.2

Para saber qué hay dentro del objeto meta_points:

meta_points <- as(meta_points, 'Spatial')

Obtener área de municipios:

head(meta_points)

Cargar librería:

install.packages("lwgeom")

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into 㤼㸱D:/Documentos/R/win-library/4.0㤼㸲
(as 㤼㸱lib㤼㸲 is unspecified)
probando la URL 'https://cran.rstudio.com/bin/windows/contrib/4.0/lwgeom_0.2-5.zip'
Content type 'application/zip' length 5100534 bytes (4.9 MB)
downloaded 4.9 MB

package ‘lwgeom’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\prueba\AppData\Local\Temp\Rtmpq6FmKP\downloaded_packages

Calculemos el área de cada municipio (en metros cuadrados):

library(lwgeom)

package 㤼㸱lwgeom㤼㸲 was built under R version 4.0.2Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.8.0, PROJ 6.3.1

Ahora, creemos un nuevo campo para almacenar el área en kilómetros cuadrados:

mun_meta$area <- st_area(mun_meta) #Cuidar unidades

Verifique la salida:

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

Nuevamente, necesitamos una conversión de características simples a polígonos espaciales:

mun_meta$km2

Units: [m^2]
 [1]   975.4030   424.6030   808.6122
 [4]   561.5467   379.1041   197.7810
 [7]   719.5302   105.7460   754.1116
[10]   229.6312   764.3242 11025.7283
[13]  6298.7609   969.2610 12455.3204
[16]  2119.5518   918.9151 16788.4326
[19]  5576.2151  2084.8685  4061.3990
[22]   566.5750  1437.4442   834.0937
[25]   202.6045  1686.2457  6511.0136
[28]  1369.7743  4498.9691

meta_mun <- as(mun_meta, 'Spatial')

Preparar la representación:

head(meta_mun)

Luego, crea la representación (tener en cuenta las coordenadas):

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


labels <- mun_meta$NAME_2

labels

 [1] "Acacías"             
 [2] "Barranca de Upía"    
 [3] "Cabuyaro"            
 [4] "Castilla la Nueva"   
 [5] "Cumaral"             
 [6] "El Calvario"         
 [7] "El Castillo"         
 [8] "El Dorado"           
 [9] "Fuente de Oro"       
[10] "Granada"             
[11] "Guamal"              
[12] "La Macarena"         
[13] "La Uribe"            
[14] "Lejanías"            
[15] "Mapiripán"           
[16] "Mesetas"             
[17] "Puerto Concordia"    
[18] "Puerto Gaitán"       
[19] "Puerto López"        
[20] "Puerto Lleras"       
[21] "Puerto Rico"         
[22] "Restrepo"            
[23] "San Carlos de Guaroa"
[24] "San Juan de Arama"   
[25] "San Juanito"         
[26] "San Luis de Cubarral"
[27] "San Martín"          
[28] "Villavicencio"       
[29] "Vista Hermosa"       

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

Otra forma de trazar. Puede ser más simple:

m

Puede probar diferentes proveedores para mejorar su mapa. ¡Aproveche la función de completar pestañas para seleccionar el mapa base preferido con solo desplazarse por la lista de 110 proveedores!

Descomente el siguiente fragmento, complete el código y cree una visualización hermosa para la capital de su departamento. En caso de problemas, use la ayuda de R.

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

leaflet() %>%
  addTiles %>% #Add deafault OpenStreetMap map ties 
  addMarkers(lng=-73.628177, lat=4.148539, popup="Villavicencio")

leaflet() %>%

There were 26 warnings (use warnings() to see them)

  addTiles %>% #Add deafault OpenStreetMap map ties 
  addMarkers(lng=-73.628177, lat=4.148539, popup="Villavicencio")

foo <- leaflet() %>% 
    setView(lng =-73.628177, lat =4.148539, zoom = 10) %>%
    addTiles() 

foo %>% addProviderTiles("NASAGIBS.ModisTerraTRUEColorCR")

Referencias:

https://www.diva-gis.org/

https://ials.github.io/geomatica/

https://rpubs.com/geo2/vectordata