Introducción.

El objetivo del presente trabajo es presentar el proceso realizado para la representación espacial del departamento del Meta en Colombia, a traves de las funcionalidades geoespaciales de R. Teniedo en cuenta lo anterior, las herraminetas necesarias para el desarrollo de este trabajo que se encuentran compiladas en la nube en forma de librería son: Simple Features (sf) y Tidyverse, entre otras, por lo que es necesario desacrgar las mismas, así:

1. Instalando Tidyverse y sf:

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

2. De igual forma es necesario cargar en R las librerías recien descargadas:

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 el vector con los datos.

Una vez habiendo descargado el Shape File con los datos de los departamentos de Colombia, es necesario leer el mismo en el cuaderno, así:

departamentos<- read_sf("C:/Users/EQUIPO/Documents/2020-2/Geomática/Trabajo_1/MGN_DPTO_POLITICO.shp")

Según esto, leyendo los datos almacenados en la variable “departamentos”:

head(departamentos)

Es importante identifcar cual es el sistema de coordenadas de referencia del vector “departamentos”, para ello, se usa de las funciones de sf library conocida como st_crs, mostarda a continuación:

st_crs(departamentos)
## 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]]

El resultado del anterior comando permite observar que el sistema geodésico de coordenadas para “departamentos” es EPSG 4326 que corresponde al convencional sistema WGS 84, esto es, un patrón matemático que representa a la tierra por medio de un elipsiode de parámetros definidos y que, además, es uno de los modelos de representación de la tierra más acogidos mundialmente, por lo que, según el objetivo de este informe, resulta ser muy adecuado para este caso. Sin embargo, así como se verá en el siguiente numeral de este cuaderno, el sistema de referencia se puede cambiar acorde con las necesidades del usuario.

Usando ggplot para visualizar los datos geoespaciales.

Así como se mencionó en la introducción de este cuaderno-R, para poder representar gráficamente los datos espaciales, se puede usa un paquete conocido como **ggplot2*, por lo que, a continuación, se procede a llamar a la interfaz de R al mismo:

library(ggplot2)

De modo que, mediante las funcionalidades de ggplot2, se procede a graficar los datos almacenados en el vector departamentos.

#ggplot() + geom_sf(data = departamentos, fill="azure2")

Mapa de Colombia

De igual forma, se puede uasr cualquier sistema de coordenadas para graficar los datos, sin embargo, no es correcto usar otro sitema de cooredenads que ha sido expilictamente definido para otra región, se puede encontrar más información aquí. A continuación un ejemplo:

# El CRS 5940 es usado en Rusia
#ggplot() + geom_sf(data = departamentos, fill="bisque2") + coord_sf(crs=st_crs(5940))

Mapa de Colombia

Otro posibilidad es convertir todos los datos para, no solo la gráfica, a otro CRS de código EPSG, por ejemplo: 3112, especifico para Australia, de la siguiente forma:

deptos_new <- st_transform(departamentos, crs = st_crs(3112))

Llamando a “deptos_new”:

deptos_new

Graficando a “deptos_new”:

#ggplot() + geom_sf(data = deptos_new, fill="azure")

Mapa de Colombia

Con los dos ejemplos anteriores se logró observar una de la características más importantes que tiene ggplot al momento de tratar datos geoespaciales.

Filtrando datos geoespaciales basados en atributos.

Puesto que el interés de este trabajo se basa específicamente en el departamento del Meta, se pueden filtar los datos a partir del siguiente comando:

meta<- departamentos %>% filter(DPTO_CNMBR=="META")

Realizando el gráfico de “meta”:

#ggplot() + geom_sf(data = meta, fill="lemonchiffon2")

Mapa del Meta

Ahora bien, es propicio filtrar los municipios de Colombia que se en cuentran en el departamento del Meta, para ello, se procede a crear una variable que almacene a todos los muncicipos del Meta:

mun_meta<- read_sf("C:/Users/EQUIPO/Documents/2020-2/Geomática/Trabajo_1/50_META/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Llamando a “mun_meta”:

mun_meta

Graficando el mapa del departamento del Meta con sus respectivas divisiones políticas por municipios:

#ggplot() + geom_sf(data = mun_meta, fill="darkseagreen1")

Mapa del Meta

Mediante el comando cbind, se añade una nueva comlumna que permita numerar la cantidad de municipios del departamento del Meta.

mun_meta<-cbind(mun_meta, n_cpi=rep(seq(1, 29, 1)))

Llamando nuevamente a mun_meta, para revisar el número de municipios:

mun_meta

Una vez teniendo numerados los municipios, se procede a configurar el punto centroide de mun_meta, limtado por su extensión tanto en X como en Y:

meta_points <- st_centroid(mun_meta)
## Warning in st_centroid.sf(mun_meta): 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
meta_points <- cbind(mun_meta, st_coordinates(st_centroid(mun_meta$geometry)))
## Warning in st_centroid.sfc(mun_meta$geometry): st_centroid does not give correct
## centroids for longitude/latitude data
head(meta_points)

Se procede a graficar el mapa del departamento del Meta con sus respectivas divisiones políticas numeradas:

#ggplot(meta) + 
# geom_sf() +  
# geom_sf(data = meta_points, fill= "lightskyblue") +
# geom_text(data = meta_points, aes(x=X, y=Y, label= n_cpi), size = 3) +
# coord_sf(xlim = c(-75, -71), ylim = c(1.5, 5), expand = FALSE)

Mapa del Meta

Ahora, se procede a llamar a R el paquete conocido como scales. Este permite realizar mapas con un nivel de gráficas más estético, entre otras características.

library(scales)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
#ggplot(meta) + 
#  geom_sf(data=meta_points, aes(x=X, y=Y, fill =
#                                      n_cpi), color = "black", size = 0.25) +
#  geom_text(data = meta_points, aes(x=X, y=Y,label = n_cpi), size = 2) +
#  theme(aspect.ratio=1)+
#  scale_fill_distiller(name="N_MPIO", palette = "YlOrRd", breaks = pretty_breaks(n = 5))+
#  labs(title="Mapa del Meta")

Mapa del Meta

Este mapa puede ser guardado como un archivo .PDF o PNG, así:

ggsave("meta_municipios.pdf")
## Saving 7 x 5 in image
ggsave("meta_municipios_2.png", width = 6, height = 6, dpi = "screen")

Usando Leaflet para visualizar datos

Para realizar esta y última parte de este cuaderno, es necesario instalar los paquetes leaflet y lwgeom.

#install.packages("leaflet")
#  install.packages("lwgeom")

Llamando a la interfaz de R a leaflet y Lwgeom:

library(leaflet)
library(lwgeom)
## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.8.0, PROJ 6.3.1

Para poder usar la librería leaflet, es necesario convertir los datos meta_points de simple features a spatial points:

met_points <- as(meta_points, 'Spatial')

Se procede a calcular el area en metros cuadrados de cada municipio:

mun_meta$area <- st_area(mun_meta)

Ahora se crea un campo para almacenar el area en kilometros cuadrados:

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

Llamando a la interfaz la salida mun_meta$km2

mun_meta$km2
## Units: [m^2]
##  [1]  1285.7843  1123.6871   403.6742   913.3079   514.6387  1157.8556
##  [7]   624.9786   277.8402   568.2108   117.3389   570.9498   348.2999
## [13]   595.6985  2280.6239  6436.9091   815.4243  1242.2325  6908.8026
## [19]  2535.2038  3427.1841   365.0736   807.6292  1169.0729   237.3381
## [25]  5931.2293  4818.8950 10839.0233 11931.3843 17217.9846

Nuevamente, es necesario convertir los datos de simple features a spatial plygons

met_mun <- as(mun_meta, 'Spatial')

Ahora, se procede a preparar la gráfica:

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


labels <- mun_meta$MPIO_CNMBR

labels
##  [1] "VILLAVICENCIO"        "ACACÍAS"              "BARRANCA DE UPIA"    
##  [4] "CABUYARO"             "CASTILLA LA NUEVA"    "CUBARRAL"            
##  [7] "CUMARAL"              "EL CALVARIO"          "EL CASTILLO"         
## [10] "EL DORADO"            "FUENTE DE ORO"        "GRANADA"             
## [13] "GUAMAL"               "MESETAS"              "URIBE"               
## [16] "LEJANÍAS"             "PUERTO CONCORDIA"     "PUERTO LÓPEZ"        
## [19] "PUERTO LLERAS"        "PUERTO RICO"          "RESTREPO"            
## [22] "SAN CARLOS DE GUAROA" "SAN JUAN DE ARAMA"    "SAN JUANITO"         
## [25] "SAN MARTÍN"           "VISTAHERMOSA"         "LA MACARENA"         
## [28] "MAPIRIPÁN"            "PUERTO GAITÁN"

Ahora, aprovechando las características de las librerías decargadas, se procede a graficar el mapa que permite saber el nombre del municipio interactivamente:

m <- leaflet(met_mun) %>%
  setView(-73, 3, 7.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")

Llamando a m para ver el gráfico:

#m

Mapa del Meta

De igual forma, usando diferentes proveedores de mapas, se pueden graficar objetos con características interesantes, a continuación se presenta una serie de ejemplos: Ejemplo de representación sencilla, de proveedor Esri.WorldImagery.

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

Mapa del Meta

Ejemplo de mapa con el proveedor Esri.NatGeoWorldMap

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

Mapa del Meta

Ejemplo de mapa con el proveedor Stamen.Toner.

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

Mapa del Meta

Ejemplo de mapa con el proveedor CartoDB.Positron.

#leaflet() %>%
#  addProviderTiles(providers$CartoDB.Positron, options = providerTileOptions( opacity = 0.99)) %>%
#  addPolygons(data = met_mun, popup = met_mun$MPIO_CNMBR, fillColor = ~pal(km2) ,
#    weight = 2,
#    opacity = 2,
#    color = "skyblue",
#    dashArray = "3",
#    fillOpacity = 0.2,
#    highlight = highlightOptions(
#      weight = 5,
#      color = "green",
#      dashArray = "",
#      fillOpacity = 0.7,
#      bringToFront = TRUE),
#    label = labels, stroke = TRUE, smoothFactor = 0.25)

Mapa del Meta

Un último ejemplo con el mapa del proveedor Wikimedia.

#leaflet() %>%
#  addProviderTiles(providers$Wikimedia, options = providerTileOptions( opacity = 0.99)) %>%
#  addPolygons(data = met_mun, popup = met_mun$MPIO_CNMBR, fillColor = ~pal(km2) ,
#    weight = 2,
#    opacity = 2,
#    color = "lightyellow",
#    dashArray = "3",
#    fillOpacity = 0.2,
#    highlight = highlightOptions(
#      weight = 5,
#      color = "red",
#      dashArray = "",
#      fillOpacity = 0.7,
#      bringToFront = TRUE),
#    label = labels, stroke = TRUE, smoothFactor = 0.25)

Mapa del Meta

Con los ejemplos anteriores, se intentó demostrar la versatilidad de programa R-studio y su utilidad en análisis de datos geoespaciales, apliccados en este caso, para el departamento del Meta en Colombia.