CARTOGRAFÍA TEMÁTICA

By Valentina Valle Velasco

19 de Mayo, 2020

Este es un cuaderno de R Markdown que nos muestra carotgrafía temática en el departamento del Magdalena, Colombia.
Tiene el objetivo de exponer, ante los estudiantes de la asignatura Geomática Básica de la Universidad Nacional de Colombia, cómo hacer mapas temáticos en R.

Cartografía Temática

Amanda Briney indica que un mapa temático se enfoca en mostrar temas de interés, como la distribución promedio de lluvia en un área o la densidad de población en un municipio. Estos mapas difieren de los mapas comunes, puesto que en estos se pueden visualizar características naturales, como ríos y montañas, o creaciones del hombre, como las divisiones políticas. Sin embargo, en los mapas temáticos se pueden apreciar elementos que ayudan a mejorar el propósito del tema y lo que quiere mostrar el autor. Como se mencionó en la lectura de Olaya, los mapas deben hacerse teniendo en cuenta el público a quien va dirigido y sus objetivos particulares. En este cuaderno se aprenderán a realizar varios tipos de mapas, junto con las características que los identifican.

Información

Primero se descargará la información que se necesita para hacer los mapas, esta información se encuentra en la página del Dane. En ese link se encontrarán los datos sobre Necesidades Básicas Insatisfechas (NBI), del Censo Nacional de Población y Vivienda 2018. Se descargará la información en formato .xlsx, que es de excel y allá se modificará para poder usar los datos necesarios del documento. Se guardó como NBI_Magdalen.xlsx.

Adecuación

Ahora se limpiará la memoria, se instalaran las bibliotecas necesarias para realizar los mapas y se cargaran en R.

rm(list=ls())                                           
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)                                         
library(rgeos)                                          
library(raster)                                         
library(sf)                                             
library(cartography)                                    
library(SpatialPosition)

Lectura de la información de NBI

Ahora se leerá el archivo de “Necesidades Básicas Insatisfechas” para el departamento del Magdalena.

nbi <- read_excel("C:/Users/vale_/Desktop/UNAL/6to semestre/GB/NBI_Magdalena.xlsx")

Ahora se revisarán los atributos de los datos

head(nbi)

Ahora se encontrará cuál es el municipio con el mayor porcentaje de NBI en el departamento:

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

nbi %>%                                                 
  slice(which.max(VIVIENDA)) -> max_vivienda                  
max_vivienda

Nueva Granada es un municipio ubicado al sur del departamento del Magdalena, que cuenta con una población aproximadamente de 22 mil habitantes. Y como se puede observar, el municipio es azotado por la falta de suministros de necesidades básicas. Este hecho se puede comprabar con una noticia de la revista (Semana)[https://www.semana.com/nacion/articulo/nueva-granada-magdalena-epicentro-de-trashumancia/444205-3] del 2015, en donde mencionan que este municipio ha sido decretado como zona de calamidad pública ya que sus habitantes no possen acceso a agua potable, gas domiciliario, cobertura de servicios telefónicos y en las zonas rurales solo el 8% cuenta con servicio de energía eléctrica.

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

Ahora se pondrá la tabla en orden descendente:

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

nbi %>%                                                 
  arrange(desc(INASISTENCIA)) -> desc_inas                        
desc_inas

Unir datos de NBI a municipios

Ahora se agregarán los municipios de Magdalena. Para esto se leerá la información usando la biblioteca de sf.

munic <- st_read("C:/Users/vale_/Desktop/UNAL/6to semestre/GB/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\vale_\Desktop\UNAL\6to semestre\GB\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
Simple feature collection with 30 features and 9 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -74.9466 ymin: 8.936489 xmax: -73.54184 ymax: 11.34891
CRS:            4326

Ahora se revisará que hay dentro del atributo MPIO_CCDGO:

head(munic$MPIO_CCDGO)

[1] 47541 47030 47058 47161 47170 47205
30 Levels: 47001 47030 47053 47058 47161 47170 47189 47205 47245 47258 47268 47288 47318 47460 47541 47545 47551 47555 47570 ... 47980

Se puede usar la función de left_join para unir el archivo de municipios y la información de NBI.

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

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

nbi2 <- nbi                                                                               
nbi2$CODIGO <- as.factor(nbi2$CODIGO)                                                     
head(nbi2)

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

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: -74.9466 ymin: 9.618984 xmax: -73.82743 ymax: 10.39318
CRS:            4326
          MUNICIPIO MPIO_CCDGO      NBI                       geometry
1           PEDRAZA      47541 34.40349 POLYGON ((-74.8832 10.21502...
2         ALGARROBO      47030 29.03316 POLYGON ((-74.10861 10.3915...
3          ARIGUANÍ      47058 44.00368 POLYGON ((-74.03381 9.99541...
4 CERRO SAN ANTONIO      47161 36.61789 POLYGON ((-74.85084 10.367,...
5           CHIVOLO      47170 44.52438 POLYGON ((-74.49589 10.2315...
6         CONCORDIA      47205 31.37073 POLYGON ((-74.81188 10.2748...

Note que Chivolo tiene NBI= 44.52, un valor que puede ser revisado en la tabla de NBI que se tenía antes, para asegurar que la unión quedó bien.

nbi_munic_new <- st_transform(nbi_munic, crs = 3116)            

Ejemplos de mapas temáticos

Ahora se usará el paquete de cartography con el cual se hacen mapas temáticos con la misma calidad visual que tiene el mapeo clásico o un software GIS. Sus funciones son intuitivas para los usuarios y son compatibles con R.
La biblioteca de cartografía usa objetos sf o sp para producir gráficos base, ya que el paquete se basa en funcionalidades sf.

^ Mapa base de OpenStreetMap y de símbolos proporcionales

Este tipo de mapas representa variables cuantitativas a través de símbolos cuyo tamaño está relacionado con el valor que se quiere representar en la variable.La variable visual tamaño, como lo menciona Olaya, es la única que puede representar las variables cuantitativas. Se recomienda usar siempre el mismo símbolo, el más común es el círculo, pero también están los cuadrados, e incluso las barras.
Es muy importante la relación de escalado que se debe utilizar para cada tamaño del símbolo, puesto que este define como se plasma la proporción que hay entre un valor y otro en el mapa. Teniendo cuidado de que entre esos tamaños no debe haber solapes al ponerlos sobre el mapa ya que crearía confución para el lector.
Ahora se mostrará un ejemplo de un mapa de Símbolos Proporcionales. Las funciones getTiles() y tilesLayer() descargan y muestran OpenStreetMap tiles.
propSymbolsLayer muestra símbolos con áreas proporcionales a una variable cuantitativa (por ejemplo, NBI). Hay varios símbolos que se pueden utilizar como círculos, cuadrados, barras, etc. El argumento inches es usado para editar el tamaño de los símbolos.

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

opar <- par(mar = c(0,0,1.2,10))                                 
tilesLayer(x = mun.osm)                                          
plot(st_geometry(nbi_munic_new), col = NA, border = "orange", add=TRUE)                                                        

propSymbolsLayer(                                               
  x = nbi_munic_new,                                             
  var = "NBI",                                                   
  inches = 0.18,                                                 
  col = "royalblue2",                                                
  legend.pos = "topright",                                       
  legend.title.txt = "Total NBI",                                
)                                                                
layoutLayer(title = "NBI Distribution in Magdalena",             
            sources = "Sources: DANE, 2018\n© OpenStreetMap",    
            author = "Valentina Valle",                          
            frame = TRUE, north = FALSE, tabtitle = TRUE)         

north(pos = "topleft")

De este mapa se puede concluir que el valor del índice de necesidades básicas insatisfechas en Magdalena es alto, en especial en la zona central del departamento.

^ Mapas de Coropleta

Estos mapas son utilizados para representar información gráfica de un SIG, en ellos se tiene una serie de áreas definidas (por ejemplo, un departamento) y cada una de ellas representa el valor de una variable que afecta a todo el área.
A pesar de que sirve como una buena herramienta para mapear, tiene sus fallas: 1) Puede haber cambios bruscos entre las representaciones de cada área ya que se toma un solo valor para una zona muy grande, pudiendo haber un cambio no tan brusco, sino uno continuo y este puede ser opacado por la representación antes mencionada. 2) El valor que se muestra por cada área podria representar solo una parte de la información y de esta manera hace que se pierda otra parte, que podría ser relevante para lo que quiere mostrar el autor en el mapa.
Se recomienda el uso de este tipo de mapas cuando se emplee el componente valor en vez del tono para representar información cuantitativa.
A continuación se realizará un mapa de Coropletas.

opar2 <- par(mar = c(0,0,1.2,0))                                 
par(bg="gray23")                                                 
plot(st_geometry(nbi_munic_new), col = NA, border = NA, bg = "wheat1")                                                       
choroLayer(                                                      
  x = nbi_munic_new,                                             
  var = "NBI",                                                   
  method = "geom",                                               
  nclass=10,                                                      
  col = carto.pal(pal1 = "wine.pal", n1 = 10),                    
  border = "white",                                              
  lwd = 2,                                                     
  legend.pos = "topright",                                      
  legend.title.txt = "NBI",                                     
  add = TRUE                                                     
)                                                                

layoutLayer(title = "NBI Distribution in Magdalena",             
            sources = "Source: DANE, 2018",                      
            author = "Valentina Valle",                          
            frame = TRUE, north = FALSE, tabtitle = TRUE, col="black")                                                     
north(pos = "topleft")

Con este mapa se puede visualizar mucho más fácil la falta de necesidades suplidas en este departamento, que es muy alta en aproximadamente el 90% del departamento.

^ Mapa de símbolos poporcionales 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 colorea para reflejar las modalidades de una segunda variable cualitativa. Se usa una combinación de los argumentos propSymbolsTypoLayer() y typoLayer().
Primero, se creará una variable cualitativaa que es muy necesaria para crear el mapa. Se usará la función mutate para esta tarea.

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

La nueva variable se llamó poverty. Su valor depende de los límites definidos por el usuario.

head(nbi_munic_2)

Simple feature collection with 6 features and 21 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 904750.7 ymin: 1555481 xmax: 1027439 ymax: 1641112
CRS:            EPSG:3116
  DPTO_CCDGO MPIO_CCDGO           MPIO_CNMBR                                    MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng
1         47      47541              PEDRAZA Decreto dptal 1322 del 19 dediciembre de 1908   325.4799      2017  MAGDALENA  0.9740441
2         47      47030            ALGARROBO             Ordenanza 008 de Junio 24 de 1999   406.3884      2017  MAGDALENA  1.2656800
3         47      47058             ARIGUANÍ      Ordenanza 14 Bis de Noviembre 30 de 1966  1131.7442      2017  MAGDALENA  1.9928810
4         47      47161 CERRO DE SAN ANTONIO                        Ordenanza 3038 de 1912   177.1961      2017  MAGDALENA  0.7003465
5         47      47170              CHIVOLO                Decreto 107 de Marzo 8 de 1974   536.5437      2017  MAGDALENA  1.3263733
6         47      47205            CONCORDIA         Ordenanza 007 del 24 de Junio de 1999   109.4837      2017  MAGDALENA  0.5443263
   Shape_Area COD_DEPTO     DEPTO COD_MUN         MUNICIPIO      NBI   MISERIA VIVIENDA SERVICIOS HACINAMIENTO INASISTENCIA  ECONOMIA
1 0.026847311        47 MAGDALENA     541           PEDRAZA 34.40349 10.490157 12.19448  7.213635     5.324349    10.291980 13.000396
2 0.033536941        47 MAGDALENA     030         ALGARROBO 29.03316  8.444102 10.16510  7.793480     7.065902     1.958864 13.327270
3 0.093278409        47 MAGDALENA     058          ARIGUANÍ 44.00368 17.073257 21.96648 22.588743    10.316787     2.658747  8.736388
4 0.014622408        47 MAGDALENA     161 CERRO SAN ANTONIO 36.61789 14.417344 14.03794 15.121951    11.750678     2.319783 13.658537
5 0.044254449        47 MAGDALENA     170           CHIVOLO 44.52438 18.711555 32.40334 20.672232     2.746046     1.620167  9.951670
6 0.009033178        47 MAGDALENA     205         CONCORDIA 31.37073  6.972420 10.21589  8.501188     1.497779     1.332507 18.706745
                        geometry      poverty
1 POLYGON ((911720.4 1621516,...         High
2 POLYGON ((996594.4 1640935,...         High
3 POLYGON ((1004791 1597116, ...         High
4 POLYGON ((915306.4 1638318,...         High
5 POLYGON ((954161.3 1623266,... Intermediate
6 POLYGON ((919550.3 1628117,... Intermediate

library(sf)                                                     
library(cartography)
opar3 <- par(mar = c(0,0,1.2,0))                                
plot(st_geometry(nbi_munic_2), col="#f2efe9", border="#b38e43", bg = "antiquewhite1",   
     lwd = 0.8)                                       
propSymbolsTypoLayer(                                            
  x = nbi_munic_2,                                               
  var = "NBI",                                                   
  inches = 0.3,                                                  
  symbols = "square",                                            
  border = "white",                                              
  lwd = .5,                                                      
  legend.var.pos = "topright", c(1050000, 135000),                
  legend.var.title.txt = "NBI",
  var2 = "poverty",                                              
  legend.var2.values.order = c("Extreme", "High",               
                               "Intermediate"),                  
  col = carto.pal(pal1 = "turquoise.pal", n1 = 3),                
  legend.var2.pos = "topleft", c(1050000, 1200000),             
  legend.var2.title.txt = "Poverty"
)                                                                

layoutLayer(title="NBI Distribution in Magdalena",               
            author = "Valentina Valle",                         
            sources = "Source: DANE, 2018",                      
            scale = 10, tabtitle = TRUE, frame = TRUE)           
north(pos = "topright")

Aquí se puede apreciar con mayor facilidad a través del uso de colores que encasillan el NBI en variables cualitativas. También se puede ver un poco de cambios abruptos en los colores, lo que permite ver un cambio abrupto entre municipios, así la diferencia no sea muy grande.

^ Mapa de etiquetas

Ahora se intentará combinar las funciones choroLayer y labelLayer.

library(sf)                                                      
library(cartography)                                             
opar4 <- par(mar = c(0,0,1.2,0))                                
par(bg="grey25")                                                 
plot(st_geometry(nbi_munic_2), col = "black", border = "black",                                                 
     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 = "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 Magdalena",                             
  sources = "Source: DANE, 2018",                                
  author = "Valentina Valle",                                    
  frame = TRUE,                                                  
  north = TRUE,                                                  
  tabtitle = TRUE,                                              
  theme = "taupe.pal"                                            
)

Con el mapa de Municipios queda mucho más fácil ver cuales son lo municipios en los que el índice de NBI es más alto, que en este caso son Nueva Granada y Sabanas de San Ángel.

^ Mapa de isopletas

Los mapas de isopletas son basados en la suposición de que el fenómeno a ser representado tiene una distribución continua. smoothLayer() depende en gran medida del paquete SpatialPosition. La función usa 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 un mapa de isopleta.
Ahora se usara otra información para realizar el mapa de isopleta. Se utilizará la información usada en el cuaderno de R “Estadísticas Municipales de la Agricultura en Magdalena”, enfocada a la producción de yuca de este departamento.

crops2018 <- read_excel("C:/Users/vale_/Desktop/UNAL/6to semestre/GB/EVA_MAGDALENA.xlsx")

head(crops2018)

Se filtrará para usar solo la información sobre yuca en 2018:

crops2018 %>%                                                    
  filter(CULTIVO == "YUCA") -> yuca2018        

Ahora se creará un nuevo atributo que unirá los códigos de los municipios:

yuca2018$MPIO_CCDGO <- as.character(yuca2018$COD_MUN)         

munic$MPIO_CCDGO <- as.character(munic$MPIO_CCDGO)

Y la unión:

yuca_munic = left_join(munic, yuca2018, by="MPIO_CCDGO")     

head(yuca_munic)

Simple feature collection with 6 features and 23 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -74.9466 ymin: 9.618984 xmax: -73.82743 ymax: 10.39318
CRS:            4326
  DPTO_CCDGO MPIO_CCDGO           MPIO_CNMBR                                    MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng
1         47      47541              PEDRAZA Decreto dptal 1322 del 19 dediciembre de 1908   325.4799      2017  MAGDALENA  0.9740441
2         47      47030            ALGARROBO             Ordenanza 008 de Junio 24 de 1999   406.3884      2017  MAGDALENA  1.2656800
3         47      47058             ARIGUANÍ      Ordenanza 14 Bis de Noviembre 30 de 1966  1131.7442      2017  MAGDALENA  1.9928810
4         47      47161 CERRO DE SAN ANTONIO                        Ordenanza 3038 de 1912   177.1961      2017  MAGDALENA  0.7003465
5         47      47170              CHIVOLO                Decreto 107 de Marzo 8 de 1974   536.5437      2017  MAGDALENA  1.3263733
6         47      47205            CONCORDIA         Ordenanza 007 del 24 de Junio de 1999   109.4837      2017  MAGDALENA  0.5443263
   Shape_Area COD_DEP DEPARTAMENTO COD_MUN            MUNICIPIO                 GRUPO SUBGRUPO CULTIVO YEAR Area_Siembra Area_Cosecha
1 0.026847311      47    MAGDALENA   47541              PEDRAZA TUBERCULOS Y PLATANOS     YUCA    YUCA 2018         1000         1000
2 0.033536941      47    MAGDALENA   47030            ALGARROBO TUBERCULOS Y PLATANOS     YUCA    YUCA 2018          650          650
3 0.093278409      47    MAGDALENA   47058             ARIGUANI TUBERCULOS Y PLATANOS     YUCA    YUCA 2018          890          890
4 0.014622408      47    MAGDALENA   47161 CERRO DE SAN ANTONIO TUBERCULOS Y PLATANOS     YUCA    YUCA 2018         1000         1000
5 0.044254449      NA         <NA>      NA                 <NA>                  <NA>     <NA>    <NA>   NA           NA           NA
6 0.009033178      47    MAGDALENA   47205            CONCORDIA TUBERCULOS Y PLATANOS     YUCA    YUCA 2018          800          800
  Produccion Rendimiento           ESTADO CICLO                       geometry
1       8000           8 TUBERCULO FRESCO ANUAL POLYGON ((-74.8832 10.21502...
2       5200           8 TUBERCULO FRESCO ANUAL POLYGON ((-74.10861 10.3915...
3       7120           8 TUBERCULO FRESCO ANUAL POLYGON ((-74.03381 9.99541...
4       8000           8 TUBERCULO FRESCO ANUAL POLYGON ((-74.85084 10.367,...
5         NA          NA             <NA>  <NA> POLYGON ((-74.49589 10.2315...
6       6400           8 TUBERCULO FRESCO ANUAL POLYGON ((-74.81188 10.2748...

rep_yuca <- st_transform(yuca_munic, crs = 3116)               

¡¡A mapear!!

opar5 <- par(mar = c(0,0,1.2,0))                                
plot(st_geometry(rep_yuca), col = NA, border = "black", bg = "grey75")                                                        
smoothLayer(                                                     
  x = rep_yuca,                                                 
  var = 'Produccion',                                            
  typefct = "exponential",                                      
  span = 400000,                                                
  beta = 1,                                                      
  nclass = 10,                                                   
  col = carto.pal(pal1 = 'purple.pal', n1 = 12),                 
  border = "grey",                                               
  lwd = 0.1,                                                     
  mask = rep_yuca,                                             
  legend.values.rnd = -3,                                        
  legend.title.txt = "Produccion",                              
  legend.pos = "topright",                                       
  add=TRUE                                                      
)                                                               

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 = "Distribucion de la Produccion de Yuca en Magdalena",                                              
            sources = "Sources: DANE and MADR, 2018",           
            author = "Valentina Valle",                          
            frame = FALSE, north = FALSE, tabtitle = TRUE, theme = "green.pal")                                                   

north(pos = "topleft")

Como se mencionó anteriormente, el índice de NBI se concentra en el centro del departamento, lo que se puede confirmar con el anterior mapa.

Guardar Mapas

Ahora se realizará otro mapa de la producción de yuca en 2018. Esta vez usando mapas de símbolos proporcionales y coropletas. El resultado se guardará como un archivo .png.
El siguiente chunk no muestra un mapa. En cambio, escribe el mapa en el nombre del archivo yuca_2018.png debajo del directorio de trabajo.

png("C:/Users/vale_/Desktop/UNAL/6to semestre/GB/yuca_2018.png", width = 2048, height = 1526)                                     
opar6 <- par(mar = c(0,0,5,5))                                   
plot(st_geometry(rep_yuca), col="mistyrose", border="darkseagreen4",                                          
     bg = "white", lwd = 0.6)                                    
propSymbolsChoroLayer(x = rep_yuca, var = "Produccion", var2 = "Rendimiento",                                                   
                      col = carto.pal(pal1 = "orange.pal", n1 = 3, 
                                      pal2 = "blue.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 = "Produccion de Yuca en 2018",                                                      
                      legend.var.style = "e")                    
labelLayer(                                                      
  x = rep_yuca,                                           
  txt = "MPIO_CNMBR",                                            
  col= "black",                                                  
  cex = 1.0,                                                     
  font = 4,                                                      
  halo = FALSE,                                                  
  bg = "white",                                                  
  r = 0.1,                                                       
  overlap = FALSE,                                              
  show.lines = FALSE                                             
)                                                                
layoutLayer(title="Produccion y Rendimiento de Yuca en Magdalena, 2018",                                                  
            author = "Valentina Valle",                
            sources = "Sources: MADR & DANE, 2018",       
            scale = 50, tabtitle = FALSE, frame = TRUE)          
north(pos = "topleft")                                 
title(main="Produccion y Rendimiento de Yuca en Magdalena, 2018", cex.main=3,                                               
            sub= "Source: MADR & DANE, 2018", cex.sub=2)        
graticule = TRUE                                                
par(opar6)                                                      
dev.off() 

null device 
          1 

Una vez que se guarda el mapa puede agregarlo al cuaderno usando el siguiente código, el cual incluye la imagen en este sitio:

knitr::include_graphics('C:/Users/vale_/Desktop/UNAL/6to semestre/GB/yuca_2018.png')

Con este último mapa se puede observar que, en general en el departamento, la producción de yuca es alto y el mejor rendimiento se puede encontrar en la zona centro-superior.

sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Colombia.1252  LC_CTYPE=Spanish_Colombia.1252    LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C                     
[5] LC_TIME=Spanish_Colombia.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] SpatialPosition_2.0.1 cartography_2.4.1     sf_0.9-0              raster_3.0-12         rgeos_0.5-2           readxl_1.3.1         
 [7] forcats_0.5.0         stringr_1.4.0         dplyr_0.8.5           purrr_0.3.3           readr_1.3.1           tidyr_1.0.2          
[13] tibble_2.1.3          ggplot2_3.3.0         tidyverse_1.3.0       sp_1.4-1             

loaded via a namespace (and not attached):
 [1] httr_1.4.1         jsonlite_1.6.1     modelr_0.1.6       shiny_1.4.0        assertthat_0.2.1   cellranger_1.1.0   yaml_2.2.1        
 [8] pillar_1.4.3       backports_1.1.5    lattice_0.20-38    glue_1.3.1         digest_0.6.25      promises_1.1.0     rvest_0.3.5       
[15] colorspace_1.4-1   htmltools_0.4.0    httpuv_1.5.2       slippymath_0.3.1   pkgconfig_2.0.3    broom_0.5.5        haven_2.2.0       
[22] xtable_1.8-4       scales_1.1.0       later_1.0.0        farver_2.0.3       generics_0.0.2     withr_2.1.2        cli_2.0.2         
[29] magrittr_1.5       crayon_1.3.4       mime_0.9           evaluate_0.14      fs_1.3.2           fansi_0.4.1        nlme_3.1-144      
[36] xml2_1.2.5         lwgeom_0.2-1       class_7.3-15       rsconnect_0.8.16   tools_3.6.3        hms_0.5.3          lifecycle_0.2.0   
[43] munsell_0.5.0      reprex_0.3.0       isoband_0.2.0      compiler_3.6.3     e1071_1.7-3        rlang_0.4.5        classInt_0.4-2    
[50] units_0.6-6        grid_3.6.3         rstudioapi_0.11    htmlwidgets_1.5.1  crosstalk_1.0.0    labeling_0.3       base64enc_0.1-3   
[57] rmarkdown_2.1      gtable_0.3.0       codetools_0.2-16   DBI_1.1.0          R6_2.4.1           lubridate_1.7.4    knitr_1.28        
[64] rgdal_1.4-8        fastmap_1.0.1      KernSmooth_2.23-16 stringi_1.4.6      parallel_3.6.3     Rcpp_1.0.3         vctrs_0.2.4       
[71] png_0.1-7          leaflet_2.0.3      dbplyr_1.4.2       tidyselect_1.0.0   xfun_0.12