Leyendo, filtrando y vizualizando datos vectoriales geoespaciales en R

Juan Pablo Cañon Sosa

12/09/2020

1. Introduccion

Este es un R Markdown Notebook. Mediante este sketch observaremos algunas funcionalidades geospaciales de R para lo cual usaremos algunas librerias como \(simple\) \(features\) (sf) y \(tidyverse\).

Primero instalamos las librerias nescesarias, para ello usamos el siguiente comando:

#install.package(c("tydiverse", "sf"))

luego de instalar las librerias las cargamos usando el comando:

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

2. Leyendo el vector datos

Una vez descargados los archivos nescesarios desde la pagina oficial del DANE, cargamos el shapefile que contiene los datos de los departamentos de Colombia:

deptos <- read_sf("C:/Users/FLCS/Documents/R/Geomatica Basica/GB/datos/MGN_DPTO_POLITICO.shp")

Visualizamos los datos que contiene el objeto deptos que creamos :

head(deptos)

Usamos la funcion st_crs la cual nos brinda informacion sobre el sistema de referencia de coordenadas que tiene nuestro objeto deptos:

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

3. Usamos ggplot para vizualizar los datos espaciales.

Para entender mejor los datos que estamos trabajando podemos vizualizar graficos usando las funciones de la libreria \(ggplot2\) para plotear.

Instalamos la libreria:

#install.packages("ggplot2")

Cargamos la libreria \(ggplo2\) usando:

library(ggplot2)

Ploteamos los datos de deptos usando:

ggplot() + geom_sf(data = deptos, fill="skyblue3")

Usaremos el sistema de refererencia EPSG 4326 que corresponde al sistema convencional WGS84, empleado para la representacion cartografica a nivel mundial, mas informacion aqui :

ggplot() + geom_sf(data = deptos, fill="skyblue") + coord_sf(crs=st_crs(4326))

Observando las propiedades de el CRS con EPSG codigo 32618. Este corresponde a UTM 18 N. En caso de necesitar usar tal CRS, es nescesario convertir el objeto espacial desde EPSG:4326 en EPSG:32618. usando st_crs:

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

visualizamos los datos del objeto deptos_utm:

deptos_utm

Ploteamos el objeto deptos_utm:

ggplot() + geom_sf(data = deptos_utm, fill="skyblue")

4. Filtrando datos geoespaciales basados en atributos

Como solo estamos interesadosn en el departamento del Valle del Cauca, filtramos los datos:

valle_del_cauca <- deptos %>% filter(DPTO_CNMBR == "VALLE DEL CAUCA")

Ploteamos el nuevo objeto valle_del_cauca:

ggplot() + geom_sf(data = valle_del_cauca, fill="steelblue")

Repetimos los pasos previos para cargar los municipios de Colombia y filtramos solamente los del departamento de Valle del Cauca, los almacenamos en un nuevo objeto que llamaremos muni_valle_del_cauca:

munic <- read_sf("C:/Users/FLCS/Documents/R/Geomatica Basica/GB/datos/76_VALLE_DEL_CAUCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")
mun_valle_del_cauca <- munic %>% filter(DPTO_CNMBR == "VALLE DEL CAUCA")
ggplot() + geom_sf(data = mun_valle_del_cauca, fill="yellow1")

Observamos los datos que contiene el objeto mun_valle_del_cauca:

mun_valle_del_cauca

Puesto que mun_valle_del_cauca no tiene enumerados los municipios, agregamos una columna usando la funcion “cbind”, dicha columna nos permite enumerar las filas la cual llamaremos NMBR_MPIO:

mun_valle_del_cauca=cbind(mun_valle_del_cauca, NMBR_MPIO=rep(seq(1, 42, 1)))

Usamos la funcion st-centroid para configurar el punto centroide del objeto mun_valle_del_cauca entre las extenciones x e y minima y maxima, y almacenando en un nuevo objeto llamado valle_del_cauca_points:

valle_del_cauca_points <- st_centroid(mun_valle_del_cauca)

## Warning in st_centroid.sf(mun_valle_del_cauca): 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

valle_del_cauca_points <- cbind(mun_valle_del_cauca, st_coordinates(st_centroid(mun_valle_del_cauca$geometry)))

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

Ploteamos los el objeto valle_del_cauca con sus municipios enumerados:

ggplot(valle_del_cauca) + 
  geom_sf() +  
  geom_sf(data = valle_del_cauca_points, fill= "orchid") +
  geom_text(data = valle_del_cauca_points, aes(x=X, y=Y, label= NMBR_MPIO), size = 2) +
  coord_sf(xlim = c(-77.6, -75.6), ylim = c(3, 5.1), expand = FALSE) + xlab("Longitud") + ylab("Latitud")

Para los siguientes procesos nescesitaremos una libreria llamada \(scales\).

Instalamos la libreria:

#install.packages("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

Ploteamos otro mapa en el cual se muestran los municipios numerados en diferente escala de colores:

ggplot(valle_del_cauca) +
  geom_sf(data= valle_del_cauca_points, aes(x=X, y=Y, fill= NMBR_MPIO),color="green", size = 0.25) + 
  geom_text(data= valle_del_cauca_points, aes(x=X, y=Y, label= NMBR_MPIO), size = 2) +
  theme(aspect.ratio=1) +
  scale_fill_distiller(name ="N° Mun", breaks= pretty_breaks(n=10)) +
  labs(title="Mapa del Valle del Cauca en escala") + coord_sf(xlim = c(-77.6, -75.6), ylim = c(3, 5.1), expand = FALSE) + xlab("Longitud") + ylab("Latitud")

## Warning: Ignoring unknown aesthetics: x, y

Esta visualizacion no es un mapa real, guardamos este mapa como pdf o png asi:

ggsave("valle_del_cauca_municipios.pdf")

## Saving 7 x 5 in image

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

5. Usando leaflet para vizualizar los datos

Instalamos la libreria \(leaflet\) usando el comando:

#install.packages("leaflet")

Cargamos la libreria usando:

library(leaflet)

Para usar la libreria, nescesitamos convertir el objeto “valle_del_cauca_points” desde “\(simple\) \(features\)” a “\(spatial\) \(points\)”, los almacenamos en un nuevo objeto que llamaremos “vall_points”:

vall_points <- as(valle_del_cauca_points, 'Spatial')

Usamos el siguiente comando para ver que hay dentro del objeto vall_points:

#head(vall_points)

Para ver el area de los municipios instalamos el paquete “\(lwgeom\)”:

#install.packages("lwgeom")

Cargamos el paquete:

library(lwgeom)

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

Vamos a calcular el area de cada municipio (en metros cuadrados) y almacenamos los datos en un nueva columna dentro del objeto mun_valle_del_cauca, que llamaremos “area”:

mun_valle_del_cauca$area <- st_area(mun_valle_del_cauca)

Observamos los primeros datos del objeto “valle_del_cauca”:

head(mun_valle_del_cauca)

Ahora vamos a crear una nueva columna en el objeto mun_valle_del_cauca para almacenar el area en kilometros cuadrados la cual llamaremos km2:

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

Para ver los datos de la columna km2 usamos:

mun_valle_del_cauca$km2

## Units: [m^2]
##  [1]  561.96923  110.31362  305.09656   90.67825  824.82350  396.29099
##  [7]  166.82613  791.91920  295.97445  247.96781  187.58708  213.90532
## [13]  441.27357  235.05969  401.01966  268.36611  162.25518   62.30228
## [19]  214.44782 1003.73955  356.24949  135.98096  307.59329  235.41252
## [25]  201.86058  536.64635  175.54873  307.16046  908.29470   41.82398
## [31]  200.31729  113.57507  327.26523  230.16997  367.50034  264.40976
## [37]  120.72380  254.92278  627.68223  742.38434  914.37039 6274.97442

De nuevo, nescesitamos una conversion desde \(simple\) \(features\) a \(spatial\) \(polygons\) y almacenamos los datos en un nuevo objeto que llamaremos “vall_mun”:

vall_mun <- as(mun_valle_del_cauca, 'Spatial')

Usamos el siguiente codigo para ver los primeros datos que hay dentro del objeto espacial “vall_mun”:

#head(vall_mun)

Ahora, preparamos el plot :

bins <- c(0, 50, 100, 200, 300, 500, 1000, 20000, Inf)
pal <- colorBin(palette="YlGn", domain = vall_mun$area, bins = bins)


labels <- mun_valle_del_cauca$MPIO_CNMBR

labels

##  [1] "CALI"         "ANDALUCÍA"    "ANSERMANUEVO" "ARGELIA"      "BUGA"        
##  [6] "BUGALAGRANDE" "CAICEDONIA"   "CALIMA"       "CANDELARIA"   "CARTAGO"     
## [11] "EL ÁGUILA"    "EL CAIRO"     "EL CERRITO"   "EL DOVIO"     "FLORIDA"     
## [16] "GINEBRA"      "GUACARÍ"      "ALCALÁ"       "OBANDO"       "PALMIRA"     
## [21] "PRADERA"      "RESTREPO"     "RIOFRÍO"      "ROLDANILLO"   "SAN PEDRO"   
## [26] "SEVILLA"      "TORO"         "TRUJILLO"     "TULUÁ"        "ULLOA"       
## [31] "VERSALLES"    "VIJES"        "YOTOCO"       "YUMBO"        "ZARZAL"      
## [36] "LA VICTORIA"  "LA UNIÓN"     "LA CUMBRE"    "JAMUNDÍ"      "BOLÍVAR"     
## [41] "DAGUA"        "BUENAVENTURA"

Entonces, creamos el plot interactivo en el cual podremos saber el nombre de cada municipio deslizando el cursor sobre tal:

m <- leaflet(vall_mun)  %>%
  setView(-76.5, 4, 7.5) %>% addPolygons(
    fillColor = ~pal(km2) ,
    weight = 2,
    opacity = 1,
    color = "skyblue",
    dashArray = "3",
    fillOpacity = 0.7,
    highlight = highlightOptions(
      weight = 5,
      color = "tomato",
      dashArray = "",
      fillOpacity = 0.7,
      bringToFront = TRUE),
    label = labels) %>%
  addLegend(pal= pal, values = ~km2, opacity = 0.7, title = NULL,
            position = "bottomright")

Para ver el plot:

m

Ploteamos el mapa de el departamento del Valle del Cauca usando como base mapas de diferentes proveedores asi:

leaflet() %>%
  addProviderTiles(providers$Esri.WorldImagery, options = providerTileOptions( opacity = 0.99)) %>%
  addPolygons(data = vall_mun, popup = vall_mun$MPIO_CNMBR, fillColor = ~pal(km2) ,
    weight = 2,
    opacity = 1,
    color = "skyblue",
    dashArray = "3",
    fillOpacity = 0.7,
    highlight = highlightOptions(
      weight = 5,
      color = "tomato",
      dashArray = "",
      fillOpacity = 0.7,
      bringToFront = TRUE),
    label = labels, stroke = TRUE, smoothFactor = 0.25)

leaflet() %>%
  addProviderTiles(providers$OpenStreetMap, options = providerTileOptions(opacity = 0.99)) %>%
  addPolygons(data = vall_mun, popup = vall_mun$MPIO_CNMBR,
              stroke = TRUE, weight = 2,  dashArray = "3",
      highlight = highlightOptions(
      weight = 5,
      color = "skyblue",
      dashArray = "",
      fillOpacity = 0.7,
      bringToFront = TRUE),
    label = labels, fillOpacity = 0.1, smoothFactor = 0.25, fillColor = "yellow")

leaflet() %>%
  addProviderTiles(providers$Esri.NatGeoWorldMap, options = providerTileOptions(opacity = 0.99)) %>%
  addPolygons(data = vall_mun, popup = vall_mun$MPIO_CNMBR, fillColor = ~pal(km2) ,
    weight = 2,
    opacity = 1,
    color = "skyblue",
    dashArray = "3",
    fillOpacity = 0.7,
    highlight = highlightOptions(
      weight = 5,
      color = "yellow",
      dashArray = "",
      fillOpacity = 0.7,
      bringToFront = TRUE),
    label = labels, stroke = TRUE, smoothFactor = 0.25)