Cartografía temática en el Amazonas

1. Introducción

El departamento del Amazonas es uno de los 32 departamentos de Colombia y está ubicado en la zona sur del país, cuenta con 10 municipios y su capital es Leticia. Este departamento cuenta con una de las zonas más ampliamente biodiversas del mundo gracias a que posee parte de la selva amazónica (la más grande del mundo) y adicionalmente tambien tiene una gran diversidad de culturas indígenas. Su ubicación es vecina a los países de Brasil y Perú. El Amazonas tiene una superficie de casi 110.000 km², está ampliamente cubierto de selva y posee varios ríos que son tributarios del río amazonas. La economía se basa en la extracción de diferentes cultivos como maíz, plátano, arroz, caña de azúcar, ñame, aguacate, yuca, cacao, piña, etc.

Este notebook se basa en mostrar representaciones gráficas de los diferentes mapas temáticos basados en las clases de la asignatura de Geomática básica, usando como referencia el departamento del Amazonas y su producción agrícola, por medio del programa Rstudio.

2. Cartografía temática

La cartógrafia temática es un tipo de mapa, clasificado segun su tipo de información ya que se basa en la representación de una variable espacial de carácter físico, social, político, etc. De igual forma, muestra atributos o estadísticas que permiten un mejor rendimiento de lo que pasa en diversas regiones. La cartografía temática cuenta con diversas técnicas de mapeo temática, en ese informe se mostrarán 4 de ellas:

2.1 Mapas de coropletas

Estos mapas son usados principalmente para representar información geográfica en un SIG. En estos mapas se pueden observar áreas definidas las cuales tienen un valor de una variable cada una. Además, cada área definida corresponde a una unidad espacial y el valor que posee resume la variable dentro de dicha área. Para los mapas de coropletas es muy importante tener en cuenta la correcta división de valores de las variables, usualmente, el uso de distintos colores es el método más usado.\(^1\)

2.2 Mapas de símbolos proporcionales

Este tipo de mapas representa variables cuantitativas a través de símbolos en donde los tamaños están en relación con el valor que se quiere representar para la variable otorgada. Los símbolos que tienden a usarse para este tipo de mapas son círculos, cuadrados y triángulos.\(^2\)

2.3 Mapas de isolíneas

Las isolíneas son ‘líneas’ con un valor constante asociado a todos sus puntos. Un mapa de isolíneas representa valores continuos que son propios de la variable.\(^3\) En este caso, el dato está en todos los puntos del espacio y de forma continua, pero sólo se mide en los puntos de control. Las líneas se trazan con intervalos constantes y pueden tener colores secuenciales.

2.4 Mapas de puntos

Estos mapas se emplean principalemnte para representar variables que representen cantidad. Esas cantidades son trazadas mediante la repetición de puntos en número proporcional a su tamaño. Cada punto representa un valor y el conjunto de todos esos puntos suma la totalidad a representar.

3. Datos

Los datos de este informe hacen referencia a las Necesidades Básicas insatisfechas (NBI), las cuales fueron tomadas del Geoportal DANE. Para un correcto uso de los datos, al descargarlos, se edito desde el programa de excel y se eliminaron los datos que en este documento no son relevantes conocerlos; se adecuaron los datos al departamento del Amazonas.

4. Paquetes y librerías necesarias

Primero borramos la memoria:

rm(list = ls())

Instalamos los paquetes necesarios y que no se han instalado con anterioridad:

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)

Cargamos las librerias:

library(tidyverse)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.0     v purrr   0.3.4
v tibble  3.0.0     v dplyr   0.8.5
v tidyr   1.0.2     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x tidyr::extract() masks raster::extract()
x dplyr::filter()  masks stats::filter()
x dplyr::lag()     masks stats::lag()
x dplyr::select()  masks raster::select()

library(readxl)
library(rgeos)

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)

library(sf)

Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

library(cartography)
library(SpatialPosition)

5. Leer los datos de las NBI

En este paso ejecutamos el código siguiente para leer el documento de excel mencionado previamente, y otorgarle un nombre,en este caso se le dio el nombre de nbi.

nbi <- read_excel("C:/Users/paula/Desktop/Geomatica/Cartografia_tematica/CNPV-2018-NBI.xlsx")

Ver las primeras 6 filas del documento para verificar que los datos sean correctos:

head(nbi)

Para lo siguiente, se ejecuta el código para filtrar los datos con mayores niveles de NBI y ponerlos en una nueva tabla. Al observar esos datos se puede concluir que el municipio de La Victoria tiene un NBI de casi 91, lo que significa que el nivel de pobreza es bastante alto.

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

Aquí se tiene la misma metodología con la diferencia de que se filtraron los datos con los menores niveles de NBI. Como se observa, en Leticia los niveles de NBI son de aproximadamente 30 no es un nivel muy bajo, lo que concierne a que eb general el departamento del Amazonas tiene un nivel alto de pobreza.

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

min_nbi

Al correr el siguiente código se observan los datos de la tabla nbi pero con los datos de NBI en forma descendente.

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

desc_nbi

6. Unión de datos de NBI a municipios

En el siguiente paso, se cargaran los datos de los municipios del Amazonas que fueron optenidos en la realización de un Notebook anterior.

munic <- st_read("C:/Users/paula/Desktop/Geomatica/Datos_de _elevacion_en_R/Amazonas_datos_DANE/91_AMAZONAS/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\paula\Desktop\Geomatica\Datos_de _elevacion_en_R\Amazonas_datos_DANE\91_AMAZONAS\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
Simple feature collection with 11 features and 9 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -74.39634 ymin: -4.229406 xmax: -69.3955 ymax: 0.09725686
CRS:            4326

Se ejecuta el siguiente comando para observar unicamente la columna de los códigos de los municipios MPIO_CCDGO para comprobar que esos datos son equivalentes a los códigos de la tabla nbi.

head(munic$MPIO_CCDGO)

[1] 91001 91263 91405 91407 91430 91460
11 Levels: 91001 91263 91405 91407 ... 91798

Realizamos la unión de esas dos columnas con la funcion de left_join.

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

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

Creamos una nueva tabla con las columnas de MUNICIPIO, MPIO_CCDGO y NBI; para comprobar que el join fue exitoso.

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: -73.84132 ymin: -4.229406 xmax: -69.3955 ymax: 0.09725686
CRS:            4326
              MUNICIPIO MPIO_CCDGO      NBI
1               LETICIA      91001 26.99132
2      EL ENCANTO (ANM)      91263 37.50764
3     LA CHORRERA (ANM)      91405 35.81580
4      LA PEDRERA (ANM)      91407 78.83191
5     LA VICTORIA (ANM)      91430 90.96386
6 MIRITI - PARANÁ (ANM)      91460 85.76471
                        geometry
1 POLYGON ((-69.74085 -3.0026...
2 POLYGON ((-73.24647 -1.5029...
3 POLYGON ((-73.72362 -0.3804...
4 POLYGON ((-70.45577 -0.4957...
5 POLYGON ((-71.18431 0.09705...
6 POLYGON ((-71.48046 -0.0105...

nbi_munic_new <- st_transform(nbi_munic, crs = 3116)

7. Ejemplos de mapas temáticos

Para llevar a cabo estos mapas temáticos fue necesario hacer uso del paquete cartography ya que la calidad de los mapas en esta librería es adecuada y buena. Para poder hacer uso es esta librería se necesitan objetos espaciales que en este caso es sf (Simple Features).

7.1 Símbolos proporcionales usando OpenStreetMap

En los siguientes comandos se hacen uso distintas funciones como getTiles() que descarga los datos y tiles Layer() que están encargadas de mostrar los datos. Para visualizar los simbolos con areas proporcionales de una variable cuantitativa se usa la funcion propSymbolsLayer()

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

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

propSymbolsLayer(
  x = nbi_munic_new,
  var = "NBI",
  inches = 0.25,   # inches hace referencia al tamaño de las figuras
  col = "peru",
  legend.pos = "topright",
  legend.title.txt = "NBI totales"
)
layoutLayer(title = "Distribución de las NBI en el Amazonas",
            sources = "Fuente: DANE, 2018\n© OpenStreetMap",
            author = "Paula Virgüez",
            frame = TRUE, north = FALSE, tabtitle = TRUE)

north(pos = "topleft")

7.2 Coropletas

En este caso usaremos la función choroLayer() la cual está encargada de mostrar los mapas en coropletas tal como se explicó previamente.

opar <- par(mar = c (0,0,1.2,0))
par(bg = "grey90")
plot(st_geometry(nbi_munic_new), col = NA, border = NA, bg = "#85C1E9")
choroLayer(
  x= nbi_munic_new,
  var = "NBI",
  method = "geom",
  nclass = 5,
  col = carto.pal(pal1 = "kaki.pal", n1 = 5), # carto.pal() es usado para definir la paleta de colores que se va a emplear en el mapa
  border = "white",
  lwd = 0.5,
  legend.pos = "topright",
  legend.title.txt = "NBI",
  add = TRUE
)

layoutLayer(title = "Distribución de las NBI en el Amazonas",
            sources = "Fuente: DANE, 2018.",
            author = "Paula Virgüez",
            frame = TRUE, north = FALSE, tabtitle = TRUE, col = "black")
north(pos = "topleft")

7.3 Símbolos proporcionales y mapa de tipología

En este mapa los símbolos tienen un tamaño que es proporcional a una variable cuantitativa que se quiere resaltar y coloreados en funcion del valor de una segunda variable que es cualitativa.

A continuación, se hace uso de la función mutate() del paquete dplyr para crear una nueva columna, en este caso se adiciona una columna nueva con información de la pobreza en el amazonas pero diferenciando los datos con intervalos diferentes.

nbi_munic_2 <- dplyr::mutate(nbi_munic_new, Pobreza = ifelse(MISERIA > 20, "Extremo",
                                                            ifelse(HACINAMIENTO > 5,
                                                                   "Alto", "Intermedio")))

Como podemos observar, en la tabla aparece la nueva columna de “Pobreza”

head(nbi_munic_2)

Simple feature collection with 6 features and 21 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 1026292 ymin: 22875.46 xmax: 1521684 ymax: 502535.4
CRS:            EPSG:3116
  DPTO_CCDGO MPIO_CCDGO
1         91      91001
2         91      91263
3         91      91405
4         91      91407
5         91      91430
6         91      91460
                 MPIO_CNMBR
1                   LETICIA
2                EL ENCANTO
3               LA CHORRERA
4                LA PEDRERA
5       LA VICTORIA (Pacoa)
6 MIRITÍ-PARANÁ (Campoamor)
                      MPIO_CRSLC MPIO_NAREA
1  Decreto 352 de Feb 20 de 1964   6183.129
2 Decreto 274 de Mayo 28 de 1953  10897.481
3 Decreto 274 de Mayo 28 de 1953  12726.929
4 Decreto 274 de Mayo 28 de 1953  13670.317
5      ORD 12 DE JULIO 9 DE 1996   1432.795
6 Decreto 274 de Mayo 28 de 1953  16866.120
  MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area
1      2017   AMAZONAS   3.880001  0.5007018
2      2017   AMAZONAS   8.052697  0.8853254
3      2017   AMAZONAS   9.105101  1.0335830
4      2017   AMAZONAS   7.376291  1.1051976
5      2017   AMAZONAS   2.483170  0.1160863
6      2017   AMAZONAS   6.448395  1.3668663
  COD_DEP DEPARTAMENTO COD_MUN
1      91     AMAZONAS     001
2      91     AMAZONAS     263
3      91     AMAZONAS     405
4      91     AMAZONAS     407
5      91     AMAZONAS     430
6      91     AMAZONAS     460
              MUNICIPIO      NBI   MISERIA
1               LETICIA 26.99132  6.270907
2      EL ENCANTO (ANM) 37.50764 11.178986
3     LA CHORRERA (ANM) 35.81580  9.778283
4      LA PEDRERA (ANM) 78.83191 42.364672
5     LA VICTORIA (ANM) 90.96386 66.265060
6 MIRITI - PARANÁ (ANM) 85.76471 63.058824
   VIVIENDA SERVICIOS HACINAMIENTO
1  4.275949  9.620885     14.31376
2  1.588271 21.624924     17.77642
3  4.263786 22.171688     14.66742
4 21.623932 71.709402     38.14815
5 48.795181 83.132530     22.89157
6 21.294118 78.117647     63.17647
  INASISTENCIA DEP_ECONOMICA
1     2.850633      3.994764
2     4.153940      5.192425
3     3.411029      3.865833
4     7.663818     16.039886
5    37.349398     22.891566
6     4.941176      9.647059
                        geometry Pobreza
1 POLYGON ((1482558 158792.9,...    Alto
2 POLYGON ((1092482 325567.3,...    Alto
3 POLYGON ((1039394 449700.2,...    Alto
4 POLYGON ((1403426 436838.2,... Extremo
5 POLYGON ((1322206 502513.1,... Extremo
6 POLYGON ((1289202 490602.7,... Extremo

Plotear el mapa

library(sf)
library(cartography)
opar <- par(mar = c(0,0,1.2,0))
  plot(st_geometry(nbi_munic_2), col = "#E5E8E8", border = "#B2BABB", bg = "#EBF5FB",
     lwd = 0.5)
propSymbolsTypoLayer(
  x = nbi_munic_2,
  var = "NBI", # Variable cuantitativa
  inches = 0.5,
  symbols = "square",
  border = "white",
  lwd = .5,
  legend.var.pos = "topright",
  legend.var.title.txt = "NBI",
  var2 = "Pobreza", # Variable cualitativa
  legend.var2.values.order = c("Extremo", "Alto",
                               "Intermedio"),
  col = carto.pal(pal1 = "pastel.pal", n1 = 3),
  legend.var2.pos = c(820000,100000),
  legend.var2.title.txt = "Pobreza"
)

layoutLayer(title = "Distribución de las NBI en el Amazonas",
            author = "Paula Virgüez",
            sources = "Fuente: DANE, 2018.",
            scale = 100, tabtitle = TRUE, frame = FALSE, postitle = "right",
            col = "white", coltitle = "black")
north(pos = "topleft")

Como se puede observar, principalmente hacia el norte del departamento se encuentran municipios con mayor pobreza y también hacia el sur.

7.4 Mapas de etiquetas

En este caso se combinarán las funciones de choroLayer() y labelLayer. Cambié los colores de los municipios, el fondo de la imagen, la posición del símbolo Norte, la paleta de colores, etc.

opar <- par(mar = c(0,0,1.2,0))
par (bg = "grey25")
plot(st_geometry(nbi_munic_2), col = "#e4e9de", border = "#000000",
     bg = "white", lwd = 0.5)
choroLayer(
  x= nbi_munic_new,
  var = "NBI",
  method = "geom",
  nclass = 5,
  col = carto.pal(pal1 = "blue.pal", n1 = 6),
  border = "grey10",
  lwd = 0.5,
  legend.pos = c(1600000, 100000),
  legend.title.txt = "NBI",
  add = TRUE
)

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

layoutLayer(
  title = "Municipios del Amazonas",
  sources = "Fuente: DANE, 2018.",
  author = "Paula Virgüez",
  frame = TRUE,
  north = FALSE,
  tabtitle = TRUE,
  theme = "blue.pal"
)

north(pos = "topleft")

7.5 Mapa de Isopletas

Para este mapa se usarán los datos estadísticos de la producción de yuca del año 2018 en el Amazonas. Estos datos fueron usados en un Notebook anterior, sin embargo, se modificó el archivo eliminando todos los valores no relacionados al año 2018 para no tener inconvenientes posteriormente.

Cargar los datos:

cultivos2018 <- read_excel("C:/Users/paula/Desktop/Geomatica/Cartografia_tematica/EVA_Amazonas_2018.xlsx")

head(cultivos2018)

En este caso, opté por usar este cultivo ya que, basandome en el Notebook anterior, la yuca es el cultivo que tuvo los valores más altos de producción en el año 2018 en el Amazonas, aunque solo hubieron datos en dos municipios del departamento, los valores de producción fueron los más altos de los que se tiene registro, y son ideales para poder tener una mejor representación del mapa de isopletas. Aquí se usó la función filter() justamente para filtrar unicamente los datos del cultivo de yuca:

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

Se cambian los tipos de datos de las columnas CODIGO_MUN y TEMP para poder realizar el correspondiente join. Cabe mencionar que en estos pasos no se realizó la adición o edición de datos de las columnas (agregando o quitando 0) ya que los datos de las columnas de ambas tablas coinciden perfectamente sin haber hecho ningun cambio previo.

yuca2018$TEMP <- as.character(yuca2018$CODIGO_MUN)

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

Realizar el join con las tablas munic y yuca2018 con la columna MPIO_CCDGO:

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

Column `MPIO_CCDGO` joining factors with different levels, coercing to character vector

Producto del join:

head(yuca_munic)

Simple feature collection with 6 features and 24 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -73.84132 ymin: -4.229406 xmax: -69.3955 ymax: 0.09725686
CRS:            4326
  DPTO_CCDGO MPIO_CCDGO
1         91      91001
2         91      91263
3         91      91405
4         91      91407
5         91      91430
6         91      91460
                 MPIO_CNMBR
1                   LETICIA
2                EL ENCANTO
3               LA CHORRERA
4                LA PEDRERA
5       LA VICTORIA (Pacoa)
6 MIRITÍ-PARANÁ (Campoamor)
                      MPIO_CRSLC MPIO_NAREA
1  Decreto 352 de Feb 20 de 1964   6183.129
2 Decreto 274 de Mayo 28 de 1953  10897.481
3 Decreto 274 de Mayo 28 de 1953  12726.929
4 Decreto 274 de Mayo 28 de 1953  13670.317
5      ORD 12 DE JULIO 9 DE 1996   1432.795
6 Decreto 274 de Mayo 28 de 1953  16866.120
  MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area
1      2017   AMAZONAS   3.880001  0.5007018
2      2017   AMAZONAS   8.052697  0.8853254
3      2017   AMAZONAS   9.105101  1.0335830
4      2017   AMAZONAS   7.376291  1.1051976
5      2017   AMAZONAS   2.483170  0.1160863
6      2017   AMAZONAS   6.448395  1.3668663
  CODIGO_DEP DEPARTAMENTO CODIGO_MUN MUNICIPIO
1         91     AMAZONAS      91001   LETICIA
2         NA         <NA>         NA      <NA>
3         NA         <NA>         NA      <NA>
4         NA         <NA>         NA      <NA>
5         NA         <NA>         NA      <NA>
6         NA         <NA>         NA      <NA>
                  GRUPO SUBGRUPO CULTIVO YEAR
1 TUBERCULOS Y PLATANOS     YUCA    YUCA 2018
2                  <NA>     <NA>    <NA>   NA
3                  <NA>     <NA>    <NA>   NA
4                  <NA>     <NA>    <NA>   NA
5                  <NA>     <NA>    <NA>   NA
6                  <NA>     <NA>    <NA>   NA
  AREA_SEMBRADA AREA_COSECHADA PRODUCCION
1          1100           1100       7500
2            NA             NA         NA
3            NA             NA         NA
4            NA             NA         NA
5            NA             NA         NA
6            NA             NA         NA
  RENDIMIENTO           ESTADO CICLO  TEMP
1        6.82 TUBERCULO FRESCO ANUAL 91001
2        <NA>             <NA>  <NA>  <NA>
3        <NA>             <NA>  <NA>  <NA>
4        <NA>             <NA>  <NA>  <NA>
5        <NA>             <NA>  <NA>  <NA>
6        <NA>             <NA>  <NA>  <NA>
                        geometry
1 POLYGON ((-69.74085 -3.0026...
2 POLYGON ((-73.24647 -1.5029...
3 POLYGON ((-73.72362 -0.3804...
4 POLYGON ((-70.45577 -0.4957...
5 POLYGON ((-71.18431 0.09705...
6 POLYGON ((-71.48046 -0.0105...

Como se puede observar, la gran mayoría de los datos no tienen valores asignados (NA), pero no es por ningun error, los datos están correctos. Esto ocurrió por que, como se mencinó previamente, en el departamento del amazonas en el año 2018 solo hubo datos de producción de yuca en 2 de sus departamentos.

Reporyectamos los municipios para poder realizar el mapa.

rep_yuca <- st_transform(yuca_munic, crs = 3116)

Ploteamos el mapa de isopletas:

opar <- par(mar = c(0,0,1.2,0))

plot(st_geometry(rep_yuca), col = NA, border = "black", bg = "lavender")
smoothLayer(
  x = rep_yuca,
  var = 'PRODUCCION',
  typefct = "exponential",
  span = 70000,
  beta = 2,
  nclass = 10,
  col = carto.pal(pal1 = 'kaki.pal', n1= 10),
  border = NA,
  lwd = 0.1,
  mask = rep_yuca,
  legend.values.rnd = -3,
  legend.title.txt = "Producción",
  legend.pos = "topright",
  add = TRUE
)

text(x = 900000, y = 250000, cex = 0.6, adj = 0, font = 3, labels = "Función de distancia:\n- Tipo = Exponencial\n- beta = 2\n- span = 70 km")

layoutLayer(title = "Distribución de la producción de yuca en el Amazonas",
            sources = "Fuentes: DANE & MADR, 2018",
            author = "Paula Virgüez",
            frame = FALSE, north = FALSE, tabtitle = TRUE, theme = "kaki.pal")

north(pos = "topleft")

8. Guardar mapas

Aquí se guardará el archivo en una imagen de png con la resolución que se señala en el comando y las características para crear el archivo como colores, tamaños, posiciones, etc. Como se ha venido mencionando, gran parte de los datos son sin valor NA, entonces para este caso lo que a continuación se hizo fue reemplazar las celdas vacías y otorgarles el valor 0 a cada una de ellas, para que al plotear el mapa los demás datos fueran observables.

Para esto, se uso en las columnas PRODUCCIÓN y RENDIMIENTO con la función replace_na() de la librería tidyr:

rep_yuca %>% replace_na(list(PRODUCCION = list(0)))
rep_yuca %>% replace_na(list(RENDIMIENTO = list(0)))

rep_yuca <- rep_yuca %>% mutate(PRODUCCION = replace_na(PRODUCCION, 0))
rep_yuca <- rep_yuca %>% mutate(RENDIMIENTO = replace_na(PRODUCCION, 0))

Plotear el mapa, y guardar el archivo en png:

png("C:/Users/paula/Desktop/Geomatica/Cartografia_tematica/yuca_2018.png", width = 2048, height = 1526)
opar <- par(mar = c(0,0,5,5))
plot(st_geometry(rep_yuca), col = "lightcoral", border = "mistyrose2",
     bg = "mistyrose", lwd = 0.6)

propSymbolsChoroLayer(x = rep_yuca, var = "PRODUCCION", var2 = "RENDIMIENTO",
                      col = carto.pal(pal1 = "harmo.pal", n1 = 6),
                      inches = 0.8, method = "q6",
                      border = NA, lwd = 1,
                      legend.title.cex = 1.5, 
                      legend.values.cex = 1.0,
                      legend.var.pos = "topright",
                      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 yuca en 2018",
                      legend.var.style = "e")

labelLayer(
  x = rep_yuca,
  txt = "MPIO_CNMBR",
  col = "white",
  cex = 1.9,
  font = 4,
  halo = TRUE,
  bg = "black",
  r = 0.1,
  overlap = FALSE,
  show.lines = FALSE
)

layoutLayer(title = "Producción y rendimiento de yuca en el Amazonas, 2018",
            author = "Paula Virgüez",
            sources = "Fuentes: MADR & DANE, 2018",
            scale = 50, tabtitle = FALSE, frame = TRUE)

north(pos = "topleft")

title(main = "Producción y rendimiento de yuca en el Amazonas, 2018", cex.main = 3,
      sub = "Fuente: MADR & DANE, 2018", cex.sub=2)

graticule = TRUE

par(opar)

dev.off()

null device 
          1 

Mostrar la imagen .png realizada

library(ggplot2)
library(knitr)
include_graphics("C:/Users/paula/Desktop/Geomatica/Cartografia_tematica/yuca_2018.png")

NA

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 
[2] LC_CTYPE=Spanish_Colombia.1252   
[3] LC_MONETARY=Spanish_Colombia.1252
[4] LC_NUMERIC=C                     
[5] LC_TIME=Spanish_Colombia.1252    

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

other attached packages:
 [1] knitr_1.28            SpatialPosition_2.0.1
 [3] cartography_2.4.1     sf_0.9-2             
 [5] rgeos_0.5-2           readxl_1.3.1         
 [7] forcats_0.5.0         stringr_1.4.0        
 [9] dplyr_0.8.5           purrr_0.3.4          
[11] readr_1.3.1           tidyr_1.0.2          
[13] tibble_3.0.0          ggplot2_3.3.0        
[15] tidyverse_1.3.0       raster_3.1-5         
[17] sp_1.4-1             

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4.6       lubridate_1.7.8    lattice_0.20-38   
 [4] png_0.1-7          class_7.3-15       assertthat_0.2.1  
 [7] digest_0.6.25      R6_2.4.1           cellranger_1.1.0  
[10] backports_1.1.6    reprex_0.3.0       evaluate_0.14     
[13] e1071_1.7-3        httr_1.4.1         pillar_1.4.3      
[16] rlang_0.4.5        rstudioapi_0.11    rmarkdown_2.1     
[19] rgdal_1.4-8        munsell_0.5.0      broom_0.5.6       
[22] compiler_3.6.3     modelr_0.1.6       xfun_0.13         
[25] pkgconfig_2.0.3    base64enc_0.1-3    htmltools_0.4.0   
[28] tidyselect_1.0.0   codetools_0.2-16   fansi_0.4.1       
[31] crayon_1.3.4       dbplyr_1.4.3       withr_2.2.0       
[34] grid_3.6.3         lwgeom_0.2-3       nlme_3.1-144      
[37] jsonlite_1.6.1     gtable_0.3.0       lifecycle_0.2.0   
[40] DBI_1.1.0          magrittr_1.5       units_0.6-6       
[43] scales_1.1.0       KernSmooth_2.23-16 cli_2.0.2         
[46] stringi_1.4.6      fs_1.4.1           xml2_1.3.1        
[49] ellipsis_0.3.0     generics_0.0.2     vctrs_0.2.4       
[52] tools_3.6.3        glue_1.4.0         hms_0.5.3         
[55] slippymath_0.3.1   yaml_2.2.1         colorspace_1.4-1  
[58] isoband_0.2.1      classInt_0.4-3     rvest_0.3.5       
[61] haven_2.2.0