Estadísticas municipales de la agricultura en el amazonas (2007-2018)

1. Introducción

El proposito de este Notebook es brindar información acerca de los municipios del amazonas con mejor rendimiento en los cultivos de mayor producción entre los años 2007 al 2018. Para realizarlo se necesitó de la unión de datos no espaciales con datos espaciales.

Cabe resaltar que para llevar a cabo la unión de los datos no fue necesario hacer Data Wrangling ya que los valores coincidian; unicamente fue necesario cambiar el tipo de datos.

2. Funcionalidades del SIG

Para iniciar, el primer paso para antes de cargar los datos y procesarlos fue eliminar el contenido de la memoria para no tener inconvenientes con los paquetes y datos.

rm(list = ls())

Después, instalar los paquetes que aun no están instalados:

list.of.packages <- c("here", "tidyverse", "rgeos", "maptools", "raster", "sf",  "viridis", "rnaturalearth", "GSODR", "ggrepel", "cowplot")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

Cargar las librerias:

library(here)

here() starts at G:/Geomatica/Evaluaciones_Agropecuarias_Municipales_EVA

library(tidyverse)

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.3
v tibble  2.1.3     v stringr 1.4.0
v tidyr   1.0.2     v forcats 0.5.0
v readr   1.3.1     
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

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(maptools)

Checking rgeos availability: TRUE

library(raster)

Attaching package: 㤼㸱raster㤼㸲

The following object is masked from 㤼㸱package:tidyr㤼㸲:

    extract

The following object is masked from 㤼㸱package:dplyr㤼㸲:

    select

library(sf)

Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

library(viridis)

Loading required package: viridisLite

library(rnaturalearth)
library(GSODR)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

library(ggrepel)
library(cowplot)

********************************************************
Note: As of version 1.0.0, cowplot does not change the
  default ggplot2 theme anymore. To recover the previous
  behavior, execute:
  theme_set(theme_cowplot())
********************************************************

3. Explorar estadísticas agricolas en el amazonas

Leer el archivo csv que ya fue editado para que solo tenga los datos del departamento del amazonas. Este archivo fue descargado de los datos de la página del Gobierno de Colombia.

datos <- read_csv("G:/Geomatica/Evaluaciones_Agropecuarias_Municipales_EVA/EVA_Amazonas.csv")

Parsed with column specification:
cols(
  CODIGO_DEP = col_double(),
  DEPARTAMENTO = col_character(),
  CODIGO_MUN = col_double(),
  MUNICIPIO = col_character(),
  GRUPO = col_character(),
  SUBGRUPO = col_character(),
  CULTIVO = col_character(),
  YEAR = col_double(),
  AREA_SEMBRADA = col_double(),
  AREA_COSECHADA = col_double(),
  PRODUCCION = col_double(),
  RENDIMIENTO = col_double(),
  ESTADO = col_character(),
  CICLO = col_character()
)

Ver cuales son los atributos de datos

head(datos)

tail(datos)

Con el uso de la libreria dplyr se muestran unicamente los datos promediados de los grupos y municipios a lo largo de los años.

datos %>%
  group_by(MUNICIPIO, GRUPO) %>%
  summarise(rend_prom = mean(RENDIMIENTO, na.rm = TRUE)) -> rend_resumen
head(rend_resumen)

Y también solo el rendimiento por grupo en los municipios:

datos %>%
  group_by(GRUPO) %>%
  summarise(rend_dep = mean(RENDIMIENTO, na.rm = TRUE)) -> rend_amazonas

rend_amazonas

Observando los datos anteriores se puede ver que los cultivos con mayor rendimiento son FRUTALES y TUBERCULOS Y PLATANOS.

Separar por mayor rendimiento los grupos de cultivo en el año 2018:

datos %>% 
  filter(YEAR==2018) %>% 
  group_by(GRUPO, MUNICIPIO) %>%
  summarize(max_rend = max(RENDIMIENTO, na.rm = TRUE)) %>%
    slice(which.max(max_rend)) -> rend_max_18

rend_max_18

Separar los municipios con mayor area cosechada con respecto al grupo de cultivos en el año 2018:

datos %>% 
  filter(YEAR==2018) %>% 
  group_by(GRUPO, MUNICIPIO) %>%
  summarize(max_area_cosecha = max(AREA_COSECHADA, na.rm = TRUE)) %>%
    slice(which.max(max_area_cosecha)) -> area_cosecha_max

area_cosecha_max

Como se pudo observar, en este caso el area mas grande cosechada en 2018 fue el grupo de TUBERCULOS Y PLATANOS. El cultivo que más se produce en Puerto Nariño es la Yuca.

Seleccionar la producción de yuca en Puerto Nariño para cada uno de los años:

datos %>% 
  filter(MUNICIPIO=="PUERTO NARIÑO" & SUBGRUPO=="YUCA") %>% 
  group_by(YEAR, CULTIVO) ->  puerto_yuca

puerto_yuca

Ahora, graficar los datos obtenidos para tener una mejor visualizacion de estas variables y poder realizar comparaciones:

g <- ggplot(aes(x=YEAR, y=PRODUCCION/1000), data = puerto_yuca) + geom_bar(stat='identity') + labs(y='Produccion de yuca (Ton x 1000)')

g + ggtitle("Evolucion de produccion de yuca en Puerto Nariño de 2007 a 2018") + labs(caption= "Based on EMA data (DANE, 2018)")

Para saber que cultivos tuvieron la mayor area cosechada en 2018, ejecutamos los siguientes comandos:

datos %>% 
  filter(YEAR==2018) %>% 
  group_by(GRUPO) %>%
  summarize(sum_area_cosecha = sum(AREA_COSECHADA, na.rm = TRUE)) %>%
     arrange(desc(sum_area_cosecha)) -> total_area_cosecha

total_area_cosecha

Como se puede observar, los Tuberculos y platanos tuvieron la mayor area cosechada en 2018 en el Amazonas. Ahora, veremos cuales son los cultivos pertenecientes a este grupo.

datos %>%
  filter(GRUPO=="TUBERCULOS Y PLATANOS" & YEAR==2018) %>%
  group_by(CULTIVO) %>%
  summarize(sum_cosecha = sum(AREA_COSECHADA, na.rm = TRUE)) %>%
     arrange(desc(sum_cosecha)) -> total_cosecha


total_cosecha

Como se ha podido observar, el amazonas es un departamento donde las zonas con clutivo fueron bajas antes de los años 2015-2016 aproximadamente. Segun los datos observados, desde el año 2015 se ha visto un incremento significativo en el aumento de producción de yuca, convirtiendo este cultivo en el de mayor producción en el año 2018 en el amazonas. Del mismo modo, en el año 2019 hubo un aumento de los cultivos de yuca gracias a aportes tanto de la comunidad como del Instituto Colombiano Agropecuario.

Ahora, se seleccionan los municipios que tuvieron la mayor area cosechada para cada cultivo de tuberculos y platanos en 2018:

datos %>% 
  filter(YEAR==2018 & GRUPO=="TUBERCULOS Y PLATANOS") %>% 
  group_by(CULTIVO, MUNICIPIO) %>%
  summarize(max_area2 = max(AREA_COSECHADA, na.rm = TRUE)) %>%
    slice(which.max(max_area2)) -> area_cosecha2

area_cosecha2

Agregamos un campo con abreviaturas para cada grupo de cultivos para que todo se vea en orden.

total_area_cosecha$CROP <- abbreviate(total_area_cosecha$GRUPO, 3)

Plotear:

g <- ggplot(aes(x=CROP, y=sum_area_cosecha), data = total_area_cosecha) + geom_bar(stat='identity') + labs(y='Area total cosechada (Ha)')

g+ ggtitle("Area total cosechada por grupos de cultivos en 2018 Amazonas") + theme(plot.title = element_text(hjust = 0.5)) +
   labs(caption= "Based on EMA data (DANE, 2018)")

4. Unir estadisticas agricolas a las municipales

Descargar los datos administrativos del Amazonas disponibles en el geoportal del DANE.

Abrir la libreria sf y despues leer los datos administrativos descargados:

amz_munic <- sf::st_read("G:/Geomatica/Datos_de _elevacion_en_R/Amazonas_datos_DANE/91_AMAZONAS/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Reading layer `MGN_MPIO_POLITICO' from data source `G:\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

Abrir los datos de la tabla amz_munic y verificar que tipos de datos y variables hay:

amz_munic

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
First 10 features:
   DPTO_CCDGO MPIO_CCDGO                MPIO_CNMBR
1          91      91001                   LETICIA
2          91      91263                EL ENCANTO
3          91      91405               LA CHORRERA
4          91      91407                LA PEDRERA
5          91      91430       LA VICTORIA (Pacoa)
6          91      91460 MIRITÍ-PARANÁ (Campoamor)
7          91      91530            PUERTO ALEGRÍA
8          91      91536              PUERTO ARICA
9          91      91540             PUERTO NARIÑO
10         91      91669    SANTANDER (Araracuara)
                        MPIO_CRSLC MPIO_NAREA MPIO_NANO
1    Decreto 352 de Feb 20 de 1964   6183.129      2017
2   Decreto 274 de Mayo 28 de 1953  10897.481      2017
3   Decreto 274 de Mayo 28 de 1953  12726.929      2017
4   Decreto 274 de Mayo 28 de 1953  13670.317      2017
5        ORD 12 DE JULIO 9 DE 1996   1432.795      2017
6   Decreto 274 de Mayo 28 de 1953  16866.120      2017
7      ORD 24 DE AGOSTO 01 DE 1997   8494.128      2017
8      ORD 23 DE AGOSTO 01 DE 1997  13791.174      2017
9  Decreto 106 de Enero 18 de 1984   1509.167      2017
10  Decreto 274 de Mayo 28 de 1953  14731.085      2017
   DPTO_CNMBR Shape_Leng Shape_Area
1    AMAZONAS   3.880001  0.5007018
2    AMAZONAS   8.052697  0.8853254
3    AMAZONAS   9.105101  1.0335830
4    AMAZONAS   7.376291  1.1051976
5    AMAZONAS   2.483170  0.1160863
6    AMAZONAS   6.448395  1.3668663
7    AMAZONAS   4.777815  0.6901332
8    AMAZONAS   8.677672  1.1180316
9    AMAZONAS   1.968859  0.1223627
10   AMAZONAS  10.636046  1.1952334
                         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...
7  POLYGON ((-73.78757 -0.4681...
8  POLYGON ((-71.18864 -1.3893...
9  POLYGON ((-70.36439 -3.6903...
10 POLYGON ((-71.99609 -0.2763...

Ahora, se elige una variable comun, en este caso se elige MPIO_CCDGO de amz_munic ya que junto con la tabla de datos los valores coinciden con los de CODIGO_MUN.

Verificar los datos y tipo de datos que hay en un solo municipio, en este caso se elige Leticia:

datos %>% filter (MUNICIPIO=="LETICIA") -> ltc_datos

ltc_datos

class(ltc_datos$CODIGO_MUN)

[1] "numeric"

Ahora, adicionamos una nueva columna a la tabla datos llamada MPIO_CCDGO que tendrá los mismos valores que la columna CODIGO_MUN, esto para evitar confusiones al momento de realizar la unión de los datos con la tabla amz_munic que tiene una columna con el mismo nombre (MPIO_CCDGO).

En este paso no se realizo la edición de los valores de la nueva columna MPIO_CCDGO ya que los valores coinciden con los de la tabla amz_munic, no hay necesidad de editar los datos ni realizar copia de la tabla datos.

Adicionalmente, se cambia el tipo de datos de esta columna para que más adelante cuando se lleve a cabo el join el tipo de datos y los valores sean iguales.

datos$MPIO_CCDGO <- as.character(datos$CODIGO_MUN)

class(datos$MPIO_CCDGO)

[1] "character"

Para asegurar que efectivamente se adiciono una nueva columna con los datos requeridos, ejecutamos el siguiente codigo que muestra las primeras seis filas de la tabla:

head(datos)

El siguiente paso es filtrar los datos en una nueva tabla (datos2) de un solo cultivo, en este caso la yuca:

datos %>% filter(CULTIVO == "YUCA") -> datos2

head(datos2)

class(datos2)

[1] "spec_tbl_df" "tbl_df"      "tbl"         "data.frame" 

Nuevamente filtramos, pero en este caso solo se muestran los atributos de MUNICIPIO, MPIO_CCDGO, YEAR, PRODUCCION y RENDIMIENTO en la nueva tabla datos3, especificamente para el cultivo de yuca.

datos3 <- datos2 %>% dplyr::select(MUNICIPIO, MPIO_CCDGO, YEAR, PRODUCCION, RENDIMIENTO)

datos3

datos3  %>%
  gather("YEAR", "PRODUCCION", "RENDIMIENTO", key = variable, value = number)

head(datos3)

Lo siguiente corresponde a convertir la tabla de atributos de formato largo a formato ancho:

datos3 %>%
  group_by(MPIO_CCDGO) %>%
  mutate(Visit = 1:n()) %>%
  gather("YEAR", "PRODUCCION", "RENDIMIENTO", key = variable, value = number) %>%
  unite(combi, variable, Visit) %>%
  spread(combi, number) -> datos4

head(datos4)

tail(datos4)

Realizamos una copia de la tabla amz_munic en el caso de que se presente algun error ejecutando los comandos.

amz_munic2 <- amz_munic

En este paso lo que se hace es cambiar el tipo de datos de la tabla amz_munic2 de factor a caracter, para que el join entre las tablas amz_munic2 y datos4 sea exitoso.

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

Realizar el join.

amz_munic_stat = left_join(amz_munic2, datos4, by="MPIO_CCDGO")

Mostrar el resumen de la unión realizada y verificar que los datos sean correctos.

summary(amz_munic_stat)

 DPTO_CCDGO  MPIO_CCDGO                            MPIO_CNMBR                           MPIO_CRSLC   MPIO_NAREA      MPIO_NANO   
 91:11      Length:11          EL ENCANTO               :1    Decreto 106 de Enero 18 de 1984:1    Min.   : 1433   Min.   :2017  
            Class :character   LA CHORRERA              :1    Decreto 274 de Mayo 28 de 1953 :6    1st Qu.: 7339   1st Qu.:2017  
            Mode  :character   LA PEDRERA               :1    Decreto 352 de Feb 20 de 1964  :1    Median :10897   Median :2017  
                               LA VICTORIA (Pacoa)      :1    ORD 12 DE JULIO 9 DE 1996      :1    Mean   : 9954   Mean   :2017  
                               LETICIA                  :1    ORD 23 DE AGOSTO 01 DE 1997    :1    3rd Qu.:13731   3rd Qu.:2017  
                               MIRITÍ-PARANÁ (Campoamor):1    ORD 24 DE AGOSTO 01 DE 1997    :1    Max.   :16866   Max.   :2017  
                               (Other)                  :5                                                                       
    DPTO_CNMBR   Shape_Leng       Shape_Area      MUNICIPIO          PRODUCCION_1    PRODUCCION_10  PRODUCCION_11  PRODUCCION_12 
 AMAZONAS:11   Min.   : 1.969   Min.   :0.1161   Length:11          Min.   : 345.0   Min.   :2195   Min.   :4332   Min.   :7500  
               1st Qu.: 4.329   1st Qu.:0.5954   Class :character   1st Qu.: 555.0   1st Qu.:2926   1st Qu.:4808   1st Qu.:7625  
               Median : 6.448   Median :0.8853   Mode  :character   Median : 606.0   Median :3658   Median :5284   Median :7750  
               Mean   : 6.238   Mean   :0.8070                      Mean   : 886.2   Mean   :3658   Mean   :5284   Mean   :7750  
               3rd Qu.: 8.365   3rd Qu.:1.1116                      3rd Qu.:1450.0   3rd Qu.:4389   3rd Qu.:5759   3rd Qu.:7875  
               Max.   :10.636   Max.   :1.3669                      Max.   :1475.0   Max.   :5120   Max.   :6235   Max.   :8000  
                                                                    NA's   :6        NA's   :9      NA's   :9      NA's   :9     
  PRODUCCION_2   PRODUCCION_3     PRODUCCION_4  PRODUCCION_5    PRODUCCION_6    PRODUCCION_7    PRODUCCION_8    PRODUCCION_9  
 Min.   : 335   Min.   : 414.0   Min.   :480   Min.   :740.0   Min.   :200.0   Min.   :494.0   Min.   :220.0   Min.   :234.0  
 1st Qu.: 584   1st Qu.: 578.5   1st Qu.:560   1st Qu.:742.5   1st Qu.:202.5   1st Qu.:498.8   1st Qu.:227.5   1st Qu.:304.5  
 Median : 636   Median : 743.0   Median :640   Median :745.0   Median :205.0   Median :503.5   Median :235.0   Median :375.0  
 Mean   :1628   Mean   : 743.0   Mean   :640   Mean   :745.0   Mean   :205.0   Mean   :503.5   Mean   :235.0   Mean   :375.0  
 3rd Qu.:3185   3rd Qu.: 907.5   3rd Qu.:720   3rd Qu.:747.5   3rd Qu.:207.5   3rd Qu.:508.2   3rd Qu.:242.5   3rd Qu.:445.5  
 Max.   :3401   Max.   :1072.0   Max.   :800   Max.   :750.0   Max.   :210.0   Max.   :513.0   Max.   :250.0   Max.   :516.0  
 NA's   :6      NA's   :9        NA's   :9     NA's   :9       NA's   :9       NA's   :9       NA's   :9       NA's   :9      
 RENDIMIENTO_1   RENDIMIENTO_10   RENDIMIENTO_11  RENDIMIENTO_12  RENDIMIENTO_2   RENDIMIENTO_3   RENDIMIENTO_4 RENDIMIENTO_5
 Min.   :2.500   Min.   : 3.530   Min.   :6.540   Min.   :6.670   Min.   :3.400   Min.   :2.000   Min.   :2     Min.   :2    
 1st Qu.:2.830   1st Qu.: 5.147   1st Qu.:6.603   1st Qu.:6.707   1st Qu.:3.560   1st Qu.:2.572   1st Qu.:2     1st Qu.:2    
 Median :3.300   Median : 6.765   Median :6.665   Median :6.745   Median :3.770   Median :3.145   Median :2     Median :2    
 Mean   :3.198   Mean   : 6.765   Mean   :6.665   Mean   :6.745   Mean   :3.662   Mean   :3.145   Mean   :2     Mean   :2    
 3rd Qu.:3.520   3rd Qu.: 8.383   3rd Qu.:6.728   3rd Qu.:6.782   3rd Qu.:3.770   3rd Qu.:3.717   3rd Qu.:2     3rd Qu.:2    
 Max.   :3.840   Max.   :10.000   Max.   :6.790   Max.   :6.820   Max.   :3.810   Max.   :4.290   Max.   :2     Max.   :2    
 NA's   :6       NA's   :9        NA's   :9       NA's   :9       NA's   :6       NA's   :9       NA's   :9     NA's   :9    
 RENDIMIENTO_6 RENDIMIENTO_7 RENDIMIENTO_8 RENDIMIENTO_9     YEAR_1        YEAR_10        YEAR_11        YEAR_12         YEAR_2    
 Min.   :2     Min.   :3.8   Min.   :2     Min.   :2     Min.   :2007   Min.   :2016   Min.   :2017   Min.   :2018   Min.   :2008  
 1st Qu.:2     1st Qu.:3.8   1st Qu.:2     1st Qu.:2     1st Qu.:2007   1st Qu.:2016   1st Qu.:2017   1st Qu.:2018   1st Qu.:2008  
 Median :2     Median :3.8   Median :2     Median :2     Median :2007   Median :2016   Median :2017   Median :2018   Median :2008  
 Mean   :2     Mean   :3.8   Mean   :2     Mean   :2     Mean   :2007   Mean   :2016   Mean   :2017   Mean   :2018   Mean   :2008  
 3rd Qu.:2     3rd Qu.:3.8   3rd Qu.:2     3rd Qu.:2     3rd Qu.:2007   3rd Qu.:2016   3rd Qu.:2017   3rd Qu.:2018   3rd Qu.:2008  
 Max.   :2     Max.   :3.8   Max.   :2     Max.   :2     Max.   :2007   Max.   :2016   Max.   :2017   Max.   :2018   Max.   :2008  
 NA's   :9     NA's   :9     NA's   :9     NA's   :9     NA's   :6      NA's   :9      NA's   :9      NA's   :9      NA's   :6     
     YEAR_3         YEAR_4         YEAR_5         YEAR_6         YEAR_7         YEAR_8         YEAR_9              geometry 
 Min.   :2009   Min.   :2010   Min.   :2011   Min.   :2012   Min.   :2013   Min.   :2014   Min.   :2015   POLYGON      :11  
 1st Qu.:2009   1st Qu.:2010   1st Qu.:2011   1st Qu.:2012   1st Qu.:2013   1st Qu.:2014   1st Qu.:2015   epsg:4326    : 0  
 Median :2009   Median :2010   Median :2011   Median :2012   Median :2013   Median :2014   Median :2015   +proj=long...: 0  
 Mean   :2009   Mean   :2010   Mean   :2011   Mean   :2012   Mean   :2013   Mean   :2014   Mean   :2015                     
 3rd Qu.:2009   3rd Qu.:2010   3rd Qu.:2011   3rd Qu.:2012   3rd Qu.:2013   3rd Qu.:2014   3rd Qu.:2015                     
 Max.   :2009   Max.   :2010   Max.   :2011   Max.   :2012   Max.   :2013   Max.   :2014   Max.   :2015                     
 NA's   :9      NA's   :9      NA's   :9      NA's   :9      NA's   :9      NA's   :9      NA's   :9                        

5. Plotting

En este paso se requiere instalar el paquete “RColorBrewer” para llevar a cabo algunas de las caracteristicas del mapa.

Abrir las librerias:

library(RColorBrewer)
library(leaflet)

Por ultimo, plotear los municipios con su correspondiente producción de yuca específicamente en el año 2018.

En este caso fue modificada la escala de valores y los colores usados para representar dicha escala. También se cambió el color y grosor de la silueta de cada municipio. Observe que se seleccionó la columna PRODUCCION_12 de la nueva tabla amz_munic_stat ya que corresponde a los datos de producción de yuca del año 2018 en Toneladas.

bins <- c(0,200, 400, 800, 1600, 3200, 6400, 9000)
pal <- colorBin("viridis", domain = amz_munic_stat$PRODUCCION_12, bins = bins)

mapa <- leaflet(data = amz_munic_stat) %>%
  addTiles()  %>%
  addPolygons(label = ~PRODUCCION_12,
              popup = ~MPIO_CCDGO,
              fillColor = ~pal(PRODUCCION_12),
              color = "#444444",
              weight = 1,
              smoothFactor = 0.5,
              opacity = 2,
              fillOpacity = 0.5,
              highlightOptions = highlightOptions(color = "black", weight = 3, bringToFront = TRUE)
              ) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend("bottomright", pal = pal, values = ~PRODUCCION_12,
            title = "Produccion de yuca en el Amazonas [Ton] (2018)",
            opacity = 1
            )

mapa

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] devtools_2.3.0      usethis_1.6.0       leaflet_2.0.3      
 [4] RColorBrewer_1.1-2  cowplot_1.0.0       ggrepel_0.8.2      
 [7] GSODR_2.0.1         rnaturalearth_0.1.0 viridis_0.5.1      
[10] viridisLite_0.3.0   sf_0.9-0            raster_3.0-12      
[13] maptools_0.9-9      rgeos_0.5-2         sp_1.4-1           
[16] forcats_0.5.0       stringr_1.4.0       purrr_0.3.3        
[19] readr_1.3.1         tidyr_1.0.2         tibble_2.1.3       
[22] ggplot2_3.3.0       tidyverse_1.3.0     here_0.1           
[25] dplyr_0.8.5        

loaded via a namespace (and not attached):
 [1] nlme_3.1-144            fs_1.3.2               
 [3] lubridate_1.7.4         httr_1.4.1             
 [5] rprojroot_1.3-2         tools_3.6.3            
 [7] backports_1.1.5         R6_2.4.1               
 [9] KernSmooth_2.23-16      DBI_1.1.0              
[11] colorspace_1.4-1        withr_2.1.2            
[13] prettyunits_1.1.1       processx_3.4.2         
[15] tidyselect_1.0.0        gridExtra_2.3          
[17] compiler_3.6.3          cli_2.0.2              
[19] rvest_0.3.5             xml2_1.2.5             
[21] desc_1.2.0              scales_1.1.0           
[23] classInt_0.4-2          callr_3.4.3            
[25] digest_0.6.25           foreign_0.8-75         
[27] rmarkdown_2.1           base64enc_0.1-3        
[29] pkgconfig_2.0.3         htmltools_0.4.0        
[31] sessioninfo_1.1.1       dbplyr_1.4.2           
[33] htmlwidgets_1.5.1       rlang_0.4.5            
[35] readxl_1.3.1            rstudioapi_0.11        
[37] generics_0.0.2          farver_2.0.3           
[39] jsonlite_1.6.1          crosstalk_1.1.0.1      
[41] magrittr_1.5            Rcpp_1.0.3             
[43] munsell_0.5.0           fansi_0.4.1            
[45] lifecycle_0.2.0         stringi_1.4.6          
[47] yaml_2.2.1              pkgbuild_1.0.6         
[49] grid_3.6.3              crayon_1.3.4           
[51] lattice_0.20-38         haven_2.2.0            
[53] hms_0.5.3               ps_1.3.2               
[55] knitr_1.28              pillar_1.4.3           
[57] pkgload_1.0.2           codetools_0.2-16       
[59] reprex_0.3.0            glue_1.3.1             
[61] evaluate_0.14           leaflet.providers_1.9.0
[63] remotes_2.1.1           data.table_1.12.8      
[65] modelr_0.1.6            vctrs_0.2.3            
[67] testthat_2.3.2          cellranger_1.1.0       
[69] gtable_0.3.0            assertthat_0.2.1       
[71] xfun_0.12               broom_0.5.5            
[73] e1071_1.7-3             class_7.3-15           
[75] memoise_1.1.0           units_0.6-5            
[77] ellipsis_0.3.0