1. Introducción

La cartografía temática actualmente ha sido más explotada debido a su potencial en diferentes áreas del conocimiento, implementada desde científicos, diseñadores gráficos, y hasta por expertos en sistemas de información geográfica, debido a que ayuda a comprender la organización y el funcionamiento en relaciones socioterritoriales (Fenómenos naturales, sociales, ecoómicos, culturales, históricos o políticos). Así, es posible analizar aspectos generales de un modelo, junto con particularidades de una zona específica (Régnauld, H. y Lefort, I., 2015).

En los mapas tempaticos se presenta de forma convencional, por medio de signos y símbolos, información tanto cualitativa como cuantitativa, correspondiente a la iformación georgráfica de interés, aquellas relaciones socioculturales (Gómez, 2004).

Por eso, en este cuaderno de R se ilustra la cartografía temática para el departamento de Casanare - Colombia, para poder comprender diferentes dinámicas sociales que se dan en el territorio. Para esto, se han planteado diversos mapas, según sea el análisis que se quiera realizar o la forma en la que se quiera representar una información:

Mapa coroplético: Estos mapas contienen áreas que están sombreadas o con patrones en proporción a la variable estadística que se muestra en el mapa. Los datos se agregan sobre unidades de área predefinidas (área definida políticamente o administrativa - censo o zip).

Mapa isoplético: Un mapa de curvas de nivel que muestra información fluida y continua (como datos meteorológicos o de contaminación); los datos se representan mediante líneas que conectan puntos de igual valor numérico.

Mapa dasimétrico: Un mapa temático que utiliza símbolos para clasificar espacialmente datos volumétricos; sirve como una alternativa a los mapas coropléticos.

2. Datos

Se utilizarán datos de Necesidades Básicas Insatisfechas (NBI), provenientes del censo del DANE (2018)

En este caso, los mapas temáticos son sumamente útiles para el análisis y la representación de información demofráfica.

Una vez descargados los datos NBI, se remueven los datos de aquellos municipios que no se encuentran dentro del departamento de interés, en este caso, Casanare.

3. Preparación

Inicialmente, es necesario realizar una limpieza al ambiente de trabajo de R, mediante la función “rm(list=ls())”, y seleccionar la carpeta en el ordenador en la cual se almacenarán todos los archivos del proyecto, mediante la función “setwd”.

###
rm(list=ls())
setwd("C:/Users/yarya/Desktop/Geomatica/ProyectoCasanare")
getwd()
## [1] "C:/Users/yarya/Desktop/Geomatica/ProyectoCasanare"
##

Una vez realizado esto, es necesraio instalar las librerías que se necesitan para ejecutar las funciones necesarias. La sigueinte función permite ejecutar el código, sin que se instalen siempre las mismas librerías, por lo que se pueden escribir aquellas que son necesarias para el trabajo práctico, y el programa instalará únicamente las que hacen falta.

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

Ahora si es posible cargar las librerías necesarias.

###
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(readxl)
library(rgeos)
## Loading required package: sp
## rgeos version: 0.5-5, (SVN revision 640)
##  GEOS runtime version: 3.8.0-CAPI-1.13.1 
##  Linking to sp version: 1.4-2 
##  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.8.0, GDAL 3.0.4, PROJ 6.3.1
library(cartography)
## Warning: package 'cartography' was built under R version 4.0.3
library(SpatialPosition)
## Warning: package 'SpatialPosition' was built under R version 4.0.3
##

4. Leer datos NBI - Casanare

Inicialmente, es necesario cargar el archivo descargado previamente, provenientes del DANE, y observar los atributos de estos datos.

###
nbi <- read_excel("C:/Users/yarya/Desktop/Geomatica/ProyectoCasanare/NBI_Casanare.xlsx")
nbi$CODIGO<-as.character(nbi$CODIGO)
nbi
##

Vamos a encontrar el municipio con mayor porcentaje de NBI, mediante el código “slice(which.max(NBI))”

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

max_nbi
##

Esto puede relacionarse en gran medida a que en el municipio de Tamará, a diferencia de otros municipios aledaños, la disponibilidad de computadores, redes de internet, y básicamente una infrasetrura la cual aparenta no haber evolucionado en los últimos 25 años (Botero, 2015)

Y ahora el municipio con el menor porcentaje de NBI, utilizando el código “which.min”

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

min_nbi
##

Por otro lado, Monterrey presenta el NB1 más bajo de la región (7,97%), lo cual puede relacionarse con la incorporación en los útimos años de la ganadería intensiva y el doble proósito, junto con las actividades de turismo las cuales son muy deseadas por los visitantes, permitiendo así un sustento económico favorable. Sin embargo, el sustento de muchos campesinos con agricultura tradicional, más la producción de arroz mecanizado en gran escala, permite el sustento de pequeños, medianos, y grandes agricultores.

Finalmente, vamos a organizar los datos en orden descendente, de mayor a menor porcentraje de NBI, para así lograr analizar los datos de forma más eficiente. Esto puede realizarse mediante el código “arrange(desc(NBI))”

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

desc_nbi
##

5. Unir datos NBI a municipios

Para relacionar los datos de NBI con los difernetes municipios, vamos a utilizar los datos políticos proporcionados por el DANE.

###
munic <- st_read("C:/Users/yarya/Desktop/Geomatica/Clase5/DatosCasanare2017/85_CASANARE/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\yarya\Desktop\Geomatica\Clase5\DatosCasanare2017\85_CASANARE\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 19 features and 9 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -73.07777 ymin: 4.287476 xmax: -69.83591 ymax: 6.346111
## geographic CRS: WGS 84
##

Ahora, es necesario verificar los datos presnetes con relación a los municipios, dentro del atributo “MPIO_CNMBR”

De igual forma, es posible verificar la información correspondiente al códogo del municipio, para asegurar la adecuada unión posterior

munic$MPIO_CNMBR
##  [1] "YOPAL"                "AGUAZUL"              "CHAMEZA"             
##  [4] "LA SALINA"            "MANÍ"                 "MONTERREY"           
##  [7] "NUNCHIA"              "OROCUÉ"               "PORE"                
## [10] "RECETOR"              "SABANALARGA"          "SACAMA"              
## [13] "SAN LUIS DE PALENQUE" "TAMARA"               "TAURAMENA"           
## [16] "TRINIDAD"             "VILLANUEVA"           "HATO COROZAL"        
## [19] "PAZ DE ARIPORO"
munic$MPIO_CCDGO
##  [1] "85001" "85010" "85015" "85136" "85139" "85162" "85225" "85230" "85263"
## [10] "85279" "85300" "85315" "85325" "85400" "85410" "85430" "85440" "85125"
## [19] "85250"
##

Se observa que los cógigos si corresponden a los presnetes en el archivo de NBI proporcionado por la DIAN, por lo que se puede realizar la unión de ambos archivos con base al códogo del municipio.

###
nbi_munic = left_join(munic, nbi, by=c("MPIO_CCDGO"="CODIGO"))
nbi_munic %>%
  dplyr::select(MUNICIPIO, MPIO_CCDGO, NBI)  ->  check_nbi_munic

check_nbi_munic
##

Se observa dos valores de N/A relacionados con los códogos 85015 y 85136, correspondientes a La Salina y Chameza, ya que estos municipios no se encuentran dentro de los datos proporcionados por el DANE.

Ahora, es necesario reproyectar los datos de municipios, ya que se encuentran en el sistema de referencia WGS 84.

###
nbi_munic_new <- st_transform(nbi_munic, crs = 3116)
nbi_munic_new
##

Notese que ahora los datos se encuentran ubicados en el sistema de referencia MAGNA-SIRGAS, correspondiente a 3116.

6. Mapas temáticos

Ahora, con la información cargada al sistema, es posible generar diferentes mápas temáticos. Estas funciones requeren el uso obligatorio de las librerias “sf” y “sp”, para producir las gráficos base.

6.1. OpenStreetMap Basemap y Símbolos Proporcionales

Este mapa permite relacionar el tamaño de los datos, es decir, en una escala, aquellos más grandes y más pequeños, con una forma geométrica de tamaño directamente proporcional. Así, se observa que las zonas con valores de NBI más bajos presnetan círculos anaranjados de tamaño menor, en comparación con auqellas zonas de NBI más alto (obsérvese el mapa 4 para relacionar los municipios).

###
mun.osm <- getTiles(
x = nbi_munic_new, 
type = "OpenStreetMap", 
zoom = 8,
cachedir = TRUE,
crop = FALSE
)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
opar <- par(mar = c(0,0,1.2,0))
tilesLayer(x = mun.osm)
plot(st_geometry(nbi_munic_new), col = NA, border = "grey", add=TRUE)
propSymbolsLayer(
  x = nbi_munic_new, 
  var = "NBI", 
  inches = 0.15, 
  col = "Orange",
  legend.pos = "topright",  
  legend.title.txt = "NBI Total",
  legend.frame = T,
  legend.style = "e"
  
)
layoutLayer(title = " Distribución del NBI en Casanare",
            sources = "Fuentes: DANE, 2018\n© OpenStreetMap",
            author = " David Agudelo ",
            frame = TRUE, north = FALSE, tabtitle = TRUE)
north(pos = "topleft")

##

6.2. Mapa coroplético

En este caso, la represnetación gráfica de la variable en cuestión se realiza mediante el uso de colores con diferentes tonos, relacionados con una escala. En este caso, es fundamentakl el uso del comando “choroLayer”.

En el presente caso, también se representarpan los valores de NBI para cada departamento, donde, en lo personal, la representación es mucho más sencilla de comprender y estéticamente más pulida. Se observan los municipios con tonalidades más oscuras aquellos con un NBI mayor, indicando una mayor cantidad de necesidades básicas insatisfechas.

###
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=5,
  col = carto.pal(pal1 = "sand.pal", n1 = 5),
  border = "white", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "NBI",
  add = TRUE
) 
layoutLayer(title = "Distribución del NBI en Casanare", 
            sources = "Fuente: DANE, 2018",
            author = "David Agudelo", 
            frame = TRUE, north = TRUE, tabtitle = TRUE, col="black") 
north(pos = "topleft")

##

En relación al cuaderno realizado anteriormente, se observa que aquellos municipios con mayor producción de cultivos de palma y arroz mecanizado, presentan un menor valor de NBI.

6.3. Simbolos proporcionales y mapa tipológico

Ahora, vamos a crear mapas de símbolos que son proporcionales a los valores de la primera variable, y el color refleja la modalidad de la segunda variable cualitativa.

Por lo anterior, lo primero que se debe hacer es crear una variable cualitativa y se usa la función “mutate” para realizar la acción.

Así, se puede visualizar el nuevo atributo llamado “Pobreza”, y su valor depende de los valores del umbral definidos por el usuario.

###
nbi_munic_2 <- dplyr::mutate(nbi_munic_new, Pobreza = ifelse(MISERIA > 20,"Extreme", ifelse(HACINAMIENTO > 5, "High", "Intermediate")))

nbi_munic_2
##

Como ningún municipio presneta condición de “Miseria”, la distribución del NBI en Casanare no presentará la opción de “Extrema”

Así, en el mapa los valores de color rojo representan pobreza alta, y los cuadros verdes represnetan pobreza intermedia.

###
opar <- par(mar = c(0,0,1.2,0))
plot(st_geometry(nbi_munic_2), col= "#f2efe9", border="#b38e43", bg = "#aadaff", 
     lwd = 0.5)
propSymbolsTypoLayer(
  x = nbi_munic_2, 
  var = "NBI", 
  inches = 0.3,
  symbols = "square",
  border = "white",
  lwd = .5,
  legend.var.pos = c(1400000, 1040000), 
  legend.var.title.txt = "NBI",
  var2 = "Pobreza",
  legend.var2.values.order = c("High","Intermediate"),
  col = carto.pal(pal1 = "multi.pal", n1 = 3),
  legend.var2.pos = c(1400000, 980000), 
  legend.var2.title.txt = "Pobreza"
) 
layoutLayer(title="Distribución del NBI en Casanare", 
            author = "David Agudelo", 
            sources = "Fuente: DANE, 2018", 
            scale = 1, tabtitle = TRUE, frame = TRUE)
north(pos = "topleft")

##

Se observa que la región más al noroccidente presentar valores de pobreza alta, mientras que las zonas hacia el suroriente del departamento presentan pobreza intermedia.

6.4. Mapa de etiquetas

Con este tipo de mapas, se pueden combinar las funciones “choroLayer” y “labelLayer”, esto permite ubicar los nombres de cada municipio sobre la ubicación específica.

###
opar <- par(mar = c(0,0,1.2,0))
par(bg="grey25")
plot(st_geometry(nbi_munic_2), col = NA, border = NA, 
     bg = "grey75", lwd = 0.5)
choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "green.pal", n1 = 5),
  border = "white", 
  lwd = 0.5,
  legend.pos =  c(1400000, 980000),
  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 = "NBI - Municipios de Casanare", 
  sources = "Fuente: DANE, 2018",  
  author = "David Agudelo", 
  frame = TRUE,
  north = TRUE, 
  tabtitle = TRUE, 
  theme = "taupe.pal"
) 

##

6.5. Mapa de Isopletas

Estos mapas se basan en el supuesto de que el fenomeno a representar es de distribución continua, y utilizan un enfoque de modelado de interceptación espacial, los cuales tinenen el objetivo de calcular indicadores basados en valores de stock ponderados por distancia.ermite una representación espacial del fenómeno independiente de la heterogeneidad inicial de la división territorial.

Para esto, se debe inicialmente cargar la información relacionada con la producción de los cultivos, en etse caso me enfocaré en el cultivo de palma de aceite.

###
Cultivos2018 <-read.csv("C:/Users/yarya/Desktop/Geomatica/ProyectoCasanare/Clase10/Datos_Siembra_Casanare.csv", sep = ";")
Cultivos2018$COD_N<-as.double(Cultivos2018$COD_N)
Cultivos2018$YEAR<-as.double(Cultivos2018$YEAR)
head(Cultivos2018)
##

Una vez subidos los datos, se reliza el filtro correspondiente, en este caso por cultivo, pero también se realiza el filtro por año para evitar problemas en el desarrollo del código.

Posterior a esto, se debe realizar la unión de ambas fuentes de datos, asegurándose de que un par de atributos sean idénticos. Esto es posible ya sea señalando que los atributos “COD_N” y “MPIO_CCDGO” son los mismos, o nombrandolos de la misma forma.

###
Cultivos2018 %>%
  filter(CULTIVO == "PALMA DE ACEITE", YEAR==2018)-> Palma2018
Palma2018$TEMP <-  as.character(Palma2018$COD_N)
Palma2018$MPIO_CCDGO <- as.factor(paste(Palma2018$TEMP))
Palma_munic = left_join(munic, Palma2018, by="MPIO_CCDGO")

Palma_munic
##
###
opar <- par(mar = c(0,0,1.2,0))
rep_Palma<- st_transform(Palma_munic, crs = 3116)
plot(st_geometry(rep_Palma),col = NA, border = "black", bg = "grey75")
smoothLayer(
  x = rep_Palma,
  var = 'PRODUCCION_TON',
  typefct = "exponential",
  span = 25000,
  beta = 2,
  nclass = 10,
  col = carto.pal(pal1 = 'green.pal', n1 = 10),
  border = "grey",
  lwd = 0.1, 
  mask = rep_Palma,
  legend.title.txt = "Producción",
  legend.pos = c(1400000, 980000), 
  add=TRUE
)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
text(x = 650000, y = 1200000, cex = 0.6, adj = 0, font = 3,  labels = 
       "Distance function:\n- type = exponential\n- beta = 2\n- span = 20 km")
layoutLayer(title = "Distribución de producción de Palma en Casanare",
            sources = "Fuente: DANE and MADR, 2018",
            author = "David Agudelo",
            frame = FALSE, north = FALSE, tabtitle = TRUE, theme = "green.pal")
north(pos = "topright")

##

7. Guardar mapas

Por último, vamos a crear un mapa relacionado con la producción de palma, para evidenciar el proceso de guardado como objeto png.

###
png("C:/Users/yarya/Desktop/Geomatica/ProyectoCasanare/Palma_2018.png", width = 2048, height = 1526)
opar <- par(mar = c(0,0,5,5))
plot(st_geometry(rep_Palma), col="darkseagreen3",
     bg = "white", lwd = 0.6)
propSymbolsChoroLayer(x = rep_Palma, var = "PRODUCCION_TON", var2 = "RENDIMIENTO_TON.HA",
                      col = carto.pal(pal1 = "green.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(in Ton/Ha)",
                      legend.var.title.txt = "Coffe Production in 2018",
                      legend.var.style = "e")
labelLayer(
  x = rep_Palma,
  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 de palma y rendimiento en Casanare, 2018",
            author = "David Agudelo", 
            sources = "Fuente: MADR & DANE, 2018", 
            scale = 50, tabtitle = FALSE, frame = TRUE)
north(pos = "topleft")
#
title(main="Coffee Production & Yield in Antioquia, 2018", cex.main=3,
      sub= "Source: MADR & DANE, 2018", cex.sub=2)
#
graticule = TRUE
#
par(opar)
### close the plot
dev.off()
## png 
##   2
knitr::include_graphics("C:/Users/yarya/Desktop/Geomatica/ProyectoCasanare/Palma_2018.png")

##

8. Bibliografía

Botero, J. 2015. Tamará, su historia reciente (Colombia).

DANE - Censo Nacional de Población y Vivienda (CNPV) 2018.

El Tiempo. Monterrey, polo de desarrollo. Disponible en: https://www.eltiempo.com/archivo/documento/MAM-166319

Gómez, M. 2004. Metodologías y técnicas de la cartografía temática. Universidad Nacional Autónoma de México.

Régnauld, H. y Lefort, I. (2015): L’image et la géographie : la progressive elaboration d’un nouveau régime épistémique, L’Information Géographique, 2015/4 (79): 8-12.