Leer, filtrar y visualizar datos vectoriales geoespaciales en R
Laura Juliana Peñuela Mojica
14 de septiembre de 2020
#1.Introducción
¡Hola!, mi nombre es Laura Peñuela. En este cuaderno de R mostraré las funciones básicas que se encuentran en las librerías “simple feature” (sf) y tidyverse. El objetivo es obtener diferentes mapas de Colombia, Cundinamraca y Bogotá por medio de la manipulación de archivos shape.
El primer pasa es instalar y cargar las librerías requeridas
#install.packages("tidyverse")
#install.packages("sf")library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.2 v purrr 0.3.4
## v tibble 3.0.2 v dplyr 1.0.0
## v tidyr 1.1.0 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.Leer los datos vectoriales Al tener las librerías cargadas ya se pueden leer archivos tipo Shape. Se deben descargar los datos de Colombia de la división administrativa en el sigueinte link: https://www.diva-gis.org/gdata.Posterior a ello, se cargan con la función “read_sf”. Se debe poner un nombre al objeto, en este caso es llamado “deptos”
deptos <- read_sf("C:/Users/JUAN PABLO/Documents/geo maps/COL_adm1.shp")Aplicamos la función “head”, para mirar las primeras filas y de este modo conocer las categorías que presenta el archivo Shape
head(deptos)Como los datos corresponden a datos espaciales es útil ver el sistema de referencia en el que ellos se encuentran. Para ello se debe aplicar la función “st_crs” al objeto
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.Visualizar los datos espaciales por medio de ggplot
La función ggplot sirve para vincular un gráfico a un conjunto de datos. El mapa siguiente corresponde a la representación de los datos del objeto “deptos”. Como se puede observar el mapa de Colombia se ubica en unas coordenadas de acuerdo al sistema de referencia.
ggplot() + geom_sf(data = deptos, fill="darkseagreen")
Los datos pueden ser representados en cualquier sistema de referencia. En este caso se usó el sistema de referencia CRS3978 de Canadá con el objetivo de ejemplificar lo mencionado anteriormente. Sin embargo, no se recomienda usar un sistema de referencia que se haya definido explícitamene para otro país o región
ggplot() + geom_sf(data = deptos, fill="palegreen") + coord_sf(crs=st_crs(3978))
Es necesario convertir el objeto “deptos” al sistema coordenado de referencias de EPSG:4326 a EPSG:32618 para que coincida con el sistema de coordenadas universal transversal de Mercator (UTM)
deptos_utm <- st_transform(deptos, crs = st_crs(32618))
deptos_utmAhora se visualiza el nuevo mapa con el sistema de coordenadas UTM por medio de la función ggplot
ggplot() + geom_sf(data = deptos_utm, fill="palegreen3")
#4. Filtrar datos espaciales según los atributos
En este caso el interés es conocer los datos exclusivos para el departamento de Cundinamarca. Esto se realiza al usar el conjunto de símbolos %>%, como se observa en el siguiente cuadro de código. Posterior a filtrar se visualiza el mapa con ggplot
cundinamarca <- deptos %>% filter(NAME_1 == "Cundinamarca")
ggplot() + geom_sf(data = cundinamarca, fill="plum1") 
Se pueden usar los pasos anteriores para filtrar la división administrativa de los municipios de Cundinamrca. Para ello necesitamos leer los datos sf llamados COL_adm2.shp que contienen esta información. Nuevamente visualizamos con ggplot y así podemos ver los municipios.
municipios <- read_sf("C:/Users/JUAN PABLO/Documents/geo maps/COL_adm2.shp")
municipiosmun_cundi <- municipios %>% filter(NAME_1 == "Cundinamarca")
ggplot() + geom_sf(data = mun_cundi, fill="orchid1") 
Ahora se lee objeto para ver las catogorías que contienen los datos sf
mun_cundiSe crea el objeto “cundinamrca_ponits” y se aplica la función “st_centrodide” que provee el centroide en una geometría
cundinamarca_points<- st_centroid(mun_cundi)## Warning in st_centroid.sf(mun_cundi): 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 datacundinamarca_points <- cbind(mun_cundi, st_coordinates(st_centroid(mun_cundi$geometry)))## Warning in st_centroid.sfc(mun_cundi$geometry): st_centroid does not give
## correct centroids for longitude/latitude dataCon el siguiente cuadro de código se puede obtener un mapa más detallado de Cundinamarca con sus respectivos municipios.
ggplot(cundinamarca) +
geom_sf() +
geom_sf(data = cundinamarca_points, fill = "lightpink") +
geom_text(data = cundinamarca_points, aes(x=X, y=Y,label = ID_2), size = 1) +
coord_sf(xlim = c(-75.2, -72.7), ylim = c(3.5, 6), expand = FALSE)
Instalar y cargar la librería “scales”. Este paquete permite implementar escalas de una manera independete al sistema de gráficos.
#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_factorPor medio de este cuadro de código se puede realizar un mapa que contenga la información de la columna ID_2 y a su vez asingar una escala de color de acuerdo con el número asignado a cada municipio.
ggplot(cundinamarca) +
geom_sf(data=cundinamarca_points, aes(x=X, y=Y, fill =
ID_2), color = "black", size = 0.25) +
geom_text(data = cundinamarca_points, aes(x=X, y=Y,label = ID_2), size = 1) +
theme(aspect.ratio=1)+
scale_fill_distiller(name="ID_2", palette = "BuPu", breaks = pretty_breaks(n = 5))+
labs(title="Mapa de Cundinamarca")## Warning: Ignoring unknown aesthetics: x, y
Esta visualización puede ser guardada como pdf o png por medio de la función “ggsave”
ggsave("cundinamarca_municipios.pdf")## Saving 7 x 5 in imageggsave("map_cundinamarca.png", width = 6, height = 6, dpi = "screen")#5.Usar la librería “leaflet” para visualizar datos.
Primero se debe instalar y cargar la librería
#install.packages("leaflet")library(leaflet)Para usar la biblioteca es necesario convertir el objeto “cundi_points” a puntos espaciales.
#cun_points <- as(cundi_points, 'Spatial')Para obtener el área de los municipios es necesario cargar la librería “lwgeom”
#install.packages("lwgeom")library(lwgeom)## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.8.0, PROJ 6.3.1El siguiente cuadro de código calcula el área de cada municipio en metros cuadrados
mun_cundi$area <- st_area(mun_cundi)Ahora se obtiene el área de cada municipio pero en kilométros cuadrados
mun_cundi$km2 <- mun_cundi$area/(1000000)Observar la salida del objeto creado
mun_cundi$km2## Units: [m^2]
## [1] 76.15363 66.44340 119.21639 119.85758 107.47458 184.38331
## [7] 211.57878 52.30997 75.99254 117.53693 457.78629 59.40923
## [13] 49.09098 592.53649 306.07664 81.08696 173.19669 119.07608
## [19] 215.35620 319.86611 127.49569 53.16071 91.75662 120.19914
## [25] 119.12324 476.79685 155.70691 74.20352 116.43546 69.93551
## [31] 307.46709 422.85426 43.60324 353.11446 102.47574 106.21432
## [37] 179.98855 785.82995 324.85095 119.43901 139.57971 62.74291
## [43] 211.58357 607.42953 222.64887 335.18189 348.70941 157.17983
## [49] 209.69651 135.08022 161.57786 145.41077 220.72771 119.14988
## [55] 101.99121 1633.38987 94.86542 117.17884 99.77867 205.90696
## [61] 74.27546 60.67551 393.49430 171.02753 100.00624 776.64209
## [67] 146.68829 567.06836 186.86541 77.23026 202.57553 130.78296
## [73] 147.32023 89.47211 193.42735 264.71363 123.48160 316.43054
## [79] 1650.84698 100.16249 143.14552 105.48010 174.79363 106.28521
## [85] 252.18718 119.75797 288.11260 187.33163 129.79540 115.06902
## [91] 59.40725 73.51989 214.98681 82.08100 100.60445 92.41850
## [97] 64.98824 262.08657 66.39905 154.53974 318.62965 112.95061
## [103] 114.96535 89.05010 135.46022 152.43358 77.24548 111.08833
## [109] 75.71577 198.60633 147.27553 198.06351 921.84588 55.61708
## [115] 188.60388Se debe convertir el objeto cun_mun a poliginos espaciales
cun_mun <- as(mun_cundi, 'Spatial')Seguido a ello, se debe preparar el plot
bins <- c(0, 50, 100, 200, 300, 500, 1000, 2000, Inf)
pal <- colorBin("YlOrRd", domain = cun_mun$km2, bins = bins)
labels <- mun_cundi$NAME_2
labels## [1] "Agua de Dios" "Albán"
## [3] "Anapoima" "Anolaima"
## [5] "Apulo" "Arbeláez"
## [7] "Beltrán" "Bituima"
## [9] "Bojacá" "Cáqueza"
## [11] "Cabrera" "Cachipay"
## [13] "Cajicá" "Caparrapí"
## [15] "Carmen de Carupa" "Chía"
## [17] "Chaguaní" "Chipaque"
## [19] "Choachí" "Chocontá"
## [21] "Cogua" "Cota"
## [23] "Cucunubá" "El Colegio"
## [25] "El Peñon" "Fómeque"
## [27] "Facatativá" "Fúquene"
## [29] "Fosca" "Funza"
## [31] "Fusagasugá" "Gachalá"
## [33] "Gachancipá" "Gachetá"
## [35] "Gama" "Girardot"
## [37] "Guachetá" "Guaduas"
## [39] "Guasca" "Guataquí"
## [41] "Guatavita" "Guayabal de Síquima"
## [43] "Guayabetal" "Gutiérrez"
## [45] "Jerusalén" "Junín"
## [47] "La Calera" "La Mesa"
## [49] "La Palma" "La Peña"
## [51] "La Vega" "Lenguazaque"
## [53] "Machetá" "Madrid"
## [55] "Manta" "Medina"
## [57] "Mosquera" "Nariño"
## [59] "Nemocón" "Nilo"
## [61] "Nimaima" "Nocaima"
## [63] "Pacho" "Paime"
## [65] "Pandi" "Paratebueno"
## [67] "Pasca" "Puerto Salgar"
## [69] "Pulí" "Quebradanegra"
## [71] "Quetame" "Quipile"
## [73] "Ricaurte" "San Antonio del Tequendama"
## [75] "San Bernardo" "San Cayetano"
## [77] "San Francisco" "San Juan de Río Seco"
## [79] "Santafé de Bogotá" "Sasaima"
## [81] "Sesquilé" "Sibaté"
## [83] "Silvania" "Simijaca"
## [85] "Soacha" "Sopó"
## [87] "Subachoque" "Suesca"
## [89] "Supatá" "Susa"
## [91] "Sutatausa" "Tabio"
## [93] "Tausa" "Tena"
## [95] "Tenjo" "Tibacuy"
## [97] "Tibirita" "Tocaima"
## [99] "Tocancipá" "Topaipí"
## [101] "Ubalá" "Ubaque"
## [103] "Une" "Utica"
## [105] "Venecia" "Vergara"
## [107] "Vianí" "Villa de San Diego de Ubaté"
## [109] "Villagómez" "Villapinzón"
## [111] "Villeta" "Viotá"
## [113] "Yacopí" "Zipacón"
## [115] "Zipaquirá"Ahora se crea el plot
m <- leaflet(cun_mun) %>%
setView(-75.5, 7, 8) %>% addPolygons(
fillColor = ~pal(km2),
weight = 2,
opacity = 1,
color = "black",
dashArray = "3",
fillOpacity = 0.7,
highlight = highlightOptions(
weight = 5,
color = "PuRd",
dashArray = "",
fillOpacity = 0.7,
bringToFront = TRUE),
label = labels) %>%
addLegend(pal = pal, values = ~km2, opacity = 0.7, title = "Cundimarca_Área",
position = "bottomright")Se llama al objeto que da lugar al siguiente mapa donde se puede observar la división administrativa de Cundinamrca y a su vez el tamaño aproximado de los municipios en kilómetros cuadrados gracias a la escala de color
mleaflet() %>%
addProviderTiles(providers$Esri.WorldImagery, options= providerTileOptions(opacity = 0.99)) %>%
addPolygons(data = cun_mun, popup= cun_mun$NAME_2,
stroke = TRUE, fillOpacity = 0.25, smoothFactor = 0.25
)leaflet() %>%
addProviderTiles(providers$Esri.NatGeoWorldMap, options= providerTileOptions(opacity = 1)) %>%
addPolygons(data = cun_mun, popup= cun_mun$NAME_2,
stroke = TRUE, fillOpacity = 0.25, smoothFactor = 1, weight = 2
)Otra visualización similar
leaflet() %>%
addProviderTiles(providers$Esri.WorldStreetMap, options= providerTileOptions(opacity = 1)) %>%
addPolygons(data = cun_mun, popup= cun_mun$NAME_2,
stroke = TRUE, fillOpacity = 0.25, smoothFactor = 1, weight = 2
)A paritr de este código se puede observar una vista de la ciduad capital de Cundimarca (Bogotá). Se realiza un procedimiento similiar al realizado cuando se visualizó Cundinamarca
municipios <- read_sf("C:/Users/JUAN PABLO/Documents/geo maps/COL_adm2.shp")
head(municipios)bogota <- municipios %>% filter(NAME_2 == "Santafé de Bogotá")
ggplot() + geom_sf(data = bogota, fill="orchid")
leaflet() %>%
addProviderTiles(providers$Esri.NatGeoWorldMap, options= providerTileOptions(opacity = 0.99)) %>%
addPolygons(data = bogota, popup= bogota$NAME_2,
stroke = TRUE, fillOpacity = 0.25, smoothFactor = 0.25, weight = 2
)rpubs.upload.method="internal"