Lectura, filtrado y visualización de datos vectoriales geoespaciales en R para el departamento de Huila

Introducción

Este es un R Markdown Notebook que elabora un mapa apoyado en las funciones básicas proporcionadas por las funciones simples (sf) y las bibliotecas tidyverse.

Para empezar es necesario instalar las librerías apropiadas:

#install.packages(c("tidyverse", "sf"))

Ahora cargamos las librerias instaladas previamente:

library(tidyverse)

## -- Attaching packages ---------------------------------------- tidyverse 1.3.0 --

## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.2
## v tidyr   1.1.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()

library(sf)

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

Leyendo los datos vectoriales

Ahora, vamos a leer un shapefile de datos del departamento descargados previamente.

munpios <- read_sf("C:/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Para ver el contenido del objeto definido atrás utilizamos:

head(munpios)

## Simple feature collection with 6 features and 9 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -76.4261 ymin: 1.853989 xmax: -74.94464 ymax: 3.284017
## geographic CRS: WGS 84
## # A tibble: 6 x 10
##   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR
##   <chr>      <chr>      <chr>      <chr>           <dbl>     <int> <chr>     
## 1 41         41001      NEIVA      1612            1270.      2017 HUILA     
## 2 41         41306      GIGANTE    1789             504.      2017 HUILA     
## 3 41         41319      GUADALUPE  Ordenanza~       250.      2017 HUILA     
## 4 41         41349      HOBO       1884             194.      2017 HUILA     
## 5 41         41357      ÍQUIRA     1887             357.      2017 HUILA     
## 6 41         41359      ISNOS      Ordenanza~       371.      2017 HUILA     
## # ... with 3 more variables: Shape_Leng <dbl>, Shape_Area <dbl>,
## #   geometry <POLYGON [°]>

Ahora, es útil conocer las características de los datos para conocer cuál es el sistema de referencia de coordenadas de los datos vectoriales almacenados en el objeto. Esto se consigue a través de la función st_crs de la biblioteca sf.

st_crs(munpios)

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

Usando ggplot para visualizar datos geoespaciales:

Para el primer plot que veremos, podemos usar las funcionalidades de ggplot:

ggplot() + geom_sf(data = munpios)

Para utilizar un sistema de referencia de coordenadas es necesario conocer el apropiado para la región de interés.

# Advierta la infuencia en el plot cuando se utiliza el CRS manejado en Canada
ggplot() + geom_sf(data = munpios) + coord_sf(crs=st_crs(3978))

Las propiedades del CRS con código EPSG 32618. Corresponde a UTM 18 N. Si se necesitára dicho CRS, es necesario convertir el objeto espacial de EPSG4326 en EPSG: 32618.

munpios_utm <- st_transform(munpios, crs = st_crs(32618))

munpios_utm

## Simple feature collection with 37 features and 9 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 319297.7 ymin: 171589.9 xmax: 565179 ymax: 424810.4
## projected CRS:  WGS 84 / UTM zone 18N
## # A tibble: 37 x 10
##    DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR
##  * <chr>      <chr>      <chr>      <chr>           <dbl>     <int> <chr>     
##  1 41         41001      NEIVA      1612            1270.      2017 HUILA     
##  2 41         41306      GIGANTE    1789             504.      2017 HUILA     
##  3 41         41319      GUADALUPE  Ordenanza~       250.      2017 HUILA     
##  4 41         41349      HOBO       1884             194.      2017 HUILA     
##  5 41         41357      ÍQUIRA     1887             357.      2017 HUILA     
##  6 41         41359      ISNOS      Ordenanza~       371.      2017 HUILA     
##  7 41         41378      LA ARGENT~ Ordenanza~       356.      2017 HUILA     
##  8 41         41396      LA PLATA   Ordenanza~       815.      2017 HUILA     
##  9 41         41483      NÁTAGA     Ordenanza~       132.      2017 HUILA     
## 10 41         41503      OPORAPA    Ordenanza~       156.      2017 HUILA     
## # ... with 27 more rows, and 3 more variables: Shape_Leng <dbl>,
## #   Shape_Area <dbl>, geometry <POLYGON [m]>

ggplot() + geom_sf(data = munpios_utm)

Filtrado de datos geoespaciales basados en atributos

La funcion st-centroid permite modificar el punto centroide que tiene definida el objeto y lo remplaza por otro valor.

huila_points<- st_centroid(munpios)

## Warning in st_centroid.sf(munpios): st_centroid assumes attributes are constant
## over geometries of x

## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for longitude/
## latitude data

huila_points <- cbind(munpios, st_coordinates(st_centroid(munpios$geometry)))

## Warning in st_centroid.sfc(munpios$geometry): st_centroid does not give correct
## centroids for longitude/latitude data

Para asignar su respectivo código a cada municipio y cambiar el estilo del mapa:

ggplot(munpios) +
  geom_sf() +
  geom_sf(data = huila_points, fill = "antiquewhite") + 
  geom_text(data = huila_points, aes(x=X, y=Y,label = MPIO_CCDGO), size = 2) +
  coord_sf(xlim = c(-76.7, -74.3), ylim = c(1.49, 4), expand = FALSE)

library(scales)

## 
## Attaching package: 'scales'

## The following object is masked from 'package:purrr':
## 
##     discard

## The following object is masked from 'package:readr':
## 
##     col_factor

De nuevo, ploteamos un mapa con mejores características:

ggplot(munpios) + 
  geom_sf(data=huila_points, aes(x=X, y=Y, fill =
                                       MPIO_CCDGO), color = "black", size = 0.25) +
  geom_text(data = huila_points, aes(x=X, y=Y,label = MPIO_CCDGO), size = 2) +
  theme(aspect.ratio=1)+
  labs(title="Mapa de Huila")

## Warning: Ignoring unknown aesthetics: x, y

Para guardar los mapas visualizados lo puede hacer en pdf o como png:

ggsave("munpios.pdf")

## Saving 7 x 5 in image

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

Usando leaflet para visualizar datos

Intalamos la librería necesaria:

# install.packages("leaflet")

Cargamos la librería:

library(leaflet)

Esta biblioteca requiere convertir de características simples a puntos espaciales.

hui_points <- as(huila_points, 'Spatial')

Para ver qué hay dentro del objeto tecleamos:

#head(hui_points)

Para conocer el area de los munucipios instalamos el paquete:

# install.packages("lwgeom")

Cargamos la librería:

library(lwgeom)

## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.8.0, PROJ 6.3.1

El cálculo del área de cada municipio (en metros cuadrados) se obtiene asi:

munpios$area <- st_area(munpios)

Y para el área en kilómetros cuadrados:

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

Para ver los datos colocados en kilometros cuadrados:

munpios$km2

## Units: [m^2]
##  [1] 1269.26174  503.37245  249.40676  194.12171  356.23768  370.11715
##  [7]  355.96522  813.72650  131.82205  155.62859  278.07835  883.76179
## [13]  177.98341  193.47749  629.85959  251.12693  465.98246 1388.09923
## [19]  339.12792  429.27276  362.45241  366.73566  530.23318  471.34100
## [25]  185.53442  543.73941  332.70453  606.11657   80.35526 1584.78955
## [31]  462.57691  785.81395  180.81700  589.42196  795.23340  273.50649
## [37]  539.90610

Ahora volvemos a caracteristicas simples con el comando:

hui_mun <- as(munpios, 'Spatial')

Para revisar los datos del comando anterior:

#head(hui_mun)

Para elaborar el plot:

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


labels <- munpios$MPIO_CNMBR

labels

##  [1] "NEIVA"        "GIGANTE"      "GUADALUPE"    "HOBO"         "ÍQUIRA"      
##  [6] "ISNOS"        "LA ARGENTINA" "LA PLATA"     "NÁTAGA"       "OPORAPA"     
## [11] "PAICOL"       "PALERMO"      "PALESTINA"    "PITAL"        "PITALITO"    
## [16] "RIVERA"       "SALADOBLANCO" "SAN AGUSTÍN"  "SANTA MARÍA"  "SUAZA"       
## [21] "TARQUI"       "TESALIA"      "TELLO"        "TERUEL"       "TIMANÁ"      
## [26] "VILLAVIEJA"   "YAGUARÁ"      "GARZÓN"       "ELÍAS"        "COLOMBIA"    
## [31] "CAMPOALEGRE"  "BARAYA"       "ALTAMIRA"     "ALGECIRAS"    "AIPE"        
## [36] "AGRADO"       "ACEVEDO"

Elaboración del plot. Posibilidad para asignar características de diseño:

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

Para ver el mapa:

Está no es la única manera de plotear. Para proporcionar mejoras de ploteado:

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

Lectura, filtrado y visualización de datos vectoriales geoespaciales en R para el departamento de Huila

Sergio García

14/9/2020