1. Introduccion

En este cuaderno de R Markdown se ilustra las estadísticas agricolas para el departamento de Cundinamarca en Colombia.

2. Funcionalidades SIG

La exploración de estadísticas no espaciales es esencial para comprender lo que está sucediendo en los territorios. Varias bibliotecas R, en particular dplyr y tidyverse, son muy útiles para explorar y resumir estadísticas.

Por otro lado, las operaciones geoespaciales pueden mejorar nuestra comprensión de varios problemas que afectan a las regiones geográficas. Por ejemplo, desea averiguar cuál es la ubicación de aquellos municipios cuyos rendimientos de cosecha son sobresalientes (o, alternativamente, más bajos que el promedio). Para realizar dicha exploración, necesitamos unir datos no espaciales con datos espaciales. Lo haremos con R.

Además, también podríamos explorar uniones espaciales. Estas operaciones se basan en la intersección entre dos objetos espaciales, a menudo puntos y los polígonos. Hay muchas formas en que podemos unir objetos, que pueden incluir opciones específicas como cruces, cerca, dentro, toques, etc. El punto es no presionarlo, ¡podemos hacerlo en otro R Notebook!

Primero eliminaremos el contenido de la memoria:

rm(list=ls())

Ahora, instalaremos las librerias que necesitamos. Tenga en cuenta que, en el siguiente fragmento, cualquier paquete se instala solo si no se ha instalado previamente.

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)
library(here)
library(tidyverse)
library(rgeos)
library(maptools)
library(raster)
library(sf)
library(viridis)
library(rnaturalearth)
library(GSODR)
library(ggrepel)
library(cowplot)

3. Explorando estadísticas agrícolas de Cundinamarca

Revisaremos el directorio de trabajo

getwd()
## [1] "C:/Users/mario/Documents/Ing. Agronomica/9 Semestre/Geomatica/Informes tecnicos/Estadisticas Agricolas"

Utilizando la funcion “read_csv2()” podremos leer los datos de excel que son formato “csv” separados por punto y coma:

datos <- read_csv2("./EVA_Cund.csv")
## Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
## Parsed with column specification:
## cols(
##   COD_DEP = col_double(),
##   DEPARTAMENTO = col_character(),
##   COD_MUN = col_double(),
##   MUNICIPIO = col_character(),
##   GRUPO = col_character(),
##   SUBGRUPO = col_character(),
##   CULTIVO = col_character(),
##   DESAGREGACION_REGIONAL_SISTEMA_PRODUCTIVO = col_character(),
##   YEAR = col_double(),
##   PERIODO = col_character(),
##   Area_Siembra = col_double(),
##   Area_Cosecha = col_double(),
##   Produccion = col_double(),
##   Rendimiento = col_double(),
##   ESTADO = col_character(),
##   N_CIENTIFICO = col_character(),
##   CICLO_CULTIVO = col_character()
## )

Ahora revisaremos que los datos esten bien, mirando los primeros datos “head()” y ultimos “tail()” de la base de datos del departamento de Cundinamarca.

head(datos)
## # A tibble: 6 x 17
##   COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO SUBGRUPO CULTIVO DESAGREGACION_R~
##     <dbl> <chr>          <dbl> <chr>     <chr> <chr>    <chr>   <chr>           
## 1      25 CUNDINAMARCA   25754 SOACHA    HORT~ ACELGA   ACELGA  ACELGA          
## 2      25 CUNDINAMARCA   25214 COTA      HORT~ ACELGA   ACELGA  ACELGA          
## 3      25 CUNDINAMARCA   25754 SOACHA    HORT~ ACELGA   ACELGA  ACELGA          
## 4      25 CUNDINAMARCA   25214 COTA      HORT~ ACELGA   ACELGA  ACELGA          
## 5      25 CUNDINAMARCA   25754 SOACHA    HORT~ ACELGA   ACELGA  ACELGA          
## 6      25 CUNDINAMARCA   25126 CAJICA    HORT~ ACELGA   ACELGA  ACELGA          
## # ... with 9 more variables: YEAR <dbl>, PERIODO <chr>, Area_Siembra <dbl>,
## #   Area_Cosecha <dbl>, Produccion <dbl>, Rendimiento <dbl>, ESTADO <chr>,
## #   N_CIENTIFICO <chr>, CICLO_CULTIVO <chr>
tail(datos)
## # A tibble: 6 x 17
##   COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO SUBGRUPO CULTIVO DESAGREGACION_R~
##     <dbl> <chr>          <dbl> <chr>     <chr> <chr>    <chr>   <chr>           
## 1      25 CUNDINAMARCA   25524 PANDI     HORT~ CALABAC~ CALABA~ ZUCCHINI        
## 2      25 CUNDINAMARCA   25436 MANTA     HORT~ CALABAC~ CALABA~ ZUCCHINI        
## 3      25 CUNDINAMARCA   25524 PANDI     HORT~ CALABAC~ CALABA~ ZUCCHINI        
## 4      25 CUNDINAMARCA   25436 MANTA     HORT~ CALABAC~ CALABA~ ZUCCHINI        
## 5      25 CUNDINAMARCA   25807 TIBIRITA  HORT~ CALABAC~ CALABA~ ZUCCHINI        
## 6      25 CUNDINAMARCA   25524 PANDI     HORT~ CALABAC~ CALABA~ ZUCCHINI        
## # ... with 9 more variables: YEAR <dbl>, PERIODO <chr>, Area_Siembra <dbl>,
## #   Area_Cosecha <dbl>, Produccion <dbl>, Rendimiento <dbl>, ESTADO <chr>,
## #   N_CIENTIFICO <chr>, CICLO_CULTIVO <chr>

Estos datos estadisticos nos brindan informacion sobre el área sembrada, área cosechada, y rendimiento para diferentes cultivos en diferentes años. Los atributos SUBGRUPO y CULTIVO parecen referirse a la misma cosa (es decir, un cultivo). Los cultivos se clasifican además en un GRUPO dado.

En esta tabla, no tenemos unidades. Sin embargo, si verificamos el archivo csv original, encontramos que las unidades de área son (Hectáreas) y que las unidades de rendimiento son (Ton/ha)

Utilizaremos la biblioteca dplyr para explorar el contenido del objeto de datos.

Primero, obtengamos un resumen del rendimiento (es decir, el rendimiento promedio durante varios años) por grupo y municipio:

datos %>%
  group_by(MUNICIPIO, GRUPO) %>%
  summarise(rend_prom = mean(Rendimiento, na.rm = TRUE)) -> rend_resumen
### Let's visualize only the first six records
head(rend_resumen)
## # A tibble: 6 x 3
## # Groups:   MUNICIPIO [1]
##   MUNICIPIO    GRUPO                                            rend_prom
##   <chr>        <chr>                                                <dbl>
## 1 AGUA DE DIOS CEREALES                                              44.3
## 2 AGUA DE DIOS FIBRAS                                               154. 
## 3 AGUA DE DIOS FRUTALES                                             314. 
## 4 AGUA DE DIOS OTROS PERMANENTES                                     11.4
## 5 AGUA DE DIOS PLANTAS AROMATICAS, CONDIMENTARIAS Y MEDICINALES      28.9
## 6 AGUA DE DIOS TUBERCULOS Y PLATANOS                                109.
tail(rend_resumen)
## # A tibble: 6 x 3
## # Groups:   MUNICIPIO [1]
##   MUNICIPIO GRUPO                                            rend_prom
##   <chr>     <chr>                                                <dbl>
## 1 ZIPAQUIRA CEREALES                                              42.2
## 2 ZIPAQUIRA FRUTALES                                              20.2
## 3 ZIPAQUIRA HORTALIZAS                                            36.6
## 4 ZIPAQUIRA LEGUMINOSAS                                           65.3
## 5 ZIPAQUIRA PLANTAS AROMATICAS, CONDIMENTARIAS Y MEDICINALES      29  
## 6 ZIPAQUIRA TUBERCULOS Y PLATANOS                                204.

Tambien podemos calcular el rendimiento promedio por GRUPO en los municipios de Cundinamarca:

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

rend_cundinamarca
## # A tibble: 12 x 2
##    GRUPO                                            rend_dep
##    <chr>                                               <dbl>
##  1 CEREALES                                             77.5
##  2 FIBRAS                                              163. 
##  3 FLORES Y FOLLAJES                                   264. 
##  4 FORESTALES                                           19.5
##  5 FRUTALES                                            131. 
##  6 HONGOS                                                3.2
##  7 HORTALIZAS                                          166. 
##  8 LEGUMINOSAS                                          82.5
##  9 OLEAGINOSAS                                         223. 
## 10 OTROS PERMANENTES                                    63.7
## 11 PLANTAS AROMATICAS, CONDIMENTARIAS Y MEDICINALES     61.4
## 12 TUBERCULOS Y PLATANOS                               282.

Los rendimientos promedio más altos corresponden a FLORES Y FOLLAJES, OLEAGINOSAS, HORTALIZAS, FIBRAS Y FORESTALES.

Utilizando el paquete “dplyr” podremos ver cuáles son los municipios con mayor rendimiento para cada grupo de cultivos en 2018:

  1. “filter()” con esta funcion filtramos los datos de solamente el año 2018, que son los que queremos utilizar. 2.“group_by()” toma una tbl o data ya existente y lo convierte en un tbl o data agrupada donde las operaciones se realizan “por grupo”.
  2. “n.rm” nos indica si los valores faltantes se deben eliminar, en este caso al colocar “TRUE”, esta funcion eliminara si hay datos faltantes para evitar errores.
  3. “summarise()” Crea una o más variables escalares que resuman las variables de una tbl existente. Tbls con grupos creados por group_by () dará como resultado una fila en la salida para cada grupo. Tbls sin grupos dará como resultado una fila.
  4. “slice()” Elije filas por su posición ordinal en la tabla. Los tbls agrupados usan la posición ordinal dentro del grupo.
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
## # A tibble: 12 x 3
## # Groups:   GRUPO [12]
##    GRUPO                                            MUNICIPIO   max_rend
##    <chr>                                            <chr>          <dbl>
##  1 CEREALES                                         CHOCONTA        3033
##  2 FIBRAS                                           BELTRAN          296
##  3 FLORES Y FOLLAJES                                SUESCA          2528
##  4 FORESTALES                                       PARATEBUENO       15
##  5 FRUTALES                                         VIOTA           4101
##  6 HONGOS                                           GRANADA            4
##  7 HORTALIZAS                                       FOMEQUE        10084
##  8 LEGUMINOSAS                                      LA PALMA         305
##  9 OLEAGINOSAS                                      MEDINA           357
## 10 OTROS PERMANENTES                                SASAIMA          571
## 11 PLANTAS AROMATICAS, CONDIMENTARIAS Y MEDICINALES CHIPAQUE         456
## 12 TUBERCULOS Y PLATANOS                            GUATAVITA       3265

Ahora, veamos cuáles son los municipios con mayor área cosechada para cada grupo de cultivos en 2018:

datos %>% 
  filter(YEAR==2018) %>% 
  group_by(GRUPO, MUNICIPIO) %>%
  summarize(max_area_cosecha = max(Area_Cosecha, na.rm = TRUE)) %>%
    slice(which.max(max_area_cosecha)) -> area_cosecha_max
area_cosecha_max
## # A tibble: 12 x 3
## # Groups:   GRUPO [12]
##    GRUPO                                      MUNICIPIO         max_area_cosecha
##    <chr>                                      <chr>                        <dbl>
##  1 CEREALES                                   GACHALA                        765
##  2 FIBRAS                                     BELTRAN                        390
##  3 FLORES Y FOLLAJES                          NEMOCON                        400
##  4 FORESTALES                                 SAN JUAN DE RIOS~               91
##  5 FRUTALES                                   SAN JUAN DE RIOS~             2221
##  6 HONGOS                                     GRANADA                          1
##  7 HORTALIZAS                                 MOSQUERA                       600
##  8 LEGUMINOSAS                                UBAQUE                         856
##  9 OLEAGINOSAS                                PARATEBUENO                   5273
## 10 OTROS PERMANENTES                          CAPARRAPI                     9223
## 11 PLANTAS AROMATICAS, CONDIMENTARIAS Y MEDI~ CACHIPAY                       108
## 12 TUBERCULOS Y PLATANOS                      TAUSA                         8200

Tenga en cuenta que el área más grande cosechada en 2018 correspondió a CULTIVOS PERMANENTES (Caparrapi), para el cultivo de caña panelera.

Primero, seleccionemos la producción de OTROS CULTIVOS PERMANENTES (Toneladas) en Caparrapi para cada año:

datos %>% 
  filter(MUNICIPIO=="CAPARRAPI" & SUBGRUPO=="CAÑA") %>% 
  group_by(YEAR, CULTIVO) ->  caparrapi_cana
caparrapi_cana
## # A tibble: 12 x 17
## # Groups:   YEAR, CULTIVO [12]
##    COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO SUBGRUPO CULTIVO
##      <dbl> <chr>          <dbl> <chr>     <chr> <chr>    <chr>  
##  1      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
##  2      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
##  3      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
##  4      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
##  5      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
##  6      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
##  7      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
##  8      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
##  9      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
## 10      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
## 11      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
## 12      25 CUNDINAMARCA   25148 CAPARRAPI OTRO~ CAÑA     CAÑA P~
## # ... with 10 more variables: DESAGREGACION_REGIONAL_SISTEMA_PRODUCTIVO <chr>,
## #   YEAR <dbl>, PERIODO <chr>, Area_Siembra <dbl>, Area_Cosecha <dbl>,
## #   Produccion <dbl>, Rendimiento <dbl>, ESTADO <chr>, N_CIENTIFICO <chr>,
## #   CICLO_CULTIVO <chr>

Hagamos una exploración gráfica básica:

# Usaremos la libreria ggplo2
g <- ggplot(aes(x=YEAR, y=Produccion/1000), data = caparrapi_cana) + geom_bar(stat='identity') + labs(y='Produccion de Caña (Ton/Ha)')

Utilizando la funcion “geom_bar()” hace que la altura de la barra sea proporcional al número de casos en cada grupo.

g + ggtitle("Evolución de la producción de caña en Caparrapi de 2007 a 2018") + labs(caption= "Basado en EMA data (DANE, 2018)")

Ahora, investiguemos qué cultivos tuvieron la mayor área cosechada en 2018:

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

total_area_cosecha
## # A tibble: 12 x 2
##    GRUPO                                            sum_area_cosecha
##    <chr>                                                       <dbl>
##  1 OTROS PERMANENTES                                           75475
##  2 TUBERCULOS Y PLATANOS                                       56695
##  3 FRUTALES                                                    33575
##  4 CEREALES                                                    13686
##  5 HORTALIZAS                                                   8689
##  6 LEGUMINOSAS                                                  7124
##  7 OLEAGINOSAS                                                  5394
##  8 FLORES Y FOLLAJES                                            1817
##  9 PLANTAS AROMATICAS, CONDIMENTARIAS Y MEDICINALES              617
## 10 FIBRAS                                                        390
## 11 FORESTALES                                                    328
## 12 HONGOS                                                          1

Podemos ver que OTROS CULTIVOS PERMANENTES tuvieron la mayor parte del área cosechada en 2018 para Antioquia. ¿Qué cultivos pertenecen a este GRUPO? Puede buscar más información en DANE, la agencia gubernamental a cargo de las estadísticas nacionales en Colombia o también podemos buscar esta información en los datos mismos:

datos %>%
  filter(GRUPO=="OTROS PERMANENTES" & YEAR==2018) %>%
  group_by(CULTIVO) %>%
  summarize(sum_cosecha = sum(Area_Cosecha, na.rm = TRUE)) %>%
     arrange(desc(sum_cosecha)) -> total_cosecha


total_cosecha
## # A tibble: 3 x 2
##   CULTIVO       sum_cosecha
##   <chr>               <dbl>
## 1 CAÑA PANELERA       40481
## 2 CAFE                29087
## 3 CACAO                5907

Como se puede observar la CAÑA PANELERA es el cultivo con mayor area cosechada en el año 2018.

Ahora, veamos cuáles son los municipios con mayor área cosechada para cada cultivo permanente en 2018:

datos %>% 
  filter(YEAR==2018 & GRUPO=="OTROS PERMANENTES") %>% 
  group_by(CULTIVO, MUNICIPIO) %>%
  summarize(max_area2 = max(Area_Cosecha, na.rm = TRUE)) %>%
    slice(which.max(max_area2)) -> area_cosecha2

area_cosecha2
## # A tibble: 3 x 3
## # Groups:   CULTIVO [3]
##   CULTIVO       MUNICIPIO max_area2
##   <chr>         <chr>         <dbl>
## 1 CACAO         YACOPI         2371
## 2 CAFE          VIOTA          3179
## 3 CAÑA PANELERA CAPARRAPI      9223

Volvamos a los datos de cultivos permanentes. Antes de graficar, necesitaremos agregar, al objeto total_area_cosecha, un nuevo campo con abreviaturas para cada GRUPO de cultivos. De lo contrario, la grafica se verá desordenada.

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

Ahora, vamos a gráficar:

# Usaremos la libreria ggplot 2
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("Superficie total cosechada por grupos de cultivos en 2018 para Cundinamarca") + theme(plot.title = element_text(hjust = 0.5)) +
   labs(caption= "Based on EMA data (DANE, 2018)")

Ahora intentaremos obtener informacion sobre el rendimiento de diferentes cultivos del departamento de Cundinamarca.

  1. Filtraremos los datos para el año 2018, y traeremos los datos de rendimiento totales, por grupo de cultivo.
datos %>% 
  filter(YEAR==2018) %>% 
  group_by(GRUPO) %>%
  summarize(sum_rend = sum(Rendimiento, na.rm = TRUE)) %>%
     arrange(desc(sum_rend)) -> total_rend

total_rend
## # A tibble: 12 x 2
##    GRUPO                                            sum_rend
##    <chr>                                               <dbl>
##  1 HORTALIZAS                                          71166
##  2 TUBERCULOS Y PLATANOS                               50278
##  3 FRUTALES                                            22839
##  4 CEREALES                                            13955
##  5 LEGUMINOSAS                                         11285
##  6 OTROS PERMANENTES                                    9874
##  7 FLORES Y FOLLAJES                                    4658
##  8 OLEAGINOSAS                                          1029
##  9 PLANTAS AROMATICAS, CONDIMENTARIAS Y MEDICINALES      934
## 10 FIBRAS                                                296
## 11 FORESTALES                                             54
## 12 HONGOS                                                  4

Para realizar la grafica, pondremos abreviatura de 3 letras del total de datos de rendimiento de los grupos de cultivo para no visualizar la grafica desordenada

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

Ahora realizaremos la grafica:

# Usaremos la libreria ggplot 2
g <- ggplot(aes(x=CROP, y=sum_rend), data = total_rend) + geom_bar(stat='identity') + labs(y='Rendimiento [ton/Ha]', x='Grupo de cultivo')
g+ ggtitle("Rendimiento por grupos de cultivos en 2018 para Cundinamarca") + theme(plot.title = element_text(hjust = 0.5)) +
   labs(caption= "Based on EMA data (DANE, 2018)")

4. Uniendo estadísticas agrícolas a municipios

Ahora utilizaremos los datos administrativos de Cundinamarca, utilizando el archivo de forma correspondiente a Marco Geoestadistico Departamental que está disponible en el Geoportal del DANE.

Comencemos leyendo los datos usando la libreria sf:

ant_munic <- sf::st_read("C:/Users/mario/Documents/Ing. Agronomica/9 Semestre/Geomatica/Cundinamarca/DANE/25_CUNDINAMARCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\mario\Documents\Ing. Agronomica\9 Semestre\Geomatica\Cundinamarca\DANE\25_CUNDINAMARCA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 116 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -74.89063 ymin: 3.730129 xmax: -73.05256 ymax: 5.837258
## CRS:            4326
ant_munic
## Simple feature collection with 116 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -74.89063 ymin: 3.730129 xmax: -73.05256 ymax: 5.837258
## CRS:            4326
## First 10 features:
##    DPTO_CCDGO MPIO_CCDGO    MPIO_CNMBR                               MPIO_CRSLC
## 1          25      25483        NARIÑO                                     1899
## 2          25      25513         PACHO                                     1604
## 3          25      25506       VENECIA         Decreto 727 Septiembre 5 de 1951
## 4          25      25491       NOCAIMA                                     1735
## 5          25      25489       NIMAIMA            Ordenanza 30 de Julio de 1904
## 6          25      25488          NILO                                     1899
## 7          25      25486       NEMOCÓN                                     1778
## 8          25      25530   PARATEBUENO Ordenanza 36 del 30 de Noviembre de 1981
## 9          25      25535         PASCA                                     1887
## 10         25      25572 PUERTO SALGAR    Ordenanza 47 del 14 de Agosto de 1935
##    MPIO_NAREA MPIO_NANO   DPTO_CNMBR Shape_Leng  Shape_Area
## 1    55.16263      2017 CUNDINAMARCA  0.5219117 0.004493674
## 2   402.62678      2017 CUNDINAMARCA  1.2021630 0.032839624
## 3   122.20332      2017 CUNDINAMARCA  0.4873683 0.009951825
## 4    68.93823      2017 CUNDINAMARCA  0.3987992 0.005621814
## 5    59.49982      2017 CUNDINAMARCA  0.4990524 0.004852617
## 6   224.70968      2017 CUNDINAMARCA  0.8442993 0.018305852
## 7    98.36504      2017 CUNDINAMARCA  0.5509410 0.008022009
## 8   885.52351      2017 CUNDINAMARCA  1.7923439 0.072136135
## 9   271.71161      2017 CUNDINAMARCA  0.7143557 0.022134302
## 10  515.94145      2017 CUNDINAMARCA  1.3299992 0.042109370
##                          geometry
## 1  MULTIPOLYGON (((-74.7413 4....
## 2  MULTIPOLYGON (((-74.21183 5...
## 3  MULTIPOLYGON (((-74.44812 4...
## 4  MULTIPOLYGON (((-74.39011 5...
## 5  MULTIPOLYGON (((-74.34091 5...
## 6  MULTIPOLYGON (((-74.56068 4...
## 7  MULTIPOLYGON (((-73.88945 5...
## 8  MULTIPOLYGON (((-73.06503 4...
## 9  MULTIPOLYGON (((-74.24153 4...
## 10 MULTIPOLYGON (((-74.51278 5...

Tenga en cuenta que ant_munic es una colección de características simple y los datos utilizan el sistema de referencia de coordenadas geográficas WGS1984 (es decir, código 4326 epsg).

Podemos usar la función left_join para unir los municipios y las estadísticas agrícolas seleccionadas. Busquemos ayuda para comprender cómo funciona esta función, utilizando la funcion “help()” o escribiendo en la consola “?left_join”

help("left_join")
## starting httpd help server ... done
ant_munic
## Simple feature collection with 116 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -74.89063 ymin: 3.730129 xmax: -73.05256 ymax: 5.837258
## CRS:            4326
## First 10 features:
##    DPTO_CCDGO MPIO_CCDGO    MPIO_CNMBR                               MPIO_CRSLC
## 1          25      25483        NARIÑO                                     1899
## 2          25      25513         PACHO                                     1604
## 3          25      25506       VENECIA         Decreto 727 Septiembre 5 de 1951
## 4          25      25491       NOCAIMA                                     1735
## 5          25      25489       NIMAIMA            Ordenanza 30 de Julio de 1904
## 6          25      25488          NILO                                     1899
## 7          25      25486       NEMOCÓN                                     1778
## 8          25      25530   PARATEBUENO Ordenanza 36 del 30 de Noviembre de 1981
## 9          25      25535         PASCA                                     1887
## 10         25      25572 PUERTO SALGAR    Ordenanza 47 del 14 de Agosto de 1935
##    MPIO_NAREA MPIO_NANO   DPTO_CNMBR Shape_Leng  Shape_Area
## 1    55.16263      2017 CUNDINAMARCA  0.5219117 0.004493674
## 2   402.62678      2017 CUNDINAMARCA  1.2021630 0.032839624
## 3   122.20332      2017 CUNDINAMARCA  0.4873683 0.009951825
## 4    68.93823      2017 CUNDINAMARCA  0.3987992 0.005621814
## 5    59.49982      2017 CUNDINAMARCA  0.4990524 0.004852617
## 6   224.70968      2017 CUNDINAMARCA  0.8442993 0.018305852
## 7    98.36504      2017 CUNDINAMARCA  0.5509410 0.008022009
## 8   885.52351      2017 CUNDINAMARCA  1.7923439 0.072136135
## 9   271.71161      2017 CUNDINAMARCA  0.7143557 0.022134302
## 10  515.94145      2017 CUNDINAMARCA  1.3299992 0.042109370
##                          geometry
## 1  MULTIPOLYGON (((-74.7413 4....
## 2  MULTIPOLYGON (((-74.21183 5...
## 3  MULTIPOLYGON (((-74.44812 4...
## 4  MULTIPOLYGON (((-74.39011 5...
## 5  MULTIPOLYGON (((-74.34091 5...
## 6  MULTIPOLYGON (((-74.56068 4...
## 7  MULTIPOLYGON (((-73.88945 5...
## 8  MULTIPOLYGON (((-73.06503 4...
## 9  MULTIPOLYGON (((-74.24153 4...
## 10 MULTIPOLYGON (((-74.51278 5...

Necesitamos un atributo común (o variable compartida) para basar la unión. El mejor atributo es una identificación. En ant_munic, el atributo MPIO_CCDGO parece estar bien (se lee 25758 para SOPO). En datos, el atributo correspondiente es COD_MUN (lee 25758 para SOPO).

Podemos correr el siguiente codigo, filtrando unicamente el municipio que queremos ver para conocer el correspondiente COD_MUN

ant_munic %>% 
  filter(MPIO_CNMBR=="SOPO") %>%
  group_by(MPIO_CCDGO, MPIO_CNMBR ) ->  dats2

dats2
## Simple feature collection with 1 feature and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -74.01106 ymin: 4.812494 xmax: -73.91581 ymax: 4.974034
## CRS:            4326
## # A tibble: 1 x 10
## # Groups:   MPIO_CCDGO, MPIO_CNMBR [1]
##   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR
## * <fct>      <fct>      <fct>      <fct>           <dbl>     <int> <fct>     
## 1 25         25758      SOPO       1653             111.      2017 CUNDINAMA~
## # ... with 3 more variables: Shape_Leng <dbl>, Shape_Area <dbl>,
## #   geometry <MULTIPOLYGON [°]>

Ahora verifiquemos la última declaración:

datos %>% filter (MUNICIPIO =="SOPO") ->  med_datos
med_datos
## # A tibble: 160 x 17
##    COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO SUBGRUPO CULTIVO
##      <dbl> <chr>          <dbl> <chr>     <chr> <chr>    <chr>  
##  1      25 CUNDINAMARCA   25758 SOPO      HORT~ AJO      AJO    
##  2      25 CUNDINAMARCA   25758 SOPO      HORT~ AJO      AJO    
##  3      25 CUNDINAMARCA   25758 SOPO      HORT~ AJO      AJO    
##  4      25 CUNDINAMARCA   25758 SOPO      HORT~ AJO      AJO    
##  5      25 CUNDINAMARCA   25758 SOPO      HORT~ AJO      AJO    
##  6      25 CUNDINAMARCA   25758 SOPO      HORT~ AJO      AJO    
##  7      25 CUNDINAMARCA   25758 SOPO      HORT~ AJO      AJO    
##  8      25 CUNDINAMARCA   25758 SOPO      HORT~ AJO      AJO    
##  9      25 CUNDINAMARCA   25758 SOPO      HORT~ AJO      AJO    
## 10      25 CUNDINAMARCA   25758 SOPO      HORT~ AJO      AJO    
## # ... with 150 more rows, and 10 more variables:
## #   DESAGREGACION_REGIONAL_SISTEMA_PRODUCTIVO <chr>, YEAR <dbl>, PERIODO <chr>,
## #   Area_Siembra <dbl>, Area_Cosecha <dbl>, Produccion <dbl>,
## #   Rendimiento <dbl>, ESTADO <chr>, N_CIENTIFICO <chr>, CICLO_CULTIVO <chr>
class(med_datos$COD_MUN)
## [1] "numeric"

Para poder hacer la unión, necesitamos cambiar tanto el tipo de datos como el contenido del código que identifica a cada municipio. Para esta tarea, es una buena idea crear una copia de los datos estadísticos originales. Con este enfoque, cualquier movimiento falso no estropeará los datos originales.

Avancemos paso a paso:

datos2 <- datos
datos2$TEMP <- as.character(datos2$COD_MUN)
datos2$MPIO_CCDGO <- as.factor(paste(datos2$TEMP, sep=""))
datos2$TEMP <- as.factor(datos2$TEMP)
class(datos2$TEMP)
## [1] "factor"
class(datos2$MPIO_CCDGO)
## [1] "factor"
head(datos2)
## # A tibble: 6 x 19
##   COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO SUBGRUPO CULTIVO DESAGREGACION_R~
##     <dbl> <chr>          <dbl> <chr>     <chr> <chr>    <chr>   <chr>           
## 1      25 CUNDINAMARCA   25754 SOACHA    HORT~ ACELGA   ACELGA  ACELGA          
## 2      25 CUNDINAMARCA   25214 COTA      HORT~ ACELGA   ACELGA  ACELGA          
## 3      25 CUNDINAMARCA   25754 SOACHA    HORT~ ACELGA   ACELGA  ACELGA          
## 4      25 CUNDINAMARCA   25214 COTA      HORT~ ACELGA   ACELGA  ACELGA          
## 5      25 CUNDINAMARCA   25754 SOACHA    HORT~ ACELGA   ACELGA  ACELGA          
## 6      25 CUNDINAMARCA   25126 CAJICA    HORT~ ACELGA   ACELGA  ACELGA          
## # ... with 11 more variables: YEAR <dbl>, PERIODO <chr>, Area_Siembra <dbl>,
## #   Area_Cosecha <dbl>, Produccion <dbl>, Rendimiento <dbl>, ESTADO <chr>,
## #   N_CIENTIFICO <chr>, CICLO_CULTIVO <chr>, TEMP <fct>, MPIO_CCDGO <fct>

Asegúrese de verificar, en el objeto de datos2, las características del nuevo atributo MPIO_CCDGO.

Ahora, filtremos un solo año y seleccionemos solo dos atributos relevantes:

datos2 %>% filter(CULTIVO == "CAFE")  -> datos3
head(datos3)
## # A tibble: 6 x 19
##   COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO SUBGRUPO CULTIVO DESAGREGACION_R~
##     <dbl> <chr>          <dbl> <chr>     <chr> <chr>    <chr>   <chr>           
## 1      25 CUNDINAMARCA   25878 VIOTA     OTRO~ CAFE     CAFE    CAFE            
## 2      25 CUNDINAMARCA   25662 SAN JUAN~ OTRO~ CAFE     CAFE    CAFE            
## 3      25 CUNDINAMARCA   25245 EL COLEG~ OTRO~ CAFE     CAFE    CAFE            
## 4      25 CUNDINAMARCA   25596 QUIPILE   OTRO~ CAFE     CAFE    CAFE            
## 5      25 CUNDINAMARCA   25885 YACOPI    OTRO~ CAFE     CAFE    CAFE            
## 6      25 CUNDINAMARCA   25394 LA PALMA  OTRO~ CAFE     CAFE    CAFE            
## # ... with 11 more variables: YEAR <dbl>, PERIODO <chr>, Area_Siembra <dbl>,
## #   Area_Cosecha <dbl>, Produccion <dbl>, Rendimiento <dbl>, ESTADO <chr>,
## #   N_CIENTIFICO <chr>, CICLO_CULTIVO <chr>, TEMP <fct>, MPIO_CCDGO <fct>
class(datos3)
## [1] "spec_tbl_df" "tbl_df"      "tbl"         "data.frame"
datos4 <- datos3 %>% dplyr::select(MUNICIPIO, MPIO_CCDGO, YEAR, Produccion, Rendimiento) 
datos4
## # A tibble: 821 x 5
##    MUNICIPIO           MPIO_CCDGO  YEAR Produccion Rendimiento
##    <chr>               <fct>      <dbl>      <dbl>       <dbl>
##  1 VIOTA               25878       2007       3730          83
##  2 SAN JUAN DE RIOSECO 25662       2007       2018          74
##  3 EL COLEGIO          25245       2007       1199          64
##  4 QUIPILE             25596       2007        823          44
##  5 YACOPI              25885       2007       1831         103
##  6 LA PALMA            25394       2007       1752         113
##  7 VERGARA             25862       2007       1254         112
##  8 CHAGUANI            25168       2007       1118          92
##  9 LA MESA             25386       2007        782          64
## 10 CAPARRAPI           25148       2007       1185          98
## # ... with 811 more rows
datos4 %>% 
  gather("YEAR", "Produccion", "Rendimiento" , key = variable, value = number)
## # A tibble: 2,463 x 4
##    MUNICIPIO           MPIO_CCDGO variable number
##    <chr>               <fct>      <chr>     <dbl>
##  1 VIOTA               25878      YEAR       2007
##  2 SAN JUAN DE RIOSECO 25662      YEAR       2007
##  3 EL COLEGIO          25245      YEAR       2007
##  4 QUIPILE             25596      YEAR       2007
##  5 YACOPI              25885      YEAR       2007
##  6 LA PALMA            25394      YEAR       2007
##  7 VERGARA             25862      YEAR       2007
##  8 CHAGUANI            25168      YEAR       2007
##  9 LA MESA             25386      YEAR       2007
## 10 CAPARRAPI           25148      YEAR       2007
## # ... with 2,453 more rows
head(datos4)
## # A tibble: 6 x 5
##   MUNICIPIO           MPIO_CCDGO  YEAR Produccion Rendimiento
##   <chr>               <fct>      <dbl>      <dbl>       <dbl>
## 1 VIOTA               25878       2007       3730          83
## 2 SAN JUAN DE RIOSECO 25662       2007       2018          74
## 3 EL COLEGIO          25245       2007       1199          64
## 4 QUIPILE             25596       2007        823          44
## 5 YACOPI              25885       2007       1831         103
## 6 LA PALMA            25394       2007       1752         113

Esta parte es muy importante, debido a que implica varios pasos para poder convertir la tabla de atributos de formato largo a formato ancho, lo cual podemos realizar con el siguiente codigo:

datos4 %>% 
  group_by(MPIO_CCDGO) %>% 
  mutate(Visit = 1:n()) %>% 
  gather("YEAR", "Produccion", "Rendimiento", key = variable, value = number) %>% 
  unite(combi, variable, Visit) %>%
  spread(combi, number) -> datos5
datos5
## # A tibble: 71 x 38
## # Groups:   MPIO_CCDGO [116]
##    MUNICIPIO MPIO_CCDGO Produccion_1 Produccion_10 Produccion_11 Produccion_12
##    <chr>     <fct>             <dbl>         <dbl>         <dbl>         <dbl>
##  1 AGUA DE ~ 25001                 3            NA            NA            NA
##  2 ALBAN     25019               513           341           395           406
##  3 ANAPOIMA  25035                92            75            59            57
##  4 ANOLAIMA  25040               622           755          1280          1251
##  5 APULO     25599                50            44            43            41
##  6 ARBELAEZ  25053               230           345           334           311
##  7 BELTRAN   25086                68            53            42            38
##  8 BITUIMA   25095               264           374           344           363
##  9 CABRERA   25120                 4            21            18            19
## 10 CACHIPAY  25123               762           671           940           859
## # ... with 61 more rows, and 32 more variables: Produccion_2 <dbl>,
## #   Produccion_3 <dbl>, Produccion_4 <dbl>, Produccion_5 <dbl>,
## #   Produccion_6 <dbl>, Produccion_7 <dbl>, Produccion_8 <dbl>,
## #   Produccion_9 <dbl>, Rendimiento_1 <dbl>, Rendimiento_10 <dbl>,
## #   Rendimiento_11 <dbl>, Rendimiento_12 <dbl>, Rendimiento_2 <dbl>,
## #   Rendimiento_3 <dbl>, Rendimiento_4 <dbl>, Rendimiento_5 <dbl>,
## #   Rendimiento_6 <dbl>, Rendimiento_7 <dbl>, Rendimiento_8 <dbl>,
## #   Rendimiento_9 <dbl>, YEAR_1 <dbl>, YEAR_10 <dbl>, YEAR_11 <dbl>,
## #   YEAR_12 <dbl>, YEAR_2 <dbl>, YEAR_3 <dbl>, YEAR_4 <dbl>, YEAR_5 <dbl>,
## #   YEAR_6 <dbl>, YEAR_7 <dbl>, YEAR_8 <dbl>, YEAR_9 <dbl>

Podemos ver aca todos los datos en una ventana nueva para verificarlos:

view(datos5)

También haremos una copia de la colección de características simples (de nuevo, solo en caso de un movimiento falso):

ant_munic2 <- ant_munic

Ahora es hora de unirse:

ant_munic_stat = left_join(ant_munic2, datos5, by="MPIO_CCDGO")
summary(ant_munic_stat)
##  DPTO_CCDGO   MPIO_CCDGO         MPIO_CNMBR    MPIO_CRSLC   MPIO_NAREA     
##  25:116     25001  :  1   AGUA DE DIOS:  1   1700   : 8   Min.   :  43.06  
##             25019  :  1   ALBÁN       :  1   1799   : 7   1st Qu.:  91.80  
##             25035  :  1   ANAPOIMA    :  1   1899   : 6   Median : 131.63  
##             25040  :  1   ANOLAIMA    :  1   1593   : 4   Mean   : 192.85  
##             25053  :  1   APULO       :  1   1600   : 4   3rd Qu.: 214.78  
##             25086  :  1   ARBELÁEZ    :  1   1900   : 4   Max.   :1196.75  
##             (Other):110   (Other)     :110   (Other):83                    
##    MPIO_NANO           DPTO_CNMBR    Shape_Leng       Shape_Area      
##  Min.   :2017   CUNDINAMARCA:116   Min.   :0.2928   Min.   :0.003511  
##  1st Qu.:2017                      1st Qu.:0.4988   1st Qu.:0.007486  
##  Median :2017                      Median :0.6205   Median :0.010735  
##  Mean   :2017                      Mean   :0.7064   Mean   :0.015721  
##  3rd Qu.:2017                      3rd Qu.:0.8325   3rd Qu.:0.017503  
##  Max.   :2017                      Max.   :1.9999   Max.   :0.097502  
##                                                                       
##   MUNICIPIO          Produccion_1    Produccion_10    Produccion_11   
##  Length:116         Min.   :   0.0   Min.   :  16.0   Min.   :  15.0  
##  Class :character   1st Qu.:  11.0   1st Qu.: 113.5   1st Qu.: 118.0  
##  Mode  :character   Median : 205.0   Median : 331.0   Median : 337.0  
##                     Mean   : 475.1   Mean   : 468.5   Mean   : 518.8  
##                     3rd Qu.: 753.5   3rd Qu.: 695.0   3rd Qu.: 678.0  
##                     Max.   :3730.0   Max.   :3430.0   Max.   :4114.0  
##                     NA's   :45       NA's   :49       NA's   :51      
##  Produccion_12     Produccion_2       Produccion_3     Produccion_4    
##  Min.   :  19.0   Min.   :    0.00   Min.   :   0.0   Min.   :   0.00  
##  1st Qu.: 125.0   1st Qu.:   47.75   1st Qu.:  37.0   1st Qu.:  52.75  
##  Median : 363.0   Median :  368.00   Median : 252.5   Median : 247.00  
##  Mean   : 511.2   Mean   : 1118.26   Mean   : 530.7   Mean   : 532.10  
##  3rd Qu.: 695.5   3rd Qu.: 1092.00   3rd Qu.: 845.5   3rd Qu.: 828.25  
##  Max.   :3979.0   Max.   :11232.00   Max.   :4359.0   Max.   :4209.00  
##  NA's   :53       NA's   :46         NA's   :46       NA's   :46       
##   Produccion_5     Produccion_6   Produccion_7     Produccion_8 
##  Min.   :   2.0   Min.   :   0   Min.   :   7.0   Min.   :   7  
##  1st Qu.:  49.0   1st Qu.:  48   1st Qu.:  62.0   1st Qu.:  74  
##  Median : 130.5   Median : 185   Median : 189.0   Median : 217  
##  Mean   : 469.2   Mean   : 449   Mean   : 363.9   Mean   : 366  
##  3rd Qu.: 707.0   3rd Qu.: 600   3rd Qu.: 501.0   3rd Qu.: 505  
##  Max.   :5000.0   Max.   :4161   Max.   :2333.0   Max.   :2476  
##  NA's   :46       NA's   :47     NA's   :47       NA's   :47    
##   Produccion_9    Rendimiento_1    Rendimiento_10   Rendimiento_11  
##  Min.   :  17.0   Min.   :  1.00   Min.   :  1.00   Min.   :  1.00  
##  1st Qu.: 103.0   1st Qu.: 24.00   1st Qu.: 84.50   1st Qu.: 84.00  
##  Median : 305.0   Median : 54.00   Median : 95.00   Median :101.00  
##  Mean   : 458.7   Mean   : 56.41   Mean   : 87.67   Mean   : 97.91  
##  3rd Qu.: 687.2   3rd Qu.: 81.00   3rd Qu.:103.00   3rd Qu.:126.00  
##  Max.   :3355.0   Max.   :154.00   Max.   :121.00   Max.   :182.00  
##  NA's   :48       NA's   :47       NA's   :49       NA's   :51      
##  Rendimiento_12  Rendimiento_2     Rendimiento_3    Rendimiento_4   
##  Min.   : 62.0   Min.   :   1.00   Min.   :  1.00   Min.   :  5.00  
##  1st Qu.: 86.0   1st Qu.:  35.25   1st Qu.: 42.25   1st Qu.: 52.00  
##  Median :103.0   Median :  70.00   Median : 60.00   Median : 72.00  
##  Mean   :112.1   Mean   : 164.66   Mean   : 63.16   Mean   : 71.65  
##  3rd Qu.:129.0   3rd Qu.: 111.75   3rd Qu.: 88.75   3rd Qu.: 98.00  
##  Max.   :187.0   Max.   :1075.00   Max.   :133.00   Max.   :133.00  
##  NA's   :53      NA's   :46        NA's   :46       NA's   :47      
##  Rendimiento_5    Rendimiento_6    Rendimiento_7    Rendimiento_8   
##  Min.   :  1.00   Min.   :  1.00   Min.   :  4.00   Min.   :  9.00  
##  1st Qu.:  8.00   1st Qu.:  5.75   1st Qu.: 49.00   1st Qu.: 55.00  
##  Median : 47.00   Median : 36.50   Median : 64.00   Median : 69.00  
##  Mean   : 45.37   Mean   : 41.51   Mean   : 62.99   Mean   : 73.96  
##  3rd Qu.: 74.00   3rd Qu.: 75.00   3rd Qu.: 77.00   3rd Qu.: 87.00  
##  Max.   :136.00   Max.   :114.00   Max.   :128.00   Max.   :139.00  
##  NA's   :46       NA's   :48       NA's   :47       NA's   :47      
##  Rendimiento_9        YEAR_1        YEAR_10        YEAR_11        YEAR_12    
##  Min.   :  1.00   Min.   :2007   Min.   :2016   Min.   :2017   Min.   :2018  
##  1st Qu.: 76.75   1st Qu.:2007   1st Qu.:2016   1st Qu.:2017   1st Qu.:2018  
##  Median : 89.00   Median :2007   Median :2016   Median :2017   Median :2018  
##  Mean   : 78.40   Mean   :2007   Mean   :2016   Mean   :2017   Mean   :2018  
##  3rd Qu.: 99.00   3rd Qu.:2007   3rd Qu.:2016   3rd Qu.:2017   3rd Qu.:2018  
##  Max.   :118.00   Max.   :2012   Max.   :2018   Max.   :2018   Max.   :2018  
##  NA's   :48       NA's   :45     NA's   :49     NA's   :51     NA's   :53    
##      YEAR_2         YEAR_3         YEAR_4         YEAR_5         YEAR_6    
##  Min.   :2008   Min.   :2009   Min.   :2010   Min.   :2011   Min.   :2012  
##  1st Qu.:2008   1st Qu.:2009   1st Qu.:2010   1st Qu.:2011   1st Qu.:2012  
##  Median :2008   Median :2009   Median :2010   Median :2011   Median :2012  
##  Mean   :2008   Mean   :2009   Mean   :2010   Mean   :2011   Mean   :2012  
##  3rd Qu.:2008   3rd Qu.:2009   3rd Qu.:2010   3rd Qu.:2011   3rd Qu.:2012  
##  Max.   :2012   Max.   :2013   Max.   :2014   Max.   :2015   Max.   :2016  
##  NA's   :46     NA's   :46     NA's   :46     NA's   :46     NA's   :47    
##      YEAR_7         YEAR_8         YEAR_9              geometry  
##  Min.   :2013   Min.   :2014   Min.   :2015   MULTIPOLYGON :116  
##  1st Qu.:2013   1st Qu.:2014   1st Qu.:2015   epsg:4326    :  0  
##  Median :2013   Median :2014   Median :2015   +proj=long...:  0  
##  Mean   :2013   Mean   :2014   Mean   :2015                      
##  3rd Qu.:2013   3rd Qu.:2014   3rd Qu.:2015                      
##  Max.   :2017   Max.   :2018   Max.   :2018                      
##  NA's   :47     NA's   :47     NA's   :48

5. Graficando

#install.packages("RColorBrewer")
library(RColorBrewer)
library(leaflet)

Tracemos los municipios con su producción de café correspondiente para un solo año:

bins <- c(0, 250, 500, 1000, 2000, 5000, 10000, 15000)
pal <- colorBin("YlOrRd", domain = ant_munic_stat$Produccion_12, bins = bins)

  mapa <- leaflet(data = ant_munic_stat) %>%
  addTiles() %>%
  addPolygons(label = ~Produccion_12,
              popup = ~MPIO_CNMBR,
              fillColor = ~pal(Produccion_12),
              color = "#444444",
              weight = 1,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.5,
              highlightOptions = highlightOptions(color = "white", weight = 2, bringToFront = TRUE)
              ) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend("bottomright", pal = pal, values = ~Produccion_12,
    title = "Produccion de Cafe en Cundinamarca [Ton] (2018)",
    opacity = 1
  )
mapa
sessionInfo()
## R version 3.6.0 (2019-04-26)
## 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   
## [3] 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] leaflet_2.0.3       RColorBrewer_1.1-2  cowplot_1.0.0      
##  [4] ggrepel_0.8.2       GSODR_2.1.0         rnaturalearth_0.1.0
##  [7] viridis_0.5.1       viridisLite_0.3.0   sf_0.9-0           
## [10] raster_3.1-5        maptools_0.9-9      rgeos_0.5-2        
## [13] sp_1.4-1            forcats_0.5.0       stringr_1.4.0      
## [16] dplyr_0.8.5         purrr_0.3.3         readr_1.3.1        
## [19] tidyr_1.0.2         tibble_2.1.3        ggplot2_3.3.0      
## [22] tidyverse_1.3.0     here_0.1           
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-139            fs_1.4.0                lubridate_1.7.4        
##  [4] httr_1.4.1              rprojroot_1.3-2         tools_3.6.0            
##  [7] backports_1.1.5         utf8_1.1.4              R6_2.4.1               
## [10] KernSmooth_2.23-15      DBI_1.1.0               colorspace_1.4-1       
## [13] withr_2.1.2             tidyselect_1.0.0        gridExtra_2.3          
## [16] compiler_3.6.0          cli_2.0.2               rvest_0.3.5            
## [19] xml2_1.3.0              labeling_0.3            scales_1.1.0           
## [22] classInt_0.4-2          digest_0.6.25           foreign_0.8-71         
## [25] rmarkdown_2.1           pkgconfig_2.0.3         htmltools_0.4.0        
## [28] fastmap_1.0.1           dbplyr_1.4.2            htmlwidgets_1.5.1      
## [31] rlang_0.4.5             readxl_1.3.1            rstudioapi_0.11        
## [34] shiny_1.4.0             generics_0.0.2          farver_2.0.3           
## [37] jsonlite_1.6.1          crosstalk_1.0.0         magrittr_1.5           
## [40] Rcpp_1.0.3              munsell_0.5.0           fansi_0.4.1            
## [43] lifecycle_0.2.0         stringi_1.4.6           yaml_2.2.1             
## [46] grid_3.6.0              promises_1.1.0          crayon_1.3.4           
## [49] lattice_0.20-38         haven_2.2.0             hms_0.5.3              
## [52] knitr_1.28              pillar_1.4.3            codetools_0.2-16       
## [55] reprex_0.3.0            glue_1.3.1              evaluate_0.14          
## [58] leaflet.providers_1.9.0 data.table_1.12.8       modelr_0.1.6           
## [61] vctrs_0.2.3             httpuv_1.5.2            cellranger_1.1.0       
## [64] gtable_0.3.0            assertthat_0.2.1        xfun_0.12              
## [67] mime_0.9                xtable_1.8-4            broom_0.5.5            
## [70] e1071_1.7-3             later_1.0.0             class_7.3-15           
## [73] units_0.6-6             ellipsis_0.3.0