Estadisticas municipales de la agricultura en Cauca

1. Why this notebook?

Este es un Cuaderno de R Markdown que ilustra estadisticas agricolas para el departamento de Cauca en Colombia, tiene como objetivo la comprención de los conceptos básicos de geomática útiles para la agronomía como respuesta a lo enseñado en el curso de Geomatica Básica de la Universidad Nacional de Colombia.

2. Funciones de GIS

La exploración de estadísticas no espaciales es fundamental para comprender lo que sucede en los territorios. Varias bibliotecas de 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 los municipios cuyos rendimientos de cultivos son sobresalientes (o, alternativamente, más bajos que el promedio). Para realizar dicha exploración, necesitamos unir datos no espaciales con datos espaciales.

Además, también podríamos explorar combinaciones espaciales. Estas operaciones se basan en la intersección entre dos objetos espaciales, a menudo puntos y 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!

Comencemos eliminando el contenido de la memoria:

rm(list=ls())

Ahora, instalemos las bibliotecas que necesitamos. Tenga en cuenta que, en el siguiente fragmento, cualquier paquete se instala solo si no se ha instalado previamente. Se recomienda correr uno por uno en la consola para detectar posibles errores

#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)

Ahora, carguemos las bibliotecas.

library(here)

here() starts at C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/trabajo7

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.2     v purrr   0.3.4
v tibble  3.0.3     v dplyr   1.0.2
v tidyr   1.1.2     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x tidyr::extract() masks raster::extract()
x dplyr::filter()  masks stats::filter()
x dplyr::lag()     masks stats::lag()
x dplyr::select()  masks raster::select()

library(rgeos)

rgeos version: 0.5-5, (SVN revision 640)
 GEOS runtime version: 3.8.0-CAPI-1.13.1 
 Linking to sp version: 1.4-2 
 Polygon checking: TRUE 

library(maptools)

Checking rgeos availability: TRUE

library(raster)
library(sf)

Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

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)

3. Explorando las estadísticas agrícolas en Cauca.

Previamente, he descargado datos estadísticos, en formato csv, sobre Estadisticas Municipales Agropecuarias a mi computadora. Luego, he usado Excel para eliminar filas para municipios en departamentos diferentes a Cauca, cree una nueva fila con el nombre COD donde esta registrado el ID+1 no que nos será útil más adenate. Guardé el archivo con las filas restantes como un archivo local con el nombre EVA_cauca.csv en mi computadora.

Leamos el archivo csv con “estadisticas municipales agropecuarias” para Cauca.

setwd("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/trabajo7")
datos<-read_csv("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/trabajo7/EVA_cauca.csv")

Parsed with column specification:
cols(
  ID = col_double(),
  COD = col_character(),
  MUN = col_character(),
  GRUPO = col_character(),
  SUBGRUPO = col_character(),
  CULTIVO = col_character(),
  DESAGR = col_character(),
  YEAR = col_double(),
  PERIODO = col_character(),
  HA_SIEMBRA = col_double(),
  HA_COSECHA = col_double(),
  TON_PROD = col_double(),
  REND = col_double(),
  ESTADO = col_character(),
  NOMBRE = col_character(),
  CICLO = col_character()
)

Veamos cuáles son los atributos de los datos.

head(datos)

NA

tail(datos)

Teniendo en cuenta que cada municipio dispone de estadísticas sobre superficie sembrada, superficie cosechada y rendimiento para diferentes cultivos en diferentes años. El atributo SUBGRUPO y CULTIVO parecen referirse a lo mismo (es decir, un cultivo). Los cultivos se clasifican además en un GRUPO determinado.

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

Usaremos 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(MUN, GRUPO) %>%
  summarise(rend_prom = mean(REND, na.rm = TRUE)) -> rend_resumen

`summarise()` regrouping output by 'MUN' (override with `.groups` argument)

### Revisaremos solo los 6 primeros datos.
head(rend_resumen)

También podemos calcular el rendimiento promedio por GRUPO en los municipios de Cauca:

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

`summarise()` ungrouping output (override with `.groups` argument)

rend_cauca

Nótese que los mayores rendimientos corresponden a HORTALIZAS, FLORES Y FOLLAJES y OTROS PERMANENTES.

Luego, busquemos cuáles son los municipios con mayor rendimiento para cada grupo de cultivos en 2018:

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

ningun argumento finito para max; retornando -Inf`summarise()` regrouping output by 'GRUPO' (override with `.groups` argument)

rend_max_18

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

datos %>% 
  filter(YEAR==2018) %>% 
  group_by(GRUPO, MUN) %>%
  summarize(max_HA_COSECHA = max(HA_COSECHA, na.rm = TRUE)) %>%
    slice(which.max(max_HA_COSECHA)) -> HA_COSECHA_max

`summarise()` regrouping output by 'GRUPO' (override with `.groups` argument)

HA_COSECHA_max

Nótese que la mayor superficie cosechada en 2018 correspondió a OTROS PERMANENTES en EL TAMBO, Quizás sepa que la economía de EL TAMBO está soportada principalmente por la producción cafetera intercalada con plátano.

Primero, seleccionemos la producción de cafe (toneladas) en EL TAMBO para cada año:

datos %>% 
  filter(MUN=="EL TAMBO" & SUBGRUPO=="CAFE") %>% 
  group_by(YEAR, CULTIVO) ->  TAMBO_CAFE

TAMBO_CAFE

Hagamos una exploración gráfica básica usando la biblioteca ggplot2:

g<- ggplot(aes(x=YEAR, y=TON_PROD/1000), data = TAMBO_CAFE) + geom_bar(stat='identity') + labs(y='Producción de Cafe [Ton x 1000]', x='Año')

g + ggtitle("Evolución de la producción de café en El Tambo de 2007 a 2018") + labs(caption= "Basado en datos de EMA (DANE, 2018)")

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

datos %>% 
  filter(YEAR==2018) %>% 
  group_by(GRUPO) %>%
  summarize(sum_HA_COSECHA = sum(HA_COSECHA, na.rm = TRUE)) %>%
     arrange(desc(sum_HA_COSECHA)) -> total_HA_COSECHA

`summarise()` ungrouping output (override with `.groups` argument)

total_HA_COSECHA

Podemos ver que otros cultivos permanentes también tuvieron la mayor proporción de área cosechada en 2018 para Cauca. ¿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.

También podemos buscar esta información en los propios datos:

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

`summarise()` ungrouping output (override with `.groups` argument)

total_cosecha

¡Tenga en cuenta que el departamento del Cauca hace parte del llamado nuevo eje cafetero!

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

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

`summarise()` regrouping output by 'CULTIVO' (override with `.groups` argument)

area_cosecha2

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 trama se verá desordenada.

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

Ahora es el momento de trazar usando la biblioteca ggplot2:

g <- ggplot(aes(x=CROP, y=sum_HA_COSECHA), data = total_HA_COSECHA) + geom_bar(stat='identity') + labs(y='Área total cosechada [Ha]', x= 'Cultivo')

g+ ggtitle("Superficie total cosechada por grupos de cultivos en 2018 para Cauca") + theme(plot.title = element_text(hjust = 0.5)) +
   labs(caption= "Basado en datos de EMA (DANE, 2018)") 

Ya que uno de los principales productos agricolas del Cauca es el cafe, odservaremos el rendimiento máximo de este producto en los diferentes municipios para el año 2018, como ha sido la evolución de su producción a través de los años y por ultimo revisaremos cuales han sido los municipios con mayor área cosecheada en el 2018. Primero el rendimiento máximo del cafe en los diferentes municipios para el año 2018 y cuales municipios tuvieron mayor Área de cosecha de Cafe en 2018:

datos %>% 
  filter( YEAR==2018 &GRUPO=="OTROS PERMANENTES" & SUBGRUPO=="CAFE") %>% 
   group_by(MUN,GRUPO,CULTIVO,YEAR) %>%
  summarize(max_rend = max(REND, na.rm = TRUE)) %>%
    slice(which.max(max_rend)) -> rendcafe_max_18

`summarise()` regrouping output by 'MUN', 'GRUPO', 'CULTIVO' (override with `.groups` argument)

rendcafe_max_18

Ahora miraremos la evolución de la producción del año 2007 a al 2018:

datos %>% 
  filter(SUBGRUPO=="CAFE") %>% 
  group_by(YEAR, CULTIVO) ->  CAUCA_CAFE

head(CAUCA_CAFE)

g<- ggplot(aes(x=YEAR, y=TON_PROD/1000), data = CAUCA_CAFE) + geom_bar(stat='identity') + labs(y='Producción de Cafe [Ton x 1000]', x='Año')

g + ggtitle("Evolución de la producción de café en Cauca de 2017 a 2018") + labs(caption= "Basado en datos de EMA (DANE, 2018)")

Hola veremos lo municipios de Cauca con mayor área cosechada para el 2018

datos %>% 
  filter(YEAR==2018 & SUBGRUPO=="CAFE") %>% 
  group_by(MUN, CULTIVO) %>%
  summarize(sum_cosecha = sum(HA_COSECHA, na.rm = TRUE)) %>%
  filter(sum_cosecha>4000) %>%
  arrange(desc(sum_cosecha)) -> MUN2

`summarise()` regrouping output by 'MUN' (override with `.groups` argument)

MUN2

NA

MUN2$MUN <- abbreviate(MUN2$MUN , 4)
MUN2

g <- ggplot(aes(x=MUN, y=sum_cosecha), data = MUN2) + geom_bar(stat='identity') + labs(y='Área total cosechada [Ha]', x= 'Municipio')

g+ ggtitle("Municipos con mayor área cosechada de cafe en Cauca (2018)") + theme(plot.title = element_text(hjust = 0.3)) +
   labs(caption= "Basado en datos de EMA (DANE, 2018)")

4. Incorporación de las estadísticas agrícolas a los municipios

Comencemos a leer los datos administrativos de Cauca, dados por Marco Geoestadistico Departamental que se encuentra disponible en DANE Geoportal, usando la biblioteca sf:

cau_munic <- sf::st_read("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/Departamento/19_CAUCA/Cauca/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\rojas\OneDrive\Escritorio\R estudio\Geomatica\Departamento\19_CAUCA\Cauca\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
Simple feature collection with 42 features and 9 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -77.92834 ymin: 0.9580285 xmax: -75.74782 ymax: 3.328941
geographic CRS: WGS 84

¿Qué es ant_munic?

cau_munic

Simple feature collection with 42 features and 9 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -77.92834 ymin: 0.9580285 xmax: -75.74782 ymax: 3.328941
geographic CRS: WGS 84
First 10 features:
   DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR                         MPIO_CRSLC
1          19      19001      POPAYÁN                               1537
2          19      19022     ALMAGUER                               1799
3          19      19050      ARGELIA Ordenanza 2 de Noviembre 8 de 1967
4          19      19075       BALBOA  Ordenanza 1 de Octubre 20 de 1967
5          19      19100      BOLÍVAR                               1793
6          19      19110 BUENOS AIRES                               1851
7          19      19130      CAJIBÍO                               1824
8          19      19137      CALDONO                               1746
9          19      19142       CALOTO                               1543
10         19      19212      CORINTO                               1868
   MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area
1    480.1798      2017      CAUCA  1.3729371 0.03897000
2    236.7148      2017      CAUCA  0.7102530 0.01919657
3    775.7878      2017      CAUCA  1.2241529 0.06288063
4    413.2112      2017      CAUCA  1.0163274 0.03349049
5    797.8203      2017      CAUCA  1.8027866 0.06468337
6    435.0623      2017      CAUCA  1.6333631 0.03532030
7    552.2877      2017      CAUCA  1.3822054 0.04482019
8    354.7387      2017      CAUCA  1.1029614 0.02880277
9    266.5987      2017      CAUCA  1.3621191 0.02164894
10   325.9156      2017      CAUCA  0.8276564 0.02647947
                         geometry
1  POLYGON ((-76.74145 2.57230...
2  POLYGON ((-76.80864 1.98674...
3  POLYGON ((-77.3107 2.508326...
4  POLYGON ((-77.1431 2.156948...
5  POLYGON ((-76.92984 2.11058...
6  POLYGON ((-76.76916 3.23128...
7  POLYGON ((-76.83997 2.76373...
8  POLYGON ((-76.36326 2.90734...
9  POLYGON ((-76.42866 3.20990...
10 POLYGON ((-76.28278 3.22879...

Tenga en cuenta que ant_munic es una colección de funciones simple. Asegúrese de revisar este enlace para comprender qué es una característica simple.

Tenga en cuenta también que los datos utilizan el sistema de referencia de coordenadas geográficas WGS1984 (es decir, el código 4326 epsg).

Podemos usar la función left_join para unir los municipios y las estadísticas agrícolas seleccionadas.

Necesitamos un atributo común (o variable compartida) en el que basar la unión. El mejor atributo es una identificación. En cau_munic, el atributo MPIO_CCDGO parece estar bien (lee 19001 para Popayán). En datos, el atributo correspondiente es COD_MUN (lee 19001 para Popayán).

Verifiquemos la última declaración:

datos %>% filter (MUN =="POPAYAN") ->  med_datos

med_datos

class(med_datos$COD)

[1] "character"

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.

Procedamos paso a paso:

datos2 <- datos

datos2$TEMP <-  as.character(datos2$COD)

datos2$MPIO_CCDGO <- as.factor(paste(datos2$TEMP, sep=""))

head(datos2)

Asegúrese de verificar, en el objeto datos 2, 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)

class(datos3)

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

datos4 <- datos3 %>% dplyr::select(MUN, MPIO_CCDGO, YEAR, TON_PROD, REND)

datos4

datos4 %>% 
  gather("YEAR", "TON_PROD", "REND" , key = variable, value = number)

head(datos4)

Ésta es una tarea clave. Implica varios pasos para poder convertir la tabla de atributos de formato largo a formato ancho. Más información sobre estos pasos aquí.

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

head(datos5)

tail(datos5)

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

cau_munic2 <- cau_munic

Ahora es el momento de unirse:

cau_munic_stat = left_join(cau_munic2, datos5, by="MPIO_CCDGO")

summary(cau_munic_stat)

  DPTO_CCDGO         MPIO_CCDGO         MPIO_CNMBR       
 Length:42          Length:42          Length:42         
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                         
                                                         
                                                         
                                                         
  MPIO_CRSLC          MPIO_NAREA        MPIO_NANO     DPTO_CNMBR       
 Length:42          Min.   :  56.65   Min.   :2017   Length:42         
 Class :character   1st Qu.: 214.88   1st Qu.:2017   Class :character  
 Mode  :character   Median : 457.62   Median :2017   Mode  :character  
                    Mean   : 743.88   Mean   :2017                     
                    3rd Qu.: 770.66   3rd Qu.:2017                     
                    Max.   :3619.37   Max.   :2017                     
                                                                       
   Shape_Leng       Shape_Area           MUN                REND_1      
 Min.   :0.3676   Min.   :0.004592   Length:42          Min.   :0.5200  
 1st Qu.:0.9546   1st Qu.:0.017436   Class :character   1st Qu.:0.6600  
 Median :1.2427   Median :0.037145   Mode  :character   Median :0.8300  
 Mean   :1.4087   Mean   :0.060343                      Mean   :0.8574  
 3rd Qu.:1.6100   3rd Qu.:0.062469                      3rd Qu.:0.9900  
 Max.   :3.6742   Max.   :0.293595                      Max.   :1.6400  
                                                        NA's   :7       
    REND_10         REND_11         REND_12          REND_2      
 Min.   :0.900   Min.   :0.660   Min.   :0.620   Min.   :0.4700  
 1st Qu.:1.040   1st Qu.:0.840   1st Qu.:0.860   1st Qu.:0.6150  
 Median :1.110   Median :1.100   Median :1.100   Median :0.7500  
 Mean   :1.122   Mean   :1.097   Mean   :1.106   Mean   :0.7926  
 3rd Qu.:1.190   3rd Qu.:1.260   3rd Qu.:1.290   3rd Qu.:0.9050  
 Max.   :1.520   Max.   :1.820   Max.   :1.870   Max.   :1.4600  
 NA's   :11      NA's   :11      NA's   :12      NA's   :8       
     REND_3           REND_4           REND_5           REND_6      
 Min.   :0.4400   Min.   :0.4500   Min.   :0.4500   Min.   :0.4500  
 1st Qu.:0.5950   1st Qu.:0.5875   1st Qu.:0.5750   1st Qu.:0.7450  
 Median :0.7200   Median :0.7200   Median :0.6800   Median :1.0500  
 Mean   :0.7809   Mean   :0.7872   Mean   :0.7023   Mean   :0.9061  
 3rd Qu.:0.8550   3rd Qu.:0.9325   3rd Qu.:0.8300   3rd Qu.:1.1000  
 Max.   :1.8200   Max.   :1.8700   Max.   :1.0600   Max.   :1.1000  
 NA's   :10       NA's   :10       NA's   :11       NA's   :11      
     REND_7           REND_8           REND_9        TON_PROD_1  
 Min.   :0.5500   Min.   :0.5900   Min.   :0.870   Min.   :   4  
 1st Qu.:0.7250   1st Qu.:0.7800   1st Qu.:1.010   1st Qu.: 405  
 Median :0.8100   Median :0.8900   Median :1.070   Median : 963  
 Mean   :0.7548   Mean   :0.8258   Mean   :1.075   Mean   :1472  
 3rd Qu.:0.8200   3rd Qu.:0.8900   3rd Qu.:1.150   3rd Qu.:1748  
 Max.   :0.8600   Max.   :1.1300   Max.   :1.220   Max.   :7306  
 NA's   :11       NA's   :11       NA's   :11      NA's   :7     
  TON_PROD_10     TON_PROD_11     TON_PROD_12       TON_PROD_2    
 Min.   :  270   Min.   :  187   Min.   : 228.0   Min.   :   0.0  
 1st Qu.: 1358   1st Qu.:  823   1st Qu.: 821.8   1st Qu.: 394.5  
 Median : 1864   Median : 2128   Median :2245.5   Median : 876.0  
 Mean   : 2913   Mean   : 3096   Mean   :2967.8   Mean   :1394.6  
 3rd Qu.: 3854   3rd Qu.: 4161   3rd Qu.:3994.5   3rd Qu.:1757.5  
 Max.   :11571   Max.   :11315   Max.   :9500.0   Max.   :6176.0  
 NA's   :11      NA's   :11      NA's   :12       NA's   :7       
   TON_PROD_3     TON_PROD_4       TON_PROD_5       TON_PROD_6  
 Min.   :  36   Min.   :  35.0   Min.   :   5.0   Min.   :  88  
 1st Qu.: 445   1st Qu.: 586.2   1st Qu.: 555.5   1st Qu.: 560  
 Median :1030   Median : 840.5   Median : 768.0   Median :1260  
 Mean   :1529   Mean   :1465.9   Mean   :1337.8   Mean   :1775  
 3rd Qu.:1917   3rd Qu.:1899.2   3rd Qu.:1592.5   3rd Qu.:2456  
 Max.   :6323   Max.   :6074.0   Max.   :5599.0   Max.   :7238  
 NA's   :10     NA's   :10       NA's   :11       NA's   :11    
   TON_PROD_7       TON_PROD_8     TON_PROD_9       YEAR_1    
 Min.   : 135.0   Min.   : 148   Min.   : 223   Min.   :2007  
 1st Qu.: 752.5   1st Qu.: 886   1st Qu.:1266   1st Qu.:2007  
 Median :1289.0   Median :1530   Median :1783   Median :2007  
 Mean   :1834.3   Mean   :2133   Mean   :2711   Mean   :2007  
 3rd Qu.:2417.5   3rd Qu.:2712   3rd Qu.:3807   3rd Qu.:2007  
 Max.   :5726.0   Max.   :8487   Max.   :8899   Max.   :2011  
 NA's   :11       NA's   :11     NA's   :11     NA's   :7     
    YEAR_10        YEAR_11        YEAR_12         YEAR_2    
 Min.   :2016   Min.   :2017   Min.   :2018   Min.   :2008  
 1st Qu.:2016   1st Qu.:2017   1st Qu.:2018   1st Qu.:2008  
 Median :2016   Median :2017   Median :2018   Median :2008  
 Mean   :2016   Mean   :2017   Mean   :2018   Mean   :2008  
 3rd Qu.:2016   3rd Qu.:2017   3rd Qu.:2018   3rd Qu.:2008  
 Max.   :2017   Max.   :2018   Max.   :2018   Max.   :2012  
 NA's   :11     NA's   :11     NA's   :12     NA's   :7     
     YEAR_3         YEAR_4         YEAR_5         YEAR_6    
 Min.   :2009   Min.   :2010   Min.   :2011   Min.   :2012  
 1st Qu.:2009   1st Qu.:2010   1st Qu.:2011   1st Qu.:2012  
 Median :2009   Median :2010   Median :2011   Median :2012  
 Mean   :2009   Mean   :2010   Mean   :2011   Mean   :2012  
 3rd Qu.:2009   3rd Qu.:2010   3rd Qu.:2011   3rd Qu.:2012  
 Max.   :2017   Max.   :2018   Max.   :2011   Max.   :2013  
 NA's   :10     NA's   :10     NA's   :11     NA's   :11    
     YEAR_7         YEAR_8         YEAR_9              geometry 
 Min.   :2013   Min.   :2014   Min.   :2015   POLYGON      :42  
 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.   :2014   Max.   :2015   Max.   :2016                     
 NA's   :11     NA's   :11     NA's   :11                       

5. Trazado

Para realizar el trazado necesitaremos instalar RColorBrewer:

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

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

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

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

  mapa <- leaflet(data = cau_munic_stat) %>%
  addTiles() %>%
  addPolygons(label = ~TON_PROD_12, 
              popup = ~MPIO_CNMBR,
              fillColor = ~pal(TON_PROD_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 = ~TON_PROD_12,
    title = "Producción de café en Cauca [Ton] (2018)",
    opacity = 1
  )

mapa  

sessionInfo() 

R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Mexico.1252  LC_CTYPE=Spanish_Mexico.1252   
[3] LC_MONETARY=Spanish_Mexico.1252 LC_NUMERIC=C                   
[5] LC_TIME=Spanish_Mexico.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.1.0      
 [4] ggrepel_0.8.2       GSODR_2.1.2         rnaturalearth_0.1.0
 [7] viridis_0.5.1       viridisLite_0.3.0   sf_0.9-5           
[10] maptools_1.0-2      rgeos_0.5-5