Estadisticas municipales de la agricultura en el Valle del Cuaca

1. ¿Por qué este Notebook?

Este es un R Markdown Notebook que ilustra estadísticas agrícolas para el departamento del Valle del Cauca en Colombia. El presente cuaderno de R ayuda a los estudiantes de la Universidad Nacional de Colombia que se encuentran cursando el curso de Geomática Básica a entender conceptos básicos de geomática aplicables a la agronomía.

2. Funcionalidades de GIS.

La exploración de estadísticas no espaciales es esencial para comprender que ocurre en los territorios. Muchas librerias de R, particularmente tidyverse y dplyr son muy útiles para explorar y resumir estadísticas. Por otro lado, las operaciones geoespaciales pueden mejorar nuestro entendimiento acerca de múltiples problemas que afectan las regoiones geográficas. Un buen ejemplo es cuando se desea averiguar la ubicación de ciertos municipios que tienen rendimientos de cosecha óptimos. Para realizar esta exploración se necesita unir datos no espaciales con datos espaciales. Esta labor se realizó anteriormente en QGIS, ahora se utilizará R para llevarla acabo. Adicionalmente, se puede explorar uniones espaciales. Estas opciones se basan en la interaccion de dos objetos espaciales, a menudo polígonos y puntos. Hay muchas formas de unir objetos, las cuales pueden incluir opciones específicas como: cerca, dentro, toques, etc. Comencemos removiendo el contenido de la memoria:

rm(list=ls())

Ahora, instalemos las librerias necesarias. En el siguiente chunk, se indica que solo sean instalados los paquetes que no han sido previamente 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)

Ahora, carguemos las librerias.

library(here)

library(tidyverse)

library(rgeos)

library(maptools)

library(raster)

library(sf)

library(viridis)

library(rnaturalearth)
library(GSODR)
library(ggrepel)
library(cowplot)

3. Exploración de estadísticas agrícolas en el Valle del Cauca

Previamente he descargado datos estadísticos en formato csv de la página oficial del Ministerio de Agricultura y Desarrollo Rural en mi computador. Después, en Excel, eliminé información no correspondiente al departamento del Valle del Cacua. Guardé el archivo unicamente con información correspondiente a mi departamento en un documento llamado EVA_Valle2.csv, luego lo cargué dentro de una carpeta llamada agro en RStudio Cloud. Leamos el archivo csv con “Estadísticas municipales agropecuarias” para el Valle del Cauca

datos <- read_csv("./agro/EVA_Valle2.csv")

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(),
  YEAR = col_double(),
  Area_Siembra = col_double(),
  Area_Cosecha = col_double(),
  Produccion = col_double(),
  Rendimiento = col_double(),
  ESTADO = col_character(),
  CICLO = col_character()
)

Chequeemos cuales son los atributos de datos: tanto al principio (head) como al final (tail)

head(datos)

tail(datos)

Se recalca que cada municipio cuenta con estadísticas acerca del área sembrada, área cosechada y rendimiento de los cultivos que allí se producen en diferentes años. En la tabla no hay unidades, sin embargo, si se verifica el archivio csv original se encuentra que las unidades para el área son las hectáreas y para el rendimiento son Ton/ha. Usaremos la libreria (dplyr) para explorar los contenidos de los datos del objeto. Para esto, primero obtenemos un sumario de rendimiento por el 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)

De igual forma, podemos calcular el rendimiento promedio por grupo en los municipios del Valle del Cauca:

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

rend_valledelcauca

Notemos que los rendimientos más altos corresponden a CEREALES, FIBRAS Y FLORES Y FORRAJES. Ahora veamos cuales son los municipios con mayor rendimiento para cada grupo de cultivos 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

En seguida miraremos cuales 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

Notemos que el áremas más grande cosechada en el año 2018 para Jamundí corresponde a Cereales, pues esta zona, ubicada en el sur del departamento, se siembran grandes extensiones de arroz. Seleccionemos la producción de arroz (toneladas) en Jamundí para cada año:

datos %>% 
  filter(MUNICIPIO=="JAMUNDI" & SUBGRUPO=="ARROZ") %>% 
  group_by(YEAR, CULTIVO) ->  jamundi_arroz

jamundi_arroz

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

g <- ggplot(aes(x=YEAR, y=Produccion/1000), data = jamundi_arroz) + geom_bar(stat='identity') + labs(y='Producción de arroz [Ton x 1000]')

g + ggtitle("Evolución de la produccion de arroz en Jamundi desde 2007 al 2018") + labs(caption= "Basado en datos del DANE, 2018")

Ahora miraremos que cultivos tuvieron la mayor área cosechada en el año 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

Observamos que otros cultivos permanentes corresponden a la mayor área cosechada en el 2018 en el Valle del Cauca. Para observar que cultivos pertenecen a este grupoo, puede buscar información en la página oficial del DANE, entidad responsable de la planeación, levantamiento, procesamiento, análisis y difusión de las estadísticas oficiales de Colombia. De igual forma, podemos buscar esta información en los mismos datos.

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

Resaltamos que la caña de azucar es de total importancia en la agricultura vallecaucana debido a la gran cantidad de ingenios azucareros que se encuentran en esta región. A continuación veamos cuales 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

Regresemos a los datos de cultivos permanentes. Antes de plotear, necesitamos agregar al atributo total_area_cosecha, un nuevo campo con abreviaturas para cada GRUPO de cultivos para que el plot no luzca desordenado.

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

Ahora es tiempo de 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 el 2018 para el Valle") + theme(plot.title = element_text(hjust = 0.5)) +
   labs(caption= "Basado en datos del DANE, 2018")

Se puede obtener más información correspondiente al departamento, especificamente de diferentes cultivos que son importantes principalemente para la economía de este.

4. Uniendo estadísticas agrícolas a municipios

Previamente he subido a RStudio Cloud los datos administrativos del Valle del Cauca. Empecemos leyendo los datos usando la libreria sf

ant_munic <- sf::st_read("./valledelcauca/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Reading layer `MGN_MPIO_POLITICO' from data source `/cloud/project/valledelcauca/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
Simple feature collection with 42 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -77.54977 ymin: 3.091239 xmax: -75.70724 ymax: 5.047394
CRS:            4326

¿Qué hay en ant_munic?

ant_munic

Simple feature collection with 42 features and 9 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -77.54977 ymin: 3.091239 xmax: -75.70724 ymax: 5.047394
CRS:            4326
First 10 features:
   DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR                           MPIO_CRSLC MPIO_NAREA
1          76      76001         CALI                                 1536  563.04276
2          76      76036    ANDALUCÍA     Ordenanza 38 de Abril 25 de 1921  110.46042
3          76      76041 ANSERMANUEVO                 Ordenanza 29 de 1925  305.45118
4          76      76054      ARGELIA                 Ordenanza 15 de 1956   90.79604
5          76      76111         BUGA                                 1555  825.86513
6          76      76113 BUGALAGRANDE                                 1791  396.78132
7          76      76122   CAICEDONIA Ordenanza 21 del 20 de Abril de 1923  166.98369
8          76      76126       CALIMA     Ordenanza 49 de Junio 23 de 1939  793.49323
9          76      76130   CANDELARIA                                 1797  296.46056
10         76      76147      CARTAGO                                 1863  248.16005
   MPIO_NANO      DPTO_CNMBR Shape_Leng  Shape_Area                       geometry
1       2017 VALLE DEL CAUCA  1.1636697 0.045733211 MULTIPOLYGON (((-76.59175 3...
2       2017 VALLE DEL CAUCA  0.6651004 0.008984851 MULTIPOLYGON (((-76.22406 4...
3       2017 VALLE DEL CAUCA  0.9460442 0.024871077 MULTIPOLYGON (((-76.01558 4...
4       2017 VALLE DEL CAUCA  0.4538261 0.007391017 MULTIPOLYGON (((-76.14316 4...
5       2017 VALLE DEL CAUCA  2.0063716 0.067157337 MULTIPOLYGON (((-76.31608 3...
6       2017 VALLE DEL CAUCA  1.0336063 0.032279294 MULTIPOLYGON (((-76.15131 4...
7       2017 VALLE DEL CAUCA  0.6465900 0.013590500 MULTIPOLYGON (((-75.8539 4....
8       2017 VALLE DEL CAUCA  1.5441977 0.064484981 MULTIPOLYGON (((-76.51747 4...
9       2017 VALLE DEL CAUCA  0.8709866 0.024086066 MULTIPOLYGON (((-76.30455 3...
10      2017 VALLE DEL CAUCA  0.8955580 0.020211481 MULTIPOLYGON (((-75.94518 4...

Debemos tener en cuenta que ant_munic es una colección de características simples. Tambien debemos tener en cuenta que los datos utilizan un sistema de referencia de coordenas geográficas WGS1984 (4326 epsg código ) 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) para poder hacer la unión. El mejor atributo es una identificación.

datos %>% filter (MUNICIPIO =="CALI") ->  cali_datos

cali_datos

class(cali_datos$COD_MUN)

[1] "numeric"

Para poder realizar la unión, necesitamos cambiar tanto el tipo de datos como el contenido del código que identifica a cada municipio. Para esta labor, es una buena opción crear una copia de los datos estadísticos originales con el fin de cualquier movimiento en falso no perjudique la información original. Avancemos paso a paso:

datos2 <- datos

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

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

head(datos2)

Asegurarse de verificar, en el objeto de datos 2, las caracteristicas del nuevo atributo MPIO_CCDGO Ahora, filtremos un solo año y seleccionemos 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(MUNICIPIO, MPIO_CCDGO, YEAR, Produccion, Rendimiento)

datos4

datos4 %>% 
  gather("YEAR", "Produccion", "Rendimiento" , key = variable, value = number)

head(datos4)

Esta es una tarea clave. Implica varios pasos para poder convertir la tabla de atributos de formato largo a ancho.

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

head(datos5)

tail(datos5)

Haremos una copia de la colección de características simples (en caso de un movimiento falso)

ant_munic2 <- ant_munic

Ahora es tiempo de unir:

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        MPIO_NANO   
 76:42      76001  : 1   ALCALÁ      : 1   1536   : 2   Min.   :  41.86   Min.   :2017  
            76020  : 1   ANDALUCÍA   : 1   1824   : 2   1st Qu.: 191.01   1st Qu.:2017  
            76036  : 1   ANSERMANUEVO: 1   1864   : 2   Median : 266.71   Median :2017  
            76041  : 1   ARGELIA     : 1   1555   : 1   Mean   : 492.04   Mean   :2017  
            76054  : 1   BOLÍVAR     : 1   1576   : 1   3rd Qu.: 431.81   3rd Qu.:2017  
            76100  : 1   BUENAVENTURA: 1   1630   : 1   Max.   :6292.50   Max.   :2017  
            (Other):36   (Other)     :36   (Other):33                                   
           DPTO_CNMBR   Shape_Leng       Shape_Area        MUNICIPIO          Produccion_1 
 VALLE DEL CAUCA:42   Min.   :0.4538   Min.   :0.003409   Length:42          Min.   :  76  
                      1st Qu.:0.7223   1st Qu.:0.015552   Class :character   1st Qu.: 405  
                      Median :0.9059   Median :0.021696   Mode  :character   Median : 800  
                      Mean   :1.1602   Mean   :0.039988                      Mean   :1785  
                      3rd Qu.:1.3207   3rd Qu.:0.035099                      3rd Qu.:1910  
                      Max.   :6.5953   Max.   :0.510779                      Max.   :8284  
                                                                             NA's   :3     
 Produccion_10  Produccion_11    Produccion_12     Produccion_2     Produccion_3 
 Min.   : 162   Min.   : 142.0   Min.   : 139.0   Min.   :  85.0   Min.   :  83  
 1st Qu.: 472   1st Qu.: 410.5   1st Qu.: 398.5   1st Qu.: 369.5   1st Qu.: 340  
 Median : 745   Median : 780.0   Median : 748.0   Median : 785.0   Median : 872  
 Mean   :1463   Mean   :1325.3   Mean   :1273.5   Mean   :1683.8   Mean   :1608  
 3rd Qu.:1550   3rd Qu.:1657.0   3rd Qu.:1734.5   3rd Qu.:1960.5   3rd Qu.:1553  
 Max.   :5895   Max.   :6212.0   Max.   :5906.0   Max.   :8284.0   Max.   :8764  
 NA's   :3      NA's   :3        NA's   :3        NA's   :3        NA's   :3     
  Produccion_4    Produccion_5     Produccion_6     Produccion_7   Produccion_8 
 Min.   :   90   Min.   :  86.0   Min.   : 100.0   Min.   :  75   Min.   : 100  
 1st Qu.:  414   1st Qu.: 420.5   1st Qu.: 415.5   1st Qu.: 287   1st Qu.: 351  
 Median :  861   Median : 992.0   Median : 776.0   Median : 528   Median : 670  
 Mean   : 1782   Mean   :1678.9   Mean   :1569.0   Mean   :1101   Mean   :1277  
 3rd Qu.: 1913   3rd Qu.:1908.0   3rd Qu.:1812.0   3rd Qu.:1019   3rd Qu.:1200  
 Max.   :10154   Max.   :8377.0   Max.   :8587.0   Max.   :5161   Max.   :5757  
 NA's   :3       NA's   :3        NA's   :3        NA's   :3      NA's   :3     
  Produccion_9  Rendimiento_1    Rendimiento_10  Rendimiento_11   Rendimiento_12  
 Min.   : 158   Min.   :  1.00   Min.   :  1.0   Min.   :  9.00   Min.   : 62.00  
 1st Qu.: 461   1st Qu.:  8.00   1st Qu.: 96.0   1st Qu.: 67.00   1st Qu.: 78.00  
 Median : 772   Median : 13.00   Median :104.0   Median : 84.00   Median : 86.00  
 Mean   :1476   Mean   : 41.31   Mean   : 93.1   Mean   : 90.59   Mean   : 99.46  
 3rd Qu.:1544   3rd Qu.: 76.00   3rd Qu.:107.5   3rd Qu.:116.00   3rd Qu.:119.00  
 Max.   :6108   Max.   :119.00   Max.   :127.0   Max.   :182.00   Max.   :187.00  
 NA's   :3      NA's   :3        NA's   :3       NA's   :3        NA's   :3       
 Rendimiento_2    Rendimiento_3    Rendimiento_4    Rendimiento_5    Rendimiento_6   
 Min.   :  1.00   Min.   :  1.00   Min.   :  1.00   Min.   :  1.00   Min.   :  1.00  
 1st Qu.:  9.00   1st Qu.:  6.50   1st Qu.:  7.50   1st Qu.:  7.50   1st Qu.:  5.00  
 Median : 44.00   Median : 33.00   Median : 13.00   Median : 12.00   Median :  8.00  
 Mean   : 48.03   Mean   : 47.51   Mean   : 43.54   Mean   : 43.28   Mean   : 25.41  
 3rd Qu.: 81.00   3rd Qu.: 86.00   3rd Qu.: 85.50   3rd Qu.: 86.50   3rd Qu.: 54.00  
 Max.   :143.00   Max.   :141.00   Max.   :141.00   Max.   :141.00   Max.   :102.00  
 NA's   :3        NA's   :3        NA's   :3        NA's   :3        NA's   :3       
 Rendimiento_7    Rendimiento_8    Rendimiento_9        YEAR_1        YEAR_10    
 Min.   :  4.00   Min.   :  7.00   Min.   :  1.00   Min.   :2007   Min.   :2016  
 1st Qu.: 48.00   1st Qu.: 44.50   1st Qu.: 94.00   1st Qu.:2007   1st Qu.:2016  
 Median : 64.00   Median : 78.00   Median :102.00   Median :2007   Median :2016  
 Mean   : 61.08   Mean   : 65.95   Mean   : 92.23   Mean   :2007   Mean   :2016  
 3rd Qu.: 72.00   3rd Qu.: 88.00   3rd Qu.:105.50   3rd Qu.:2007   3rd Qu.:2016  
 Max.   :121.00   Max.   :131.00   Max.   :123.00   Max.   :2007   Max.   :2016  
 NA's   :3        NA's   :3        NA's   :3        NA's   :3      NA's   :3     
    YEAR_11        YEAR_12         YEAR_2         YEAR_3         YEAR_4         YEAR_5    
 Min.   :2017   Min.   :2018   Min.   :2008   Min.   :2009   Min.   :2010   Min.   :2011  
 1st Qu.:2017   1st Qu.:2018   1st Qu.:2008   1st Qu.:2009   1st Qu.:2010   1st Qu.:2011  
 Median :2017   Median :2018   Median :2008   Median :2009   Median :2010   Median :2011  
 Mean   :2017   Mean   :2018   Mean   :2008   Mean   :2009   Mean   :2010   Mean   :2011  
 3rd Qu.:2017   3rd Qu.:2018   3rd Qu.:2008   3rd Qu.:2009   3rd Qu.:2010   3rd Qu.:2011  
 Max.   :2017   Max.   :2018   Max.   :2008   Max.   :2009   Max.   :2010   Max.   :2011  
 NA's   :3      NA's   :3      NA's   :3      NA's   :3      NA's   :3      NA's   :3     
     YEAR_6         YEAR_7         YEAR_8         YEAR_9              geometry 
 Min.   :2012   Min.   :2013   Min.   :2014   Min.   :2015   MULTIPOLYGON :42  
 1st Qu.:2012   1st Qu.:2013   1st Qu.:2014   1st Qu.:2015   epsg:4326    : 0  
 Median :2012   Median :2013   Median :2014   Median :2015   +proj=long...: 0  
 Mean   :2012   Mean   :2013   Mean   :2014   Mean   :2015                     
 3rd Qu.:2012   3rd Qu.:2013   3rd Qu.:2014   3rd Qu.:2015                     
 Max.   :2012   Max.   :2013   Max.   :2014   Max.   :2015                     
 NA's   :3      NA's   :3      NA's   :3      NA's   :3                        

5. Ploteo

Se hará un plot para 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$PProduccion_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 Café en el Valle [Ton] (2018)",
    opacity = 1
  )

mapa