Cartografia Temática

1.¿Por qué este cuaderno?

El siguiente cuaderno de R Markdown es elaborado en RStudio con el fin de ilustrar la cartografía temática en el departamento de Tolima.

2.Cartografía Temática

Los mapas temáticos son mapas basados en mapas topográficos que representan cualquier fenómeno geográfico de la superficie terrestre.Hacen referencia a la representación de ciertas características de objetos reales como lo puede ser cultivos, vegetación, suelos, entre otros; o de conceptos abstractos como indicadores de violencia, de desarrollo económico, etc.

Esto ayuda al usuario para mejorar su interpretación frente a distintas variables de manera facil y sin mayor dificultad.

2.1 Mapa coroplético

Es un tipo de mapa temático en el que las áreas se sombrean de distintos colores, frecuentemente de la misma gama cromática, que representan distintos valores de una variable anteriormente mencionada.

Un Ejemplo de mapa Coroplético es el siguiente:

Pueden encontrar más ejemplos en la pagina del DANE

2.2 Mapa de símbolos proporcionales o graduados

Un mapa de símbolos proporcionales, es un mapa de una determinada región, sobre el cual se indica determinada información cuantitativa sobre una cierta temática, asociada a determinados puntos del mapa. La información es presentada de manera gráfica utilizando símbolos cuya dimensión es proporcional a la magnitud de la cantidad representada.

Un Ejemplo de mapa es el siguiente:

Tomado de Minambiente

2.3 Mapa isarítmico o de contorno

Los mapas isarítmicos o de contorno usan isolíneas para representar valores continuos que cambian suavemente como los niveles de humedad. Estos mapas también pueden mostrar valores tridimensionales, como la elevación, en mapas topográficos. La ubicación de cada isolínea representa todos los lugares de la superficie que tienen un valor particular, también siguen la regla básica de que hay lados altos y bajos en relación con la isolínea.

Un Ejemplo de este mapa es:

Tomado de El Espectador

2.4 Mapa de puntos

Los mapas de puntos son una forma de detectar patrones espaciales o la distribución de datos sobre una región geográfica, colocando puntos de igual tamaño sobre una región geográfica. Estos se emplean especialmente para la representación de variables, tales como la población, desigualdad socioeconómica o la producción de un determinado cultivo, entre otros. Hay dos tipos de mapas de puntos, uno a uno (uno de los puntos representa un solo objeto) y uno a muchos (un punto representa una unidad en particular, por ejemplo, 1 punto = 10 personas).

Un Ejemplo de este mapa es:

3. Datos

Para realizar este cuaderno utilizamos datos de las Necesidades Básicas Insatisfechas (NBI) del Censo Nacional de Población y Vivienda 2018 que están disponibles en el Geoportal DANE para el departamento de Tolima.

Anteriormente, descargamos el archivo NBI, en formato .xlsx. Luego, utilizamos Excel para eliminar datos de municipios no ubicados dentro del Departamento de Tolima. Tambien eliminamos varias filas con imágenes institucionales o sin información. Se mantuvo las columnas que se refieren a todo el municipio (es decir, eliminar las columnas correspondientes a “cabecera” y “zona rural”). Guarde los datos resultantes, correspondientes solo a municipios de Tolima como NBI_tolima.xlsx.

4.Preparación

lo primero que haremos es limpiar la memoria

rm(list=ls())

Ahora, se instalarán las siguientes librerías

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)

A continuación, cargáremos las librerías

library(tidyverse)

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

## v ggplot2 3.3.0     v purrr   0.3.3
## v tibble  3.0.1     v dplyr   0.8.4
## v tidyr   1.0.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(readxl) 
library(rgeos)

## Loading required package: sp

## rgeos version: 0.5-2, (SVN revision 621)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 
##  Linking to sp version: 1.4-1 
##  Polygon checking: TRUE

library(raster)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:tidyr':
## 
##     extract

library(sf)

## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

library(cartography)
library(SpatialPosition)

5.Leer datos del NBI

Ahora leamos el archivo .xlsx con Necesidades Básicas Insatisfechas (NBI) del Censo Nacional de Población y Vivienda 2018 para el departamento de Tolima.

nbi <- read_excel("C:/Users/Pablo/Documents/PABLO/Geomatica/nbi/NBI_tolima.xlsx")

Vamos a ver que contiene el nbi.

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 73        TOLI~ 001     73001  IBAGUÉ     5.34   0.530    1.10      0.242
## 2 73        TOLI~ 024     73024  ALPUJARRA 12.0    0.688    0.933     1.40 
## 3 73        TOLI~ 026     73026  ALVARADO  13.9    2.03     2.89      1.07 
## 4 73        TOLI~ 030     73030  AMBALEMA  16.2    2.84     5.27      0.736
## 5 73        TOLI~ 043     73043  ANZOÁTEG~ 20.9    2.65     5.11      1.28 
## 6 73        TOLI~ 055     73055  ARMERO    14.9    2.11     6.21      0.362
## # ... with 3 more variables: HACINAMIENTO <dbl>, INACISTENCIA <dbl>,
## #   ECONOMIA <dbl>

Veamos cuál es el municipio con el mayor porcentaje de NBI:

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 73        TOLI~ 217     73217  COYAIMA    37.2    11.0     10.1      11.6
## # ... with 3 more variables: HACINAMIENTO <dbl>, INACISTENCIA <dbl>,
## #   ECONOMIA <dbl>

Veamos cuál es el municipio con el porcentaje más bajo de NBI:

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 73        TOLI~ 001     73001  IBAGUÉ     5.34   0.530     1.10     0.242
## # ... with 3 more variables: HACINAMIENTO <dbl>, INACISTENCIA <dbl>,
## #   ECONOMIA <dbl>

Con el siguiente chunck ordenaremos los municipios en orden descendente de NBI:

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

desc_nbi

## # A tibble: 47 x 12
##    COD_DEPTO DEPTO COD_MUN CODIGO MUNICIPIO   NBI MISERIA VIVIENDA SERVICIOS
##    <chr>     <chr> <chr>   <chr>  <chr>     <dbl>   <dbl>    <dbl>     <dbl>
##  1 73        TOLI~ 217     73217  COYAIMA    37.2   11.0     10.1     11.6  
##  2 73        TOLI~ 067     73067  ATACO      34.5   10.1     10.8     13.3  
##  3 73        TOLI~ 616     73616  RIOBLANCO  33.9   10.0     11.5     11.4  
##  4 73        TOLI~ 555     73555  PLANADAS   30.4    8.72    11.9     10.1  
##  5 73        TOLI~ 504     73504  ORTEGA     27.1    6.40     7.44     7.41 
##  6 73        TOLI~ 675     73675  SAN ANTO~  25.9    5.41    10.0      2.93 
##  7 73        TOLI~ 168     73168  CHAPARRAL  21.0    4.46     9.66     2.62 
##  8 73        TOLI~ 043     73043  ANZOÁTEG~  20.9    2.65     5.11     1.28 
##  9 73        TOLI~ 624     73624  ROVIRA     20.6    3.33     7.15     0.870
## 10 73        TOLI~ 483     73483  NATAGAIMA  20.5    4.74     5.58     3.10 
## # ... with 37 more rows, and 3 more variables: HACINAMIENTO <dbl>,
## #   INACISTENCIA <dbl>, ECONOMIA <dbl>

6.Unir datos de NBI a municipios

Anteriormente se tienen descargados los municipios de Tolima, los cuales se van a leer estos datos usando la librería sf:

munic <- st_read("C:/Users/Pablo/Documents/Agro/73_TOLIMA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

## Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\Pablo\Documents\Agro\73_TOLIMA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 47 features and 9 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.319342
## CRS:            4326

Miremos qué hay dentro de los atributos MPIO_CCDGO y MPIO_CNMBR:

head(munic$MPIO_CCDGO)

## [1] 73001 73024 73026 73030 73043 73055
## 47 Levels: 73001 73024 73026 73030 73043 73055 73067 73124 73148 ... 73873

head(munic$MPIO_CNMBR)

## [1] IBAGUÉ     ALPUJARRA  ALVARADO   AMBALEMA   ANZOÁTEGUI ARMERO    
## 47 Levels: ALPUJARRA ALVARADO AMBALEMA ANZOÁTEGUI ARMERO ATACO ... VILLARRICA

Para realizar la unión entre los municipios y los datos de NBI, utilizaremos la función left_join.

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

## Warning: Column `MPIO_CCDGO`/`CODIGO` joining factor and character vector,
## coercing into character vector

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:  POLYGON
## dimension:      XY
## bbox:           xmin: -75.5221 ymin: 3.276212 xmax: -74.7059 ymax: 5.147066
## CRS:            4326
##    MUNICIPIO MPIO_CCDGO       NBI                       geometry
## 1     IBAGUÉ      73001  5.344706 POLYGON ((-75.36517 4.67598...
## 2  ALPUJARRA      73024 12.008841 POLYGON ((-74.99192 3.50459...
## 3   ALVARADO      73026 13.915375 POLYGON ((-74.97668 4.70276...
## 4   AMBALEMA      73030 16.171061 POLYGON ((-74.74744 4.94928...
## 5 ANZOÁTEGUI      73043 20.879989 POLYGON ((-75.20098 4.71414...
## 6     ARMERO      73055 14.921420 POLYGON ((-74.88687 5.14696...

Para verificar que la unión fue correcta podremos ver el NBI de culquier municipio, por ejemplo el NBI de Armero es 14.921420 el cual podremos observar en las tablas anteriores.

Ahora, se vamos a reproyectar los municipios:

nbi_munic_new <- st_transform(nbi_munic, crs = 3116)

7.Ejemplos de mapas temáticos

Para realizar mapas temáticos se va a utilizar la librería cartography que tiene como objetivo obtener mapas temáticos con la calidad visual de los que se construyen en un software de mapeo o en un SIG. Esta librería usa objetos sf o sp para producir los gráficos base. Para objetos espaciales es preferible usar el formato sf.

7.1 Mapa base de OpenStreetMap y símbolos proporcionales

Las funciones getTiles() y tilesLayer() descargan y muestran los mosaicos de OpenStreetMap. La función propSymbolsLayer() muestra símbolos con áreas proporcionales a una variable cuantitativa (por ejemplo, NBI). Hay varios símbolos disponibles (círculos, cuadrados, barras). El argumento inches (pulgadas) se usa para personalizar los tamaños de los símbolos.

#install.packages("osmdata")
library(osmdata)

## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

library(rgdal)

## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/Pablo/Documents/R/R-3.6.3/library/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Pablo/Documents/R/R-3.6.3/library/rgdal/proj
##  Linking to sp version: 1.4-1

mun.osm <- getTiles(
x = nbi_munic_new, 
type = "OpenStreetMap", 
zoom = 8,
cachedir = TRUE,
crop = FALSE
)

#establecer márgenes
opar <- par(mar = c(0,0,1.2,0))

# trazar mosaicos osm
tilesLayer(x = mun.osm)

# trazar municipios (solo se trazan las fronteras)
plot(st_geometry(nbi_munic_new), col = NA, border = "grey", add=TRUE)

# trazar NBI
propSymbolsLayer(
  x = nbi_munic_new, 
  var = "NBI", 
  inches = 0.15, 
  col = "brown4",
  legend.pos = "topright",  
  legend.title.txt = "Total NBI"
)
# diseño
layoutLayer(title = " Distribución NBI en Tolima",
            sources = "Fuentes: DANE, 2018\n© OpenStreetMap",
            author = " Juan Pablo Ramírez ",
            frame = TRUE, north = FALSE, tabtitle = TRUE)

north(pos = "topleft")

En el mapa podemos observar que el circulo con mayor diametro corresponde con el municipio el cual posee mayor miseria (Coyaima) y el de menor miseria (Ibagué).

7.2 Mapas Coropléticos

En los mapas de coropletas, las áreas se sombreadas varian de acuerdo a una variable cuantitativa.Estas se usan para representar proporciones o índices. Con el comando choroLayer() muestra mapas CoropléticosLos argumentos nclass, method y breaks permiten personalizar la clasificación de la variable. El comando 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(). Tambien se puede consultar display.carto.all().

opar <- par(mar = c(0,0,1.2,0))
par(bg="grey90")
plot(st_geometry(nbi_munic_new), col = NA, border = NA, bg = "#aadaff")

# Grafica NBI
choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "harmo.pal", n1 = 5),
  border = "white", 
  lwd = 0.5,
  legend.pos = "right", 
  legend.title.txt = "NBI",
  add = TRUE
) 

layoutLayer(title = "Distribución de NBI en Tolima", 
            sources = "Fuente: DANE, 2018",
            author = " Juan Pablo Ramírez ", 
            frame = TRUE, north = TRUE, tabtitle = TRUE, col="black") 

north(pos = "topleft")

par(opar)

El mapa anterior se puede cambiar de acuerdo a gustos, por ejemplo el siguiente mapa es utilizando otras opciones de color, borde, entre otros.

opar <- par(mar = c(0,0,1.2,0))
par(bg="grey90")
plot(st_geometry(nbi_munic_new), col = NA, border = NA, bg = "#aadaff")

choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=8,
  col = carto.pal(pal1 = "sand.pal", n1 = 8),
  border = "black", 
  lwd = 1,
  legend.pos = "topleft", 
  legend.title.txt = "NBI",
  add = TRUE
) 

layoutLayer(title = "Distribución de NBI en Tolima", 
            sources = "Source: DANE, 2018",
            author = "Juan Pablo Ramírez", 
            frame = TRUE, north = TRUE, tabtitle = TRUE, col="black") 

north(pos = "topleft")

Estos mapas pueden ser más comodos para el usiario, ya que es más grafico y facil de diferenciar con la escala de colores, de este modo facilitando su interpretación y análisis de datos. En el mapa anterior se puede observar que el rango es de (5-37) el cual representa el NBI esta es una variable que mide la desigualdad económica de los distintos municipios, de este modo los municipios que posean un color marron más oscuro son los municipios con mayor desigualdad y los municipios con color marron menos oscuros son aquellos que poseen menor desigualdad.

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(). Lo primero que necesitamos es crear una variable cualitativa. Para ello usaremos la función mutate para esta tarea.

nbi_munic_2 <- dplyr::mutate(nbi_munic_new, pobreza = ifelse(MISERIA > 20, "Extrema", ifelse(HACINAMIENTO > 5, "Alta" , "Intermedia" )))

Tenga en cuenta que el nuevo atributo se llama pobreza. Su valor depende de los valores umbral definidos. Hay que asegurarse de verificar la sintaxis del comando ifelse para comprender qué está sucediendo.

head(nbi_munic_2)

## Simple feature collection with 6 features and 21 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 839613.7 ymin: 854083.5 xmax: 930308.6 ymax: 1060973
## CRS:            EPSG:3116
##   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                            MPIO_CRSLC
## 1         73      73001     IBAGUÉ                                  1906
## 2         73      73024  ALPUJARRA Decreto 650 del 13 de Octubre de 1887
## 3         73      73026   ALVARADO                                  1540
## 4         73      73030   AMBALEMA                                  1627
## 5         73      73043 ANZOÁTEGUI  Ordenanza 21 del 30 de Marzo de 1915
## 6         73      73055     ARMERO Decreto 1049 de Septiembre 29 de 1908
##   MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area COD_DEPTO  DEPTO
## 1  1377.2552      2017     TOLIMA   2.480389 0.11217138        73 TOLIMA
## 2   501.5113      2017     TOLIMA   1.054714 0.04080342        73 TOLIMA
## 3   344.3356      2017     TOLIMA   1.069589 0.02805438        73 TOLIMA
## 4   238.0448      2017     TOLIMA   1.009150 0.01940279        73 TOLIMA
## 5   469.7113      2017     TOLIMA   1.310177 0.03826692        73 TOLIMA
## 6   440.6308      2017     TOLIMA   1.444100 0.03592418        73 TOLIMA
##   COD_MUN  MUNICIPIO       NBI   MISERIA  VIVIENDA SERVICIOS HACINAMIENTO
## 1     001     IBAGUÉ  5.344706 0.5303010 1.0968326 0.2424812     1.577544
## 2     024  ALPUJARRA 12.008841 0.6876228 0.9332024 1.3998035     3.266208
## 3     026   ALVARADO 13.915375 2.0262217 2.8903456 1.0727056     3.530989
## 4     030   AMBALEMA 16.171061 2.8356836 5.2728387 0.7357449     4.659718
## 5     043 ANZOÁTEGUI 20.879989 2.6478725 5.1113160 1.2778290     5.164010
## 6     055     ARMERO 14.921420 2.1101889 6.2069574 0.3619989     3.593502
##   INACISTENCIA  ECONOMIA                       geometry    pobreza
## 1    1.2134178  1.809298 POLYGON ((857120.3 1008954,... Intermedia
## 2    1.0314342  6.164047 POLYGON ((898392.9 879339.7... Intermedia
## 3    0.9535161  7.568534 POLYGON ((900234.9 1011849,... Intermedia
## 4    1.3641937  7.648682 POLYGON ((925698.5 1039082,... Intermedia
## 5    2.0550652 10.420234 POLYGON ((875348.3 1013143,...       Alta
## 6    2.4280417  4.997351 POLYGON ((910260.5 1060962,... Intermedia

library(sf)
library(cartography)

opar <- par(mar = c(0,0,1.2,0))
plot(st_geometry(nbi_munic_2), col = "#f2efe9", border="#b38e43", bg = "#aad3df", lwd = 0.5)

propSymbolsTypoLayer(
  x= nbi_munic_2,
  var = "NBI",
  inches = 0.3,
  symbols = "square",
  border = "white",
  lwd = 0.5,
  legend.var.pos = "left",
  legend.var.title.txt = "NBI",
  var2 = "pobreza",
  legend.var2.values.order = c( "Alta", "Intermedia"),
  col = carto.pal(pal1 = "multi.pal", n1= 3),
  legend.var2.pos = "right",
  legend.var2.title.txt = "Pobreza"
)

layoutLayer(title = "Distribución de NBI en Tolima",
            author = "Juan Pablo Ramírez",
            sources = "Fuente: DANE, 2018",
            scale = 1, tabtitle = TRUE, frame = TRUE)

north(pos = "topleft")

Como en el Departamento de Tolima no hay Miseria mayor a 20 no se observan lugares de pobreza Extrema, pero si cambiamos esta variable a Miseria mayor a 11 si se observa que solo 1 municipio que posee pobreza extrema.

nbi_munic_2.1 <- dplyr::mutate(nbi_munic_new, pobreza = ifelse(MISERIA > 11, "Extrema", ifelse(HACINAMIENTO > 5, "Alta" , "Intermedia" )))

head(nbi_munic_2.1)

## Simple feature collection with 6 features and 21 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 839613.7 ymin: 854083.5 xmax: 930308.6 ymax: 1060973
## CRS:            EPSG:3116
##   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                            MPIO_CRSLC
## 1         73      73001     IBAGUÉ                                  1906
## 2         73      73024  ALPUJARRA Decreto 650 del 13 de Octubre de 1887
## 3         73      73026   ALVARADO                                  1540
## 4         73      73030   AMBALEMA                                  1627
## 5         73      73043 ANZOÁTEGUI  Ordenanza 21 del 30 de Marzo de 1915
## 6         73      73055     ARMERO Decreto 1049 de Septiembre 29 de 1908
##   MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area COD_DEPTO  DEPTO
## 1  1377.2552      2017     TOLIMA   2.480389 0.11217138        73 TOLIMA
## 2   501.5113      2017     TOLIMA   1.054714 0.04080342        73 TOLIMA
## 3   344.3356      2017     TOLIMA   1.069589 0.02805438        73 TOLIMA
## 4   238.0448      2017     TOLIMA   1.009150 0.01940279        73 TOLIMA
## 5   469.7113      2017     TOLIMA   1.310177 0.03826692        73 TOLIMA
## 6   440.6308      2017     TOLIMA   1.444100 0.03592418        73 TOLIMA
##   COD_MUN  MUNICIPIO       NBI   MISERIA  VIVIENDA SERVICIOS HACINAMIENTO
## 1     001     IBAGUÉ  5.344706 0.5303010 1.0968326 0.2424812     1.577544
## 2     024  ALPUJARRA 12.008841 0.6876228 0.9332024 1.3998035     3.266208
## 3     026   ALVARADO 13.915375 2.0262217 2.8903456 1.0727056     3.530989
## 4     030   AMBALEMA 16.171061 2.8356836 5.2728387 0.7357449     4.659718
## 5     043 ANZOÁTEGUI 20.879989 2.6478725 5.1113160 1.2778290     5.164010
## 6     055     ARMERO 14.921420 2.1101889 6.2069574 0.3619989     3.593502
##   INACISTENCIA  ECONOMIA                       geometry    pobreza
## 1    1.2134178  1.809298 POLYGON ((857120.3 1008954,... Intermedia
## 2    1.0314342  6.164047 POLYGON ((898392.9 879339.7... Intermedia
## 3    0.9535161  7.568534 POLYGON ((900234.9 1011849,... Intermedia
## 4    1.3641937  7.648682 POLYGON ((925698.5 1039082,... Intermedia
## 5    2.0550652 10.420234 POLYGON ((875348.3 1013143,...       Alta
## 6    2.4280417  4.997351 POLYGON ((910260.5 1060962,... Intermedia

library(sf)
library(cartography)

opar <- par(mar = c(0,0,1.2,0))
plot(st_geometry(nbi_munic_2.1), col = "#f2efe9", border="#b38e43", bg = "#aad3df", lwd = 0.5)

propSymbolsTypoLayer(
  x= nbi_munic_2.1,
  var = "NBI",
  inches = 0.3,
  symbols = "square",
  border = "white",
  lwd = 0.5,
  legend.var.pos = "left",
  legend.var.title.txt = "NBI",
  var2 = "pobreza",
  legend.var2.values.order = c( "Extrema" , "Alta", "Intermedia"),
  col = carto.pal(pal1 = "multi.pal", n1= 3),
  legend.var2.pos = "right",
  legend.var2.title.txt = "Pobreza"
)

layoutLayer(title = "Distribución de NBI en Tolima",
            author = "Juan Pablo Ramírez",
            sources = "Fuente: DANE, 2018",
            scale = 1, tabtitle = TRUE, frame = TRUE)

north(pos = "topleft")

7.4 Mapas de etiquetas

En este mapa se van a combinar las funciones choroLayer() y labelLayer().

library(sf)
library(cartography)

opar <- par(mar = c(0,0,1.2,0))
par(bg="grey25")
plot(st_geometry(nbi_munic_2), col = "#e4e9de", border = "darkseagreen4", 
     bg = "grey75", lwd = 0.5)

choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "sand.pal", n1 = 5),
  border = "white", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "NBI",
  add = TRUE
) 

labelLayer(
  x = nbi_munic_2, 
  txt = "MUNICIPIO", 
  col= "white", 
  cex = 0.4, 
  font = 4,
  halo = TRUE, 
  bg = "grey25", 
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)

layoutLayer(
  title = "Municipios de Tolima", 
  sources = "Fuente: DANE, 2018",  
  author = "Juan Pablo Ramírez", 
  frame = TRUE,
  north = TRUE, 
  tabtitle = TRUE, 
  theme = "taupe.pal"
) 

Este mapa es muy útli y practico, ya que se usan escalas de colores para indicar rangos de valores de NBI junto con los nombres de los municipios. Esto nos ayuda a identificar que municipio posee mayor o menor NBI de manera gráfica, algo que en los anteriores mapas no se podria identificar facilmente si no se conoce la ubicación del municipio de interes.

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. Se usa la función smoothLayer() el cual depende en gran medida del paquete SpatialPosition. Esta función 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.

Para realizar este mapa, se utilizaremos otro conjunto de datos. En este caso, se tendrán en cuenta los datos sobre la producción de café en 2018, para el departamento de Tolima ya que es el cultivo con mayor área sembrada en este departamento.

crops2018  <-  read_excel("C:/Users/Pablo/Documents/Agro/EVA_Tolima_3.xlsx")

head(crops2018) 

## # A tibble: 6 x 16
##   COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO SUBGRUPO CULTIVO  YEAR PERIODO
##     <dbl> <chr>          <dbl> <chr>     <chr> <chr>    <chr>   <dbl> <chr>  
## 1      73 TOLIMA         73055 ARMERO (~ HORT~ ACELGA   ACELGA   2012 2012A  
## 2      73 TOLIMA         73024 ALPUJARRA HORT~ ACELGA   ACELGA   2015 2015A  
## 3      73 TOLIMA         73124 CAJAMARCA TUBE~ ACHIRA   ACHIRA   2012 2012   
## 4      73 TOLIMA         73283 FRESNO    FRUT~ AGUACATE AGUACA~  2007 2007   
## 5      73 TOLIMA         73443 MARIQUITA FRUT~ AGUACATE AGUACA~  2007 2007   
## 6      73 TOLIMA         73026 ALVARADO  FRUT~ AGUACATE AGUACA~  2007 2007   
## # ... with 7 more variables: AREA_SEMBRADA <dbl>, AREA_COSECHADA <dbl>,
## #   PRODUCCION <dbl>, RENDIMIENTO <dbl>, ESTADO_FISICO_PRODUCCION <chr>,
## #   NOMBRE_CIENTIFICO <chr>, CICLO_DEL_CULTIVO <chr>

Ahora, filtraremos las filas que representan solo datos de cafe:

crops2018 %>%
  filter(CULTIVO == "CAFE") -> cafe2018

head(cafe2018)

## # A tibble: 6 x 16
##   COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO SUBGRUPO CULTIVO  YEAR PERIODO
##     <dbl> <chr>          <dbl> <chr>     <chr> <chr>    <chr>   <dbl> <chr>  
## 1      73 TOLIMA         73411 LIBANO    OTRO~ CAFE     CAFE     2007 2007   
## 2      73 TOLIMA         73555 PLANADAS  OTRO~ CAFE     CAFE     2007 2007   
## 3      73 TOLIMA         73001 IBAGUE    OTRO~ CAFE     CAFE     2007 2007   
## 4      73 TOLIMA         73283 FRESNO    OTRO~ CAFE     CAFE     2007 2007   
## 5      73 TOLIMA         73168 CHAPARRAL OTRO~ CAFE     CAFE     2007 2007   
## 6      73 TOLIMA         73067 ATACO     OTRO~ CAFE     CAFE     2007 2007   
## # ... with 7 more variables: AREA_SEMBRADA <dbl>, AREA_COSECHADA <dbl>,
## #   PRODUCCION <dbl>, RENDIMIENTO <dbl>, ESTADO_FISICO_PRODUCCION <chr>,
## #   NOMBRE_CIENTIFICO <chr>, CICLO_DEL_CULTIVO <chr>

A continuación crearemos un nuevo atributo que coincida con los códigos de los municipios:

cafe2018$TEMP <-  as.character(cafe2018$COD_MUN)

cafe2018$MPIO_CCDGO <- as.factor(cafe2018$TEMP)

Ahora, vamos a realizar la unión:

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 26 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -75.5221 ymin: 4.256457 xmax: -74.96571 ymax: 4.700691
## CRS:            4326
##   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR
## 1         73      73001     IBAGUÉ       1906   1377.255      2017     TOLIMA
## 2         73      73001     IBAGUÉ       1906   1377.255      2017     TOLIMA
## 3         73      73001     IBAGUÉ       1906   1377.255      2017     TOLIMA
## 4         73      73001     IBAGUÉ       1906   1377.255      2017     TOLIMA
## 5         73      73001     IBAGUÉ       1906   1377.255      2017     TOLIMA
## 6         73      73001     IBAGUÉ       1906   1377.255      2017     TOLIMA
##   Shape_Leng Shape_Area COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO
## 1   2.480389  0.1121714      73       TOLIMA   73001    IBAGUE
## 2   2.480389  0.1121714      73       TOLIMA   73001    IBAGUE
## 3   2.480389  0.1121714      73       TOLIMA   73001    IBAGUE
## 4   2.480389  0.1121714      73       TOLIMA   73001    IBAGUE
## 5   2.480389  0.1121714      73       TOLIMA   73001    IBAGUE
## 6   2.480389  0.1121714      73       TOLIMA   73001    IBAGUE
##               GRUPO SUBGRUPO CULTIVO YEAR PERIODO AREA_SEMBRADA AREA_COSECHADA
## 1 OTROS PERMANENTES     CAFE    CAFE 2007    2007          7945           6840
## 2 OTROS PERMANENTES     CAFE    CAFE 2008    2008          8278           7036
## 3 OTROS PERMANENTES     CAFE    CAFE 2009    2009          8313           7323
## 4 OTROS PERMANENTES     CAFE    CAFE 2010    2010          8234           7323
## 5 OTROS PERMANENTES     CAFE    CAFE 2011    2011          8798           7750
## 6 OTROS PERMANENTES     CAFE    CAFE 2012    2012          8433           7000
##   PRODUCCION RENDIMIENTO ESTADO_FISICO_PRODUCCION NOMBRE_CIENTIFICO
## 1       6156        0.90   CAFE VERDE EQUIVALENTE    COFFEA ARABICA
## 2       7921        1.13   CAFE VERDE EQUIVALENTE    COFFEA ARABICA
## 3       7323        1.00   CAFE VERDE EQUIVALENTE    COFFEA ARABICA
## 4       9300        1.27   CAFE VERDE EQUIVALENTE    COFFEA ARABICA
## 5       4257        0.55   CAFE VERDE EQUIVALENTE    COFFEA ARABICA
## 6       7000        1.00   CAFE VERDE EQUIVALENTE    COFFEA ARABICA
##   CICLO_DEL_CULTIVO  TEMP                       geometry
## 1        PERMANENTE 73001 POLYGON ((-75.36517 4.67598...
## 2        PERMANENTE 73001 POLYGON ((-75.36517 4.67598...
## 3        PERMANENTE 73001 POLYGON ((-75.36517 4.67598...
## 4        PERMANENTE 73001 POLYGON ((-75.36517 4.67598...
## 5        PERMANENTE 73001 POLYGON ((-75.36517 4.67598...
## 6        PERMANENTE 73001 POLYGON ((-75.36517 4.67598...

Reproyección de los municipios:

rep_cafe <- st_transform(cafe_munic, crs = 3116)

Ahora realizaremos el mapa:

opar <- par(mar = c(0,0,1.2,0))
plot(st_geometry(rep_cafe), col = NA, border = "black", bg = "grey75")

smoothLayer(
  x = rep_cafe, 
  var = 'PRODUCCION',
  typefct = "exponential",
  span = 20000,
  beta = 2,
  nclass = 10,
  col = carto.pal(pal1 = 'harmo.pal', n1 = 10),
  border = "grey",
  lwd = 0.1, 
  mask = rep_cafe, 
  legend.values.rnd = -3,
  legend.title.txt = "Producción",
  legend.pos = "right", 
  add=TRUE)

text(x = 650000, y = 1200000, cex = 0.6, adj = 0, font = 3, labels = "Función distancia:\n- type = exponential\n- beta = 2\n- span = 15 km")

layoutLayer(title = "Distribución de la Producción de Café en Tolima",
            sources = "Fentes: DANE y MADR, 2018",
            author = "Juan Pablo Ramírez",
            frame = FALSE, north = FALSE, tabtitle = TRUE, theme = "green.pal")

north(pos = "topleft")

8.Guardar mapas

Ahora se va a realizar otro mapa de producción de Café en 2018. Esta vez se van a usar símbolos proporcionales y mapas Coropléticos La salida se guardará como un archivo .png.Primero, algunas observaciones:

• propSymbolsChoroLaer() crea una 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(). El siguiente chunk no muestra un mapa. En cambio, escribe el mapa en el nombre de archivo cafe_2018.png debajo del directorio de trabajo.

png("C:/Users/Pablo/Documents/PABLO/Geomatica/cafe_2018.png", width = 2048, height = 1526)
opar <- par(mar = c(0,0,5,5))

plot(st_geometry(rep_cafe), col="darkseagreen3", border="darkseagreen4",  
     bg = "white", lwd = 0.6)

propSymbolsChoroLayer(x = rep_cafe, var = "PRODUCCION", var2 = "RENDIMIENTO",
                      col = carto.pal(pal1 = "harmo.pal", n1 = 3,
                                      pal2 = "red.pal", n2 = 3),
                      inches = 0.8, method = "q6",
                      border = "grey50", lwd = 1,
                      legend.title.cex = 1.5,
                      legend.values.cex = 1.0,
                      legend.var.pos = "right", 
                      legend.var2.pos = "left",
                      legend.var2.values.rnd = 2,
                      legend.var2.title.txt = "Rendimiento\n(en Ton/Ha)",
                      legend.var.title.txt = "Producción de Café en 2018",
                      legend.var.style = "e")

labelLayer(
  x = rep_cafe, 
  txt = "MPIO_CNMBR", 
  col= "white", 
  cex = 1.0, 
  font = 4,
  halo = FALSE, 
  bg = "white", 
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)

layoutLayer(title="Producción y Rendimiento de Café en Tolima, 2018",
            author = "Juan Pablo Ramírez", 
            sources = "Fuentes: MADR & DANE, 2018", 
            scale = 50, tabtitle = FALSE, frame = TRUE)

north(pos = "topleft")

title(main="Producción y Rendimiento de Café en Tolima, 2018", cex.main=3,
      sub= "Fuentes: MADR & DANE, 2018", cex.sub=2)

graticule = TRUE

par(opar)

dev.off()

## png 
##   2

Ahora, que ya hemos guardado el mapa como una imagen, se puede agregar al cuaderno de R Markdown usando la sintaxis de Markdown de la siguiente manera:

Producción de Café en Tolima, 2018.

sessionInfo()

## R version 3.6.3 (2020-02-29)
## Platform: i386-w64-mingw32/i386 (32-bit)
## Running under: Windows 10 x64 (build 18362)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252   
## [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Spain.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] rgdal_1.4-8           osmdata_0.1.3         SpatialPosition_2.0.1
##  [4] cartography_2.4.1     sf_0.9-2              raster_3.0-12        
##  [7] rgeos_0.5-2           sp_1.4-1              readxl_1.3.1         
## [10] forcats_0.5.0         stringr_1.4.0         dplyr_0.8.4          
## [13] purrr_0.3.3           readr_1.3.1           tidyr_1.0.2          
## [16] tibble_3.0.1          ggplot2_3.3.0         tidyverse_1.3.0      
## 
## loaded via a namespace (and not attached):
##  [1] jsonlite_1.6.1     rstudioapi_0.11    generics_0.0.2     magrittr_1.5      
##  [5] isoband_0.2.0      gtable_0.3.0       rmarkdown_2.1      vctrs_0.2.4       
##  [9] fs_1.4.1           class_7.3-15       hms_0.5.3          utf8_1.1.4        
## [13] xml2_1.2.2         pillar_1.4.3       htmltools_0.4.0    curl_4.3          
## [17] haven_2.2.0        broom_0.5.6        cellranger_1.1.0   lattice_0.20-40   
## [21] slippymath_0.3.1   KernSmooth_2.23-16 tidyselect_1.0.0   lubridate_1.7.8   
## [25] knitr_1.28         lifecycle_0.2.0    pkgconfig_2.0.3    R6_2.4.1          
## [29] digest_0.6.25      xfun_0.12          colorspace_1.4-1   stringi_1.4.6     
## [33] yaml_2.2.1         codetools_0.2-16   evaluate_0.14      lwgeom_0.2-3      
## [37] fansi_0.4.1        httr_1.4.1         compiler_3.6.3     cli_2.0.2         
## [41] withr_2.1.2        backports_1.1.5    munsell_0.5.0      DBI_1.1.0         
## [45] modelr_0.1.6       Rcpp_1.0.3         dbplyr_1.4.3       png_0.1-7         
## [49] ellipsis_0.3.0     assertthat_0.2.1   classInt_0.4-3     units_0.6-6       
## [53] tools_3.6.3        reprex_0.3.0       e1071_1.7-3        scales_1.1.0      
## [57] crayon_1.3.4       glue_1.3.1         rlang_0.4.5        nlme_3.1-144      
## [61] rvest_0.3.5        grid_3.6.3