1. Introduccion

En este cuaderno de Rmarkdown realizaremos diferentes tipos de mapas tematicos del departamento de Cundinamarca.

rm(list=ls())

Instalaremos los siguientes paquetes estadisticos:

list.of.packages <- c("tidyverse", "rgeos", "sf", "raster", "cartography", "SpatialPosition")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(tidyverse)
library(readxl) ## Libreria de tidyverse
library(rgeos)
library(raster)
library(sf)
library(cartography)
library(SpatialPosition)

2. Leer datos

Descargaremos la base de datos del DANE, de las Necesidades Basicas Insatisfechas (NBI) del pais, y posteriormente vamos a leer el archivo de excel de la forma “xlsx”, utilizando la funcion “read_excel()”, estos datos fueron previamente modificados eliminando filas combinadas de la base de datos que el programa de R no puede leer y filtrando solo la informacion del departamento de Cundinamarca

nbi <- read_excel("C:/Users/mario/Documents/Ing. Agronomica/9 Semestre/Geomatica/Informes tecnicos/Mapa_Tematico/NBI_CUND.xlsx")
nbi
## # A tibble: 116 x 12
##    COD_DEPTO DEPTO COD_MUN CODIGO MUNICIPIO   NBI MISERIA VIVIENDA SERVICIOS
##    <chr>     <chr> <chr>   <chr>  <chr>     <dbl>   <dbl>    <dbl>     <dbl>
##  1 25        CUND~ 001     25001  AGUA DE ~  7.86   1.13     2.44      0.825
##  2 25        CUND~ 019     25019  ALBÁN      9.80   1.11     1.89      1.49 
##  3 25        CUND~ 035     25035  ANAPOIMA  10.4    1.36     3.86      1.10 
##  4 25        CUND~ 040     25040  ANOLAIMA   9.78   1.49     2.66      1.06 
##  5 25        CUND~ 053     25053  ARBELÁEZ  10.0    1.25     3.03      0.146
##  6 25        CUND~ 086     25086  BELTRÁN    9.66   0.465    1.92      1.40 
##  7 25        CUND~ 095     25095  BITUIMA    9.67   2.02     1.98      2.70 
##  8 25        CUND~ 099     25099  BOJACÁ     5.37   0.259    0.673     0.186
##  9 25        CUND~ 120     25120  CABRERA   10.1    0.820    1.41      1.87 
## 10 25        CUND~ 123     25123  CACHIPAY   9.24   1.29     3.02      0.540
## # ... with 106 more rows, and 3 more variables: HACINAMIENTO <dbl>,
## #   INASISTENCIA <dbl>, ECONOMIA <dbl>

Visualizacion del encabezado de datos:

head(nbi)
## # A tibble: 6 x 12
##   COD_DEPTO DEPTO COD_MUN CODIGO MUNICIPIO   NBI MISERIA VIVIENDA SERVICIOS
##   <chr>     <chr> <chr>   <chr>  <chr>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 25        CUND~ 001     25001  AGUA DE ~  7.86   1.13      2.44     0.825
## 2 25        CUND~ 019     25019  ALBÁN      9.80   1.11      1.89     1.49 
## 3 25        CUND~ 035     25035  ANAPOIMA  10.4    1.36      3.86     1.10 
## 4 25        CUND~ 040     25040  ANOLAIMA   9.78   1.49      2.66     1.06 
## 5 25        CUND~ 053     25053  ARBELÁEZ  10.0    1.25      3.03     0.146
## 6 25        CUND~ 086     25086  BELTRÁN    9.66   0.465     1.92     1.40 
## # ... with 3 more variables: HACINAMIENTO <dbl>, INASISTENCIA <dbl>,
## #   ECONOMIA <dbl>

Podemos visualizar de los datos que llamamos “nbi” el valor maximo utilizando la funcion “which.max” y “slice()” para elegir filas por su posición ordinal en la tabla

nbi %>% 
    slice(which.max(NBI)) -> max_nbi

max_nbi
## # A tibble: 1 x 12
##   COD_DEPTO DEPTO COD_MUN CODIGO MUNICIPIO   NBI MISERIA VIVIENDA SERVICIOS
##   <chr>     <chr> <chr>   <chr>  <chr>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 25        CUND~ 885     25885  YACOPÍ     40.7    11.3     30.9      9.41
## # ... with 3 more variables: HACINAMIENTO <dbl>, INASISTENCIA <dbl>,
## #   ECONOMIA <dbl>

De igual forma podemos buscar el valor minimo, con “which.min()”

nbi %>% 
    slice(which.min(NBI)) -> min_nbi

min_nbi
## # A tibble: 1 x 12
##   COD_DEPTO DEPTO COD_MUN CODIGO MUNICIPIO   NBI MISERIA VIVIENDA SERVICIOS
##   <chr>     <chr> <chr>   <chr>  <chr>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 25        CUND~ 758     25758  SOPÓ       2.83  0.0807    0.125    0.0686
## # ... with 3 more variables: HACINAMIENTO <dbl>, INASISTENCIA <dbl>,
## #   ECONOMIA <dbl>

Podemos aplicar lo mismo para mirar cualquier varibale de los datos, por ejemplo “MISERIA”:

nbi %>% 
    slice(which.min(MISERIA)) -> min_nbi

min_nbi
## # A tibble: 1 x 12
##   COD_DEPTO DEPTO COD_MUN CODIGO MUNICIPIO   NBI MISERIA VIVIENDA SERVICIOS
##   <chr>     <chr> <chr>   <chr>  <chr>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 25        CUND~ 200     25200  COGUA      3.57  0.0529    0.135   0.00481
## # ... with 3 more variables: HACINAMIENTO <dbl>, INASISTENCIA <dbl>,
## #   ECONOMIA <dbl>

Ordear datos, utilizando la funcion “arrange()” y de forma descendente “desc()”, con esto podremos visualizar el municipio con mayor NBI del departamento de Cundinamarca:

nbi %>% 
  arrange(desc(NBI))  -> desc_nbi

desc_nbi
## # A tibble: 116 x 12
##    COD_DEPTO DEPTO COD_MUN CODIGO MUNICIPIO   NBI MISERIA VIVIENDA SERVICIOS
##    <chr>     <chr> <chr>   <chr>  <chr>     <dbl>   <dbl>    <dbl>     <dbl>
##  1 25        CUND~ 885     25885  YACOPÍ     40.7   11.3     30.9       9.41
##  2 25        CUND~ 518     25518  PAIME      34.7    9.13    16.1      18.0 
##  3 25        CUND~ 398     25398  LA PEÑA    29.9    8.40    11.2      18.6 
##  4 25        CUND~ 148     25148  CAPARRAPÍ  29.1    4.70    16.5       5.38
##  5 25        CUND~ 823     25823  TOPAIPÍ    27.4    7.72    18.1       9.62
##  6 25        CUND~ 258     25258  EL PEÑÓN   25.7    3.79    15.6       1.68
##  7 25        CUND~ 862     25862  VERGARA    22.3    3.05    13.8       2.67
##  8 25        CUND~ 839     25839  UBALÁ      21.3    4.17     3.71      6.96
##  9 25        CUND~ 438     25438  MEDINA     20.6    3.25    13.2       1.28
## 10 25        CUND~ 394     25394  LA PALMA   20.2    4.58    13.0       4.11
## # ... with 106 more rows, and 3 more variables: HACINAMIENTO <dbl>,
## #   INASISTENCIA <dbl>, ECONOMIA <dbl>

6. JOIN

Ahora vamos a leer los datos espaciales del formato “shapfile” del departamento Cundinamarca que se descargaron previamente de la base de datos del geoportal del DANE

munic <-st_read("C:/Users/mario/Documents/Ing. Agronomica/9 Semestre/Geomatica/Cundinamarca/DANE/25_CUNDINAMARCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\mario\Documents\Ing. Agronomica\9 Semestre\Geomatica\Cundinamarca\DANE\25_CUNDINAMARCA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 116 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -74.89063 ymin: 3.730129 xmax: -73.05256 ymax: 5.837258
## geographic CRS: WGS 84

Visualizacion del encabezado de datos del geoportal, donde podremos ver informacion importante para la elaboracion de nuestro posterior mapa tematico, como por ejemplo el codigo EPSG de coordenadas 4326 y las variables que contiene la base de datos

head(munic)
## Simple feature collection with 6 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -74.85582 ymin: 4.017652 xmax: -74.03929 ymax: 5.295559
## geographic CRS: WGS 84
##   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                       MPIO_CRSLC MPIO_NAREA
## 1         25      25483     NARIÑO                             1899   55.16263
## 2         25      25513      PACHO                             1604  402.62678
## 3         25      25506    VENECIA Decreto 727 Septiembre 5 de 1951  122.20332
## 4         25      25491    NOCAIMA                             1735   68.93823
## 5         25      25489    NIMAIMA    Ordenanza 30 de Julio de 1904   59.49982
## 6         25      25488       NILO                             1899  224.70968
##   MPIO_NANO   DPTO_CNMBR Shape_Leng  Shape_Area                       geometry
## 1      2017 CUNDINAMARCA  0.5219117 0.004493674 MULTIPOLYGON (((-74.7413 4....
## 2      2017 CUNDINAMARCA  1.2021630 0.032839624 MULTIPOLYGON (((-74.21183 5...
## 3      2017 CUNDINAMARCA  0.4873683 0.009951825 MULTIPOLYGON (((-74.44812 4...
## 4      2017 CUNDINAMARCA  0.3987992 0.005621814 MULTIPOLYGON (((-74.39011 5...
## 5      2017 CUNDINAMARCA  0.4990524 0.004852617 MULTIPOLYGON (((-74.34091 5...
## 6      2017 CUNDINAMARCA  0.8442993 0.018305852 MULTIPOLYGON (((-74.56068 4...

Ahora procederemos a revisar la estructura de los datos que tenemos en comun de los datos de NBI y los datos espaciales de DANE del departamento de Cundinamarca, para este caso utilizaremos las variables “MPIO_CCDGO” y “CODIGO” que contienen el mismo valor por el cual podremos unir las dos bases de datos posteriormente

class(nbi$CODIGO)
## [1] "character"
class(munic$MPIO_CCDGO)
## [1] "factor"

Como podemos ver la clase de las variables que vamos a unir son diferentes, la variable “CODIGO” es un character y la variable “MPIO_CCDGO” un factor, debemos transformarlas para que sean de la misma clase, lo que podemos realizar con la funcion “factor()”, o “as.factor()”

nbi$CODIGO = factor(nbi$CODIGO)

Teniendo las dos condiciones para unir las dos bases de datos: 1. Una variable en comun 2. Que las dos variables tengan la misma clase

Procederemos a realizar la union (JOIN) utilizando la funcion “left_join()” de los datos espaciales “munic” y los de NBI “nbi”, en donde la funcion “by()” es la que nos indica cual es la variable en comun a UNIR

nbi_munic = left_join(munic, nbi, by=c("MPIO_CCDGO"="CODIGO"))

Vamos a seleccionar datos con la funcion “select()” de la libreria “dplyr”

nbi_munic %>%
  dplyr::select(MUNICIPIO, MPIO_CCDGO, NBI)  ->  check_nbi_munic

head(check_nbi_munic)
## Simple feature collection with 6 features and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -74.85582 ymin: 4.017652 xmax: -74.03929 ymax: 5.295559
## geographic CRS: WGS 84
##   MUNICIPIO MPIO_CCDGO       NBI                       geometry
## 1    NARIÑO      25483  9.809524 MULTIPOLYGON (((-74.7413 4....
## 2     PACHO      25513  7.984099 MULTIPOLYGON (((-74.21183 5...
## 3   VENECIA      25506 13.102776 MULTIPOLYGON (((-74.44812 4...
## 4   NOCAIMA      25491 10.666153 MULTIPOLYGON (((-74.39011 5...
## 5   NIMAIMA      25489 12.847336 MULTIPOLYGON (((-74.34091 5...
## 6      NILO      25488  7.916667 MULTIPOLYGON (((-74.56068 4...

Para revisar que el “JOIN” se realizo correctamente revisamos que los datos esten completos y no salgan datos con la notacion (NA), para este caso se puede dar la funcion “view()” para visualizar todos los datos y se comparan con el “JOIN” que se realizo.

Ahora transformaremos las coordenadas del JOIN que estan con el codigo EPSG:4326, perteneciente al elipsoide WGS 84 – WGS84 - World Geodetic System 1984 que tiene coordenadas geograficas, a el codigo EPSG: 3116, MAGNA-SIRGAS que es el sistema de referencia del pais, cuya zona de origen es Bogota y tiene coordenadas planas. Esta transformacion la realizamos debidoa que la libreria que utilizaremos posteriormente “cartography” las requiere.

nbi_munic_new <- st_transform(nbi_munic, crs = 3116)

7. Mapa Tematico

Utilizaremos la libreria “cartography” que tiene como objetivo obtener mapas temáticos con la calidad visual de aquellos construidos con un software clásico de mapeo o SIG. Sus funciones son intuitivas para los usuarios y aseguran la compatibilidad con los flujos de trabajo comunes de R, para esta libreria el formato a utilizar para objetos espaciales es “sf”

7.1 Mapa base de OpenStreetMap y símbolos proporcionales

Las funciones “getTiles()” y “tilesLayer()” descargan y muestran los mosaicos de OpenStreetMap, y “propSymbolsLayer()” muestra símbolos con áreas proporcionales a una variable cuantitativa (por ejemplo, NBI). Hay varios símbolos disponibles (círculos, cuadrados, barras).

# Descarga
mun.osm <- getTiles(
x = nbi_munic_new, 
type = "OpenStreetMap", 
zoom = 10,
cachedir = TRUE, ### Guardar en memoria temporal
crop = FALSE
)
# Definir margenes 
opar <- par(mar = c(0,0,1.2,0))
# Plot 
tilesLayer(x = mun.osm)
# Plot Municipios
plot(st_geometry(nbi_munic_new), col = NA, border = "black", add=TRUE, bg = "white")
# Grafica NBI
propSymbolsLayer(
  x = nbi_munic_new, 
  var = "NBI", 
  inches = 0.14, 
  col = "firebrick",
  legend.pos = "topright",  
  legend.title.txt = "Total NBI",
)
# Capa de Diseño con informacion en pie de pagina, como escala, autor, fuentes, entre otros.
layoutLayer(title = " Distribucion NBI Cundinamarca", 
            sources = "Sources: DANE, 2018\n© OpenStreetMap",
            author = " Mario Celis ",
            frame = TRUE, north = FALSE, tabtitle = TRUE)
# Flecha del Norte
north(pos = "topleft")

7.2 Mapa de Coropletas

En los mapas coropléticos, las áreas se sombrean de acuerdo con la variación de una variable cuantitativa. Se utilizan para representar razones o índices. La funcion “choroLayer()” de la libreria “Cartography” muestra mapas de coropletas. Los argumentos nclass, method y breaks permiten personalizar la clasificación de la variable, la funcion “getBreaks()” permite clasificar fuera de la función misma. Las paletas de colores se definen con “col” y se puede crear un conjunto de colores con “carto.pal()”

# Establecer margenes
opar <- par(mar = c(0,0,1.2,0))
#Establecer el color de fondo de la figura
par(bg="red")
# Plot Municipios 
plot(st_geometry(nbi_munic_new), col = NA, border = NA, bg = "white")
# Plot NBI (Mapa coropletas)
choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "orange.pal", n1 = 5),
  border = "black", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "NBI",
  add = TRUE
) 
# Capa de Diseño con informacion en pie de pagina, como escala, autor, fuentes, entre otros
layoutLayer(title = "Distribucion NBI Cundinamarca", 
            sources = "Source: DANE, 2018",
            author = "Mario Alejandro Celis", 
            frame = TRUE, north = TRUE, tabtitle = TRUE, col="black") 
# Flecha del Norte
north(pos = "topleft")

7.3 Símbolos proporcionales y mapa de tipología

La función propSymbolsTypoLayer () crea un mapa de símbolos que son proporcionales a los valores de una primera variable y coloreados para reflejar las modalidades de una segunda variable cualitativa. Se utiliza una combinación de argumentos propSymbolsLayer () y typoLayer ().

Primero, necesitamos crear una variable cualitativa. Usemos la función “mutate()” para esta tarea.

nbi_munic_2 <- dplyr::mutate(nbi_munic_new, Pobreza = ifelse(MISERIA > 10, "Extrema", 
                                                         ifelse(HACINAMIENTO > 3, "Alta", "Intermedia")))

Tenga en cuenta que el nuevo atributo se llama pobreza. Su valor depende de los valores umbral definidos por el usuario. Hay que verificar el condicional “ifelse()” para comprender lo que está sucediendo con la nueva variable “Pobreza”.

head(nbi_munic_2)
## Simple feature collection with 6 features and 21 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 913602.2 ymin: 936036.9 xmax: 1004237 ymax: 1077338
## projected CRS:  MAGNA-SIRGAS / Colombia Bogota zone
##   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                       MPIO_CRSLC MPIO_NAREA
## 1         25      25483     NARIÑO                             1899   55.16263
## 2         25      25513      PACHO                             1604  402.62678
## 3         25      25506    VENECIA Decreto 727 Septiembre 5 de 1951  122.20332
## 4         25      25491    NOCAIMA                             1735   68.93823
## 5         25      25489    NIMAIMA    Ordenanza 30 de Julio de 1904   59.49982
## 6         25      25488       NILO                             1899  224.70968
##   MPIO_NANO   DPTO_CNMBR Shape_Leng  Shape_Area COD_DEPTO        DEPTO COD_MUN
## 1      2017 CUNDINAMARCA  0.5219117 0.004493674        25 CUNDINAMARCA     483
## 2      2017 CUNDINAMARCA  1.2021630 0.032839624        25 CUNDINAMARCA     513
## 3      2017 CUNDINAMARCA  0.4873683 0.009951825        25 CUNDINAMARCA     506
## 4      2017 CUNDINAMARCA  0.3987992 0.005621814        25 CUNDINAMARCA     491
## 5      2017 CUNDINAMARCA  0.4990524 0.004852617        25 CUNDINAMARCA     489
## 6      2017 CUNDINAMARCA  0.8442993 0.018305852        25 CUNDINAMARCA     488
##   MUNICIPIO       NBI   MISERIA VIVIENDA SERVICIOS HACINAMIENTO INASISTENCIA
## 1    NARIÑO  9.809524 1.4761905 1.428571 1.4761905     1.380952    0.2380952
## 2     PACHO  7.984099 0.9895547 2.600753 1.2517444     1.649258    0.8753753
## 3   VENECIA 13.102776 1.5830492 2.703361 1.3395032     4.846566    1.0716025
## 4   NOCAIMA 10.666153 2.1948402 4.062380 2.1948402     2.445129    1.3477089
## 5   NIMAIMA 12.847336 1.5691402 5.230467 1.7979732     2.615234    0.9480222
## 6      NILO  7.916667 1.2356322 2.313218 0.1867816     2.686782    1.3936782
##   ECONOMIA                       geometry    Pobreza
## 1 6.761905 MULTIPOLYGON (((926329 9864... Intermedia
## 2 2.740305 MULTIPOLYGON (((985111.1 10... Intermedia
## 3 4.724793 MULTIPOLYGON (((958849.2 94...       Alta
## 4 3.118983 MULTIPOLYGON (((965338.7 10... Intermedia
## 5 3.824779 MULTIPOLYGON (((970798 1066... Intermedia
## 6 2.643678 MULTIPOLYGON (((946369 9754... Intermedia
library(sf)
library(cartography)
#Establecer márgenes
opar <- par(mar = c(0,0,1.2,0))
par(bg="red")
# Plot de Municipios
plot(st_geometry(nbi_munic_2), col="lightcyan2", border="black", bg = "white", 
     lwd = 0.5)
# Plot simbolos con coloracion de Coropletas
propSymbolsTypoLayer(
  x = nbi_munic_2, 
  var = "NBI", 
  inches = 0.175,
  symbols = "square",
  border = "black",
  lwd = .4,
  legend.var.pos = c(1170000,1077338), 
  legend.var.title.txt = "NBI",
  var2 = "Pobreza",
  legend.var2.values.order = c("Extrema",
                               "Intermedia",
                               "Alta"),
  col = carto.pal(pal1 = "grey.pal", n1 = 3),
  legend.var2.pos = c(1170000, 977338), 
  legend.var2.title.txt = "Pobreza"
) 
# Capa de Diseño con informacion en pie de pagina, como escala, autor, fuentes, entre otros
layoutLayer(title="Distribucion NBI en Cundinamarca", 
            author = "Mario Alejandro Celis", 
            sources = "Source: DANE, 2018", 
            scale = 50, posscale = c(1175000, 905338), tabtitle = TRUE, frame = TRUE)
# Flecha del Norte
north(pos = "topleft")

7.4 Mapa con etiquetas

Este mapa utiliza la funcion “choroLayer()” para darle color a las areas municipales y la funcion “lableLayer()” para añadir las etiquetas.

library(sf)
library(cartography)
#Establecer márgenes
opar <- par(mar = c(0,0,1.2,0))
# Establecer color de fondo de la figura
par(bg ="red")
# Plot Municipios
plot(st_geometry(nbi_munic_2), col = "white", border = "darkseagreen4", 
     bg = "white", lwd = 0.5)
# Plot NBI
choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "sand.pal", n1 = 5),
  border = "black", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "NBI",
  add = TRUE
) 
# Plot de etiquetas
labelLayer(
  x = nbi_munic_2, 
  txt = "MUNICIPIO", 
  col= "white", 
  cex = 0.45, 
  font = 4,
  halo = TRUE, 
  bg = "grey25", 
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)
# Capa de Diseño con informacion en pie de pagina, como escala, autor, fuentes, entre otros
layoutLayer(
  title = "Municipios de Cundinamarca",
  sources = "Source: DANE, 2018",  
  author = "Mario Alejandro Celis", 
  frame = TRUE,
  north = TRUE, 
  tabtitle = TRUE, 
  theme = "taupe.pal",
  scale = 50
)

7.5 Mapas de Isopletas

Los mapas de Isopletas se basan en el supuesto de que el fenómeno a representar tiene una distribución continua. Estos mapas utilizan un enfoque de modelado de interacción espacial que tiene como objetivo calcular indicadores basados en valores de stock ponderados por distancia. Permite una representación espacial del fenómeno independiente de la heterogeneidad inicial de la división territorial. La funcion “smoothLayer()” depende en gran medida del paquete “SpatialPosition” y utiliza una capa de puntos marcados y un conjunto de parámetros (una función de interacción espacial y sus parámetros) y muestra una capa de mapa isopleta.

Usemos otro conjunto de datos para hacer un mapa isopleta. En este caso, subiremos datos estadísticos sobre la producción de café de 2018 en Cundinamarca. Lea el conjunto de datos

crops2018 <- read_excel("C:/Users/mario/Documents/Ing. Agronomica/9 Semestre/Geomatica/Informes tecnicos/Mapa_Tematico/EVA_Cund2018.xlsx")
head(crops2018)
## # A tibble: 6 x 17
##   COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO SUBGRUPO CULTIVO DESAGREGACION_R~
##     <dbl> <chr>          <dbl> <chr>     <chr> <chr>    <chr>   <chr>           
## 1      25 CUNDINAMARCA   25214 COTA      HORT~ ACELGA   ACELGA  ACELGA          
## 2      25 CUNDINAMARCA   25754 SOACHA    HORT~ ACELGA   ACELGA  ACELGA          
## 3      25 CUNDINAMARCA   25175 CHIA      HORT~ ACELGA   ACELGA  ACELGA          
## 4      25 CUNDINAMARCA   25126 CAJICA    HORT~ ACELGA   ACELGA  ACELGA          
## 5      25 CUNDINAMARCA   25799 TENJO     HORT~ ACELGA   ACELGA  ACELGA          
## 6      25 CUNDINAMARCA   25312 GRANADA   FRUT~ FRUTALE~ AGRAZ   AGRAZ           
## # ... with 9 more variables: YEAR <dbl>, PERIODO <chr>, Area_Siembra <dbl>,
## #   Area_Cosecha <dbl>, Produccion <dbl>, Rendimiento <dbl>, ESTADO <chr>,
## #   N_CIENTIFICO <chr>, CICLO_CULTIVO <chr>

Filtrar el cultivo de interes (Cafe)

crops2018 %>%
  filter(CULTIVO == "CAFE") -> cafe2018
cafe2018$TEMP <-  as.character(cafe2018$COD_MUN)
cafe2018$MPIO_CCDGO <- as.factor(cafe2018$TEMP)
cafe_munic = left_join(munic, cafe2018, by="MPIO_CCDGO")
## Warning: Column `MPIO_CCDGO` joining factors with different levels, coercing to
## character vector
head(cafe_munic)
## Simple feature collection with 6 features and 27 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -74.85582 ymin: 4.017652 xmax: -74.03929 ymax: 5.295559
## geographic CRS: WGS 84
##   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                       MPIO_CRSLC MPIO_NAREA
## 1         25      25483     NARIÑO                             1899   55.16263
## 2         25      25513      PACHO                             1604  402.62678
## 3         25      25506    VENECIA Decreto 727 Septiembre 5 de 1951  122.20332
## 4         25      25491    NOCAIMA                             1735   68.93823
## 5         25      25489    NIMAIMA    Ordenanza 30 de Julio de 1904   59.49982
## 6         25      25488       NILO                             1899  224.70968
##   MPIO_NANO   DPTO_CNMBR Shape_Leng  Shape_Area COD_DEP DEPARTAMENTO COD_MUN
## 1      2017 CUNDINAMARCA  0.5219117 0.004493674      NA         <NA>      NA
## 2      2017 CUNDINAMARCA  1.2021630 0.032839624      25 CUNDINAMARCA   25513
## 3      2017 CUNDINAMARCA  0.4873683 0.009951825      25 CUNDINAMARCA   25506
## 4      2017 CUNDINAMARCA  0.3987992 0.005621814      25 CUNDINAMARCA   25491
## 5      2017 CUNDINAMARCA  0.4990524 0.004852617      25 CUNDINAMARCA   25489
## 6      2017 CUNDINAMARCA  0.8442993 0.018305852      25 CUNDINAMARCA   25488
##   MUNICIPIO             GRUPO SUBGRUPO CULTIVO
## 1      <NA>              <NA>     <NA>    <NA>
## 2     PACHO OTROS PERMANENTES     CAFE    CAFE
## 3   VENECIA OTROS PERMANENTES     CAFE    CAFE
## 4   NOCAIMA OTROS PERMANENTES     CAFE    CAFE
## 5   NIMAIMA OTROS PERMANENTES     CAFE    CAFE
## 6      NILO OTROS PERMANENTES     CAFE    CAFE
##   DESAGREGACION_REGIONAL_SISTEMA_PRODUCTIVO YEAR PERIODO Area_Siembra
## 1                                      <NA>   NA    <NA>           NA
## 2                                      CAFE 2018    2018          834
## 3                                      CAFE 2018    2018          409
## 4                                      CAFE 2018    2018           41
## 5                                      CAFE 2018    2018           35
## 6                                      CAFE 2018    2018          585
##   Area_Cosecha Produccion Rendimiento                 ESTADO   N_CIENTIFICO
## 1           NA         NA          NA                   <NA>           <NA>
## 2          694        837         121 CAFE VERDE EQUIVALENTE COFFEA ARABICA
## 3          370        255          69 CAFE VERDE EQUIVALENTE COFFEA ARABICA
## 4           37         54         147 CAFE VERDE EQUIVALENTE COFFEA ARABICA
## 5           34         46         138 CAFE VERDE EQUIVALENTE COFFEA ARABICA
## 6          504        489          97 CAFE VERDE EQUIVALENTE COFFEA ARABICA
##   CICLO_CULTIVO  TEMP                       geometry
## 1          <NA>  <NA> MULTIPOLYGON (((-74.7413 4....
## 2    PERMANENTE 25513 MULTIPOLYGON (((-74.21183 5...
## 3    PERMANENTE 25506 MULTIPOLYGON (((-74.44812 4...
## 4    PERMANENTE 25491 MULTIPOLYGON (((-74.39011 5...
## 5    PERMANENTE 25489 MULTIPOLYGON (((-74.34091 5...
## 6    PERMANENTE 25488 MULTIPOLYGON (((-74.56068 4...
rep_cafe <- st_transform(cafe_munic, crs = 3116)
#Establecer márgenes
opar <- par(mar = c(0,0,1.2,0))
par(bg = "red")
# Plot Municipios 
plot(st_geometry(rep_cafe), col="yellowgreen", border = "black", bg = "white")
# Plot Mapa de Isopletas
smoothLayer(
  x = rep_cafe, 
  var = 'Rendimiento',
  typefct = "exponential",
  span = 60000,
  beta = 2,
  nclass =6,
  col = carto.pal(pal1 = 'green.pal', n1 = 6),
  border = "black",
  lwd = 0.1, 
  mask = rep_cafe, 
  legend.values.rnd = -3,
  legend.title.txt = "Produccion",
  legend.pos = "topright", 
  add=TRUE
)
# Anotacion en el mapa
text(x = 850000, y = 977338, cex = 0.6, adj = 0, font = 3,  labels = 
       "Funcion de distancia:\n- Tipo = Exponencial\n- Beta = 2\n- Span = 20 km")
# Capa de Diseño con informacion en pie de pagina, como escala, autor, fuentes, entre otros
layoutLayer(title = "Distribucion de la Produccion de Cafe en Cundinamarca",
            sources = "Sources: DANE and MADR, 2018",
            author = "Mario Alejandro Celis",
            frame = FALSE, north = FALSE, tabtitle = TRUE, theme = "white", scale = 50)
# Flecha del NORTE
north(pos = "topleft")

8. Guardando mapas

Vamos a realizar otro mapa de producción de café en 2018. Esta vez usaremos símbolos proporcionales y mapas de coropletas. La salida se guardará como un archivo .png. La funcion “propSymbolsChoroLayer()” crea un mapa de símbolos que son proporcionales a los valores de una primera variable y coloreados para reflejar la clasificación de una segunda variable. - se utiliza una combinación de argumentos “propSymbolsLayer()” y “choroLayer()”.

# Establecer margenes
opar <- par(mar = c(0,0,1.2,0))
# Plot de los Municipios
plot(st_geometry(rep_cafe), col="yellowgreen", border="darkseagreen4",  
     bg = "white", lwd = 0.6)
# Plot simbolos con coloracion coropletas
propSymbolsChoroLayer(x = rep_cafe, var = "Produccion", var2 = "Rendimiento",
                      col = carto.pal(pal1 = "green.pal", n1 = 3,
                                      pal2 = "red.pal", n2 = 3),
                      inches = 0.25, method = "q6",
                      border = "grey50", lwd = 1,
                      legend.title.cex = 0.5,
                      legend.values.cex = 0.5,
                      legend.var.pos = "right", 
                      legend.var2.pos = "left",
                      legend.var2.values.rnd = 2,
                      legend.var2.title.txt = "Rendimiento\n(Ton/Ha)",
                      legend.var.title.txt = "Produccion de Cafe en 2018",
                      legend.var.style = "e")
# Plot etiquetas
labelLayer(
  x = rep_cafe, 
  txt = "MPIO_CNMBR", 
  col= "black", 
  cex = 0.2, 
  font = 4,
  halo = FALSE,
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)
# Capa de Diseño con informacion en pie de pagina, como escala, autor, fuentes, entre otros
layoutLayer(title="Produccion y Rendimiento de Cafe en Cundinamaca, 2018",
            author = "Mario Alejandro Celis", 
            sources = "Sources: MADR & DANE, 2018", 
            scale = 50, posscale = c(1175000, 905338), tabtitle = FALSE, frame = TRUE)
# Flecha del Norte
north(pos = "topleft")

Utilizando la funcion “png()” podemos guardar el mapa en formato png en la direccion que deseemos, al realizar esto el mapa se guardara directamente y no se vera en el editor del RStudio, se debe tener muy en cuenta que el tamaño de los simbolos y letras cuando se pase a formato png se debe modificar, aumentando la gram mayoria de items para que se vean adecuadamente en el mapa final, notese estos cambios de tamaños en el siguiente codigo:

png("C:/Users/mario/Documents/Ing. Agronomica/9 Semestre/Geomatica/Informes tecnicos/Mapa_Tematico.png", width = 2048, height = 1526)
opar <- par(mar = c(0,0,5,5))
# Plot de los Municipios
plot(st_geometry(rep_cafe), col="yellowgreen", border="darkseagreen4",  
     bg = "white", lwd = 0.6)
# Plot simbolos con coloracion coropletas
propSymbolsChoroLayer(x = rep_cafe, var = "Produccion", var2 = "Rendimiento",
                      col = carto.pal(pal1 = "green.pal", n1 = 3,
                                      pal2 = "red.pal", n2 = 3),
                      inches = 1.0, method = "q6",
                      border = "grey50", lwd = 1,
                      legend.title.cex = 2.2,
                      legend.values.cex = 2.2,
                      legend.var.pos = "right", 
                      legend.var2.pos = "left",
                      legend.var2.values.rnd = 2,
                      legend.var2.title.txt = "Rendimiento\n(Ton/Ha)",
                      legend.var.title.txt = "Produccion de Cafe en 2018",
                      legend.var.style = "e")
# Plot etiquetas
labelLayer(
  x = rep_cafe, 
  txt = "MPIO_CNMBR", 
  col= "black", 
  cex =1.6, 
  font = 4,
  halo = FALSE, 
  bg = "black", 
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)
# Capa de Diseño con informacion en pie de pagina, como escala, autor, fuentes, entre otros
layoutLayer(title="Producción y Rendimiento de Cafe en Cundinamarca, 2018",
            author = "Mario Alejandro Celis", 
            sources = "Fuentes: MADR & DANE, 2018", 
            scale = 50, tabtitle = FALSE, frame = TRUE)
#
title(main="Producción y Rendimiento de Cafe en Cundinamarca, 2018", cex.main=3.5,
      sub= "Fuente: MADR & DANE, 2018", cex.sub=3)
#
graticule = TRUE
# Flecha del NORTE
north(pos = "topleft")

9. Otro tipo de mapa

Mapa de tipologias y coloreado a lapiz

Utilizando la funcion “getPencilLayer()” transformamos POLYGONS o MULTIPOLYGONS en MULTILINESTRINGS. Esta función crea una capa que imita un dibujo a mano a lápiz, y la funcion “typoLayer()” muestra un mapa de tipología de una variable cualitativa.

library(sf)
library(cartography)
opar <- par(mar = c(0,0,1.2,0))
#Establecer el color de fondo de la figura
par(bg="red")
mtq_pencil <- getPencilLayer(
  x = nbi_munic_2, 
  size = 400, 
  lefthanded = F
)
# Plot Municipios
plot(st_geometry(nbi_munic_2), col = "white", border = "black", bg = "white")
# Plot Distribucion de Pobreza
typoLayer(
  x = mtq_pencil, 
  var="Pobreza",  
  col = c("aquamarine4", "yellow3","wheat"), 
  lwd = .7,
  legend.values.order = c("Extrema",
                               "Intermedia",
                               "Alta"),
  legend.pos = "topright",
  legend.title.txt = "Pobreza", 
  add = TRUE
)
# Plot etiquetas
labelLayer(
  x = nbi_munic_2, 
  txt = "MPIO_CNMBR", 
  col= "black", 
  cex = 0.45, 
  font = 4,
  halo = FALSE,
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)
#  Plot Municipios
plot(st_geometry(nbi_munic_2), lwd = 0.5, border = "grey20", add = TRUE, lty = 3)
# labels for a few  municipalities
labelLayer(x = nbi_munic_2[nbi_munic_2$Pobreza != "Simple municipality",], txt = "LIBGEO", 
           cex = 0.9, halo = TRUE, r = 0.15)
# Capa de Diseño con informacion en pie de pagina, como escala, autor, fuentes, entre otros
layoutLayer(title="Distribucion NBI en Cundinamarca", 
            author = "Mario Alejandro Celis", 
            sources = "Source: DANE, 2018", 
            north = FALSE, tabtitle = TRUE, postitle = "left", 
            col = "black", coltitle = "white",scale = 50) 
# Flecha del Norte
north(pos = "topleft")