1. Introducción

Con el fin de identificar e ilustrar las estadisticas agricolas del departamento de Nariño en Colombia durante el año 2018 se a compilado el siguiente cuaderno.

2. Funcionalidades de GIS

Se hara uso de las librerias de R que permiten explorar y resumir estadisticas, estas librerias son: dplyr y tidyverse. Por otra parte, operaciones geoespaciales pueden mejorar el entendiendo de problemas que afectan las regiones geográficas. Para esto es necesario unir datos espaciales con no espaciales o realizar uniones espaciales. Estas operaciones están basadas en la intersección entre dos objetos especiales, usualmente puntos o polígonos. Existen muchas formas en las que se pueden unir los objetos, las cuales incluyen opciones específicas.

El primer paso es el de llamar a las librerias que se utilizarán:

library(here)
## here() starts at C:/Users/USER/Documents
library(tidyverse)
## -- 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 dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(rgeos)
## Loading required package: sp
## 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)
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:tidyr':
## 
##     extract
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)
library(ggrepel)
library(cowplot)

3. Explorando las estadisticas agrícolas en Nariño

Se realiza la descarga de los datos en la pagina oficial de Evaluaciones Agropecuarias Municipales EVA, este archivo contiene informacion de la agricultura del departamento de nariño, es por ello que se realizan las modificaciones necesarias par obtener la informacion deseada. Lo siguiente es establecer este archivo en el objeto datos:

setwd("C:/Users/USER/Desktop/Geomatica")
datos <- read_csv2("C:/Users/USER/Desktop/Geomatica/EVA_NARINO1.csv")
## Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
## Parsed with column specification:
## cols(
##   COD_DEPTO = col_double(),
##   DEPTO = col_character(),
##   COD_MUN = col_double(),
##   MUN = col_character(),
##   GRUPO = col_character(),
##   SUBGRUPO = col_character(),
##   CULTIVO = col_character(),
##   ANO = col_double(),
##   PERIODO = col_character(),
##   AREA_SEMBRADA = col_double(),
##   AREA_COSECHADA = col_double(),
##   PRODUCCION = col_double(),
##   RENDIMIENTO = col_double(),
##   ESTADO = col_character(),
##   N_CIENTIFICO = col_character(),
##   `CICLO,` = col_character()
## )
head(datos)
tail(datos)

Los datos suministrados presentan como atributos la superficie sembrada, superficie cosechada y el rendimiento de los diferentes cultivos registrados en el departamento de Nariño.

Cabe mencionar que las unidades de dichos atributos son hectareas (ha) y en el caso de el rendimiento son toneladas por hectarea (Ton/ha).

El siguente paso es hallar el promedio del rendimiento en los diferentes años (2006-2018) por grupo y municipio mediante el uso de la biblioteca dplyr el cual provee las herramientas correspondientes. Para ello el rendimiento se obtendra de acuerdo al grupo que pertenece el cultivo y al municipio donde se localiza.

datos %>%
  group_by(MUN, GRUPO) %>%
  summarise(rend_prom = mean(RENDIMIENTO, na.rm = TRUE)) -> rend_resumen
## `summarise()` regrouping output by 'MUN' (override with `.groups` argument)
head(rend_resumen)

Sumado a lo anterior cabe la posibilidad de calcular el rendimiento promedio por cada grupo de cultivos en el municipio de Nariño.

datos %>%
  group_by(GRUPO) %>%
  summarise(rend_dep = mean(RENDIMIENTO, na.rm = TRUE)) -> rend_narino
## `summarise()` ungrouping output (override with `.groups` argument)
rend_narino

De la anterior tabla caba recalcar que los mayores rendimientos corresponden a los cultivos de Hortalizas, seguido en menor medida de los frutales.

Ahora se determina los municipios que presenten un mayor rendimiento en cada uno de los grupos durante el año 2018.

datos %>% 
  filter(ANO==2018) %>% 
  group_by(GRUPO, MUN) %>%
  summarize(max_rend = max(RENDIMIENTO, na.rm = TRUE)) %>%
    slice(which.max(max_rend)) -> rend_max_18
## `summarise()` regrouping output by 'GRUPO' (override with `.groups` argument)
rend_max_18

y para complementar se determina cuales son los municipios con mayor área cosechada para cada grupo de cultivos en 2018.

datos %>% 
  filter(ANO==2018) %>% 
  group_by(GRUPO, MUN) %>%
  summarize(max_area_cosecha = max(AREA_COSECHADA, na.rm = TRUE)) %>%
    slice(which.max(max_area_cosecha)) -> area_cosecha_max
## `summarise()` regrouping output by 'GRUPO' (override with `.groups` argument)
area_cosecha_max

De acuerdo a lo determinado anteriormente, nos encontramos que el municipio de GUAITARILLA es el que presenta una mayor area de cosecha en los cultivos de CEREALES. Esto se debe a que este municipio dedica una gran extension de area para la produccion de Trigo de acuerdo a los datos obtenidos de Agronet (2018).

Si se desea ser mas especifico, frente al cultivo, es posible obtener la informacion anterior pero enfocada en el Subgrupo:

datos %>% 
  filter(MUN=="GUAITARILLA" & SUBGRUPO=="TRIGO") %>% 
  group_by(ANO, CULTIVO) ->  gua_tri

gua_tri

Sin embargo, en busca de una mejor interpretacion de los datos se realiza una representacion mediante un grafico de barras en el cual sea posible observar el comportamiento del cultivo de Trigo duarante los años de los datos recolectados en el municipio de Guaitarilla.

g <- ggplot(aes(x=ANO, y=PRODUCCION/1000), data = gua_tri) + geom_bar(stat='identity') + labs(y='Producción de Trigo en Guaitarilla [Ton x 1000]', x='AÑO')
g + ggtitle("Evolución de la producción de Trigo en Guaitarilla 2007-2018") + labs(caption= "Fuente de datos: EMA data (DANE, 2018)")

Como muestra la grafica, el cultivo de Trigo ha presentado una disminucion drastica en terminos de produccion, es por ello que es interesante realizar el analisis de que cultivos presentaron la mayor area de cosecha durante el 2018:

datos %>% 
  filter(ANO==2018) %>% 
  group_by(GRUPO) %>%
  summarize(sum_area_cosecha = sum(AREA_COSECHADA, na.rm = TRUE)) %>%
     arrange(desc(sum_area_cosecha)) -> total_area_cosecha
## `summarise()` ungrouping output (override with `.groups` argument)
total_area_cosecha

En donde se obtiene que los tuberculos son el cultivo con mayor area de cosecha durante el 2018 despues del grupo de otros permanentes en donde podemos encontrarnos con cultivos como el cacao o el café.

datos %>%
  filter(GRUPO=="OTROS PERMANENTES" & ANO==2018) %>%
  group_by(CULTIVO) %>%
  summarize(sum_cosecha = sum(AREA_COSECHADA, na.rm = TRUE)) %>%
     arrange(desc(sum_cosecha)) -> total_cosecha
## `summarise()` ungrouping output (override with `.groups` argument)
total_cosecha

Lo anteriormente mencionado se ve representado en esta ultima tabla, donde se puede observar que el cafe es el cultivo con mayor area de cosecha en el 2018 seguido de la caña panelera y el cacao.

El siguiente paso es identificar los municipios con mayor área cosechada para cada cultivo permanente en el 2018.

datos %>% 
  filter(ANO==2018 & GRUPO=="OTROS PERMANENTES") %>% 
  group_by(CULTIVO, MUN) %>%
  summarize(max_area2 = max(AREA_COSECHADA, na.rm = TRUE)) %>%
    slice(which.max(max_area2)) -> area_cosecha2
## `summarise()` regrouping output by 'CULTIVO' (override with `.groups` argument)
area_cosecha2

Para generar un mayor orden en la observacion de los datos se crea el siguiente objeto:

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

Y se raliza el plot correspondiente:

g <- ggplot(aes(x=GRUPO, y=sum_area_cosecha), data = total_area_cosecha) + geom_bar(stat='identity') + labs(y='Total Area cosechada[Ha]')
g+ ggtitle("Area total cosechada de Nariño por grupos en 2018.") + theme(plot.title = element_text(hjust = 0.5)) +
   labs(caption= "Fuente: EMA data (DANE, 2018)")

4. Incorporacion de las estadisticas agricolas a los municipios

En esta sección se tiene como requerimiento el Shapefile correspondiente al departemento de Nariño, los cuales se obtubieron del DANE (2020).

nar_munic <- sf::st_read("/Users/USER/Desktop/Geomatica/52_NARIÑO/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\USER\Desktop\Geomatica\52_NARIÑO\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 64 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -79.01021 ymin: 0.3613481 xmax: -76.83368 ymax: 2.683898
## geographic CRS: WGS 84

Es importante identificar las propiedades del archivo que se va a utilizar así como el sistema de coordenadas que se esta utilizando:

#nar_munic

Para poder realizar el enlaze entre los dos objetos es necesario un atributo común, en este caso se usara el código del municiipio. Para los datos con los que se a trabajado durante el informe el atributo corresponde a COD_MUN, en cambio para el .shp cuyo objeto se le determino como nar_munic el atributo es MPIO_CCDGO.

datos %>% filter (MUN =="PASTO") ->  pas_datos
pas_datos

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 esto se crea una copia de los datos estadísticos originales.

datos2 <- datos
datos2$TEMP <-  as.character(datos2$COD_MUN)
datos2$MPIO_CCDGO<- as.factor(datos2$TEMP)

Con lo anterior se ha realizado la union de los datos, se puede verificar:

head(datos2)

El siguiente proceso tiene como fin el filtrar solo un año y seleccionar solo dos atributos relevantes, en este caso se escoge el cafe debido a su gran importancia en el año 2018:

datos2 %>% filter(CULTIVO == "CAFE" & ANO ==2018)  -> datos3
head(datos3)
class(datos3)
## [1] "spec_tbl_df" "tbl_df"      "tbl"         "data.frame"
datos4 <- datos3 %>% dplyr::select(MUN, MPIO_CCDGO, ANO, PRODUCCION, RENDIMIENTO) 
datos4
datos4 %>% 
  gather("ANO", "PRODUCCION", "RENDIMIENTO" , key = variable, value = number)
head(datos4)

Ahora se deben realizar varios pasos con el fin de convertir la tabla de atributos de formato largo a formato ancho.

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

Se realiza una copia de los datos originales por si se da la situacion de realizar algún cambio:

nar_munic2 <- nar_munic
nar_munic_stat = left_join(nar_munic2, datos5, by="MPIO_CCDGO")
summary(nar_munic_stat)
##   DPTO_CCDGO         MPIO_CCDGO         MPIO_CNMBR         MPIO_CRSLC       
##  Length:64          Length:64          Length:64          Length:64         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##    MPIO_NAREA        MPIO_NANO     DPTO_CNMBR          Shape_Leng    
##  Min.   :  25.31   Min.   :2017   Length:64          Min.   :0.2409  
##  1st Qu.:  96.64   1st Qu.:2017   Class :character   1st Qu.:0.4636  
##  Median : 236.75   Median :2017   Mode  :character   Median :0.7825  
##  Mean   : 492.15   Mean   :2017                      Mean   :1.1660  
##  3rd Qu.: 547.73   3rd Qu.:2017                      3rd Qu.:1.4100  
##  Max.   :3611.35   Max.   :2017                      Max.   :6.4534  
##                                                                      
##    Shape_Area           MUN                ANO_1       PRODUCCION_1   
##  Min.   :0.002050   Length:64          Min.   :2018   Min.   :   3.0  
##  1st Qu.:0.007825   Class :character   1st Qu.:2018   1st Qu.: 152.0  
##  Median :0.019186   Mode  :character   Median :2018   Median : 430.0  
##  Mean   :0.039820                      Mean   :2018   Mean   : 839.3  
##  3rd Qu.:0.044323                      3rd Qu.:2018   3rd Qu.:1089.0  
##  Max.   :0.291654                      Max.   :2018   Max.   :4393.0  
##                                        NA's   :31     NA's   :31      
##  RENDIMIENTO_1            geometry 
##  Min.   :0.620   MULTIPOLYGON :64  
##  1st Qu.:0.860   epsg:4326    : 0  
##  Median :1.030   +proj=long...: 0  
##  Mean   :1.032                     
##  3rd Qu.:1.120                     
##  Max.   :1.870                     
##  NA's   :31

5. Plotting

library(RColorBrewer)
library(leaflet)

Se grafican los municipios con su correspondiente producción de café para el año 2018, para esto se hara uso de una escala de color en donde se identifique la diferencia de produccion entre municipios:

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

  mapa <- leaflet(data = nar_munic_stat) %>%
  addTiles() %>%
  addPolygons(label = ~PRODUCCION_1,
              popup = ~MPIO_CNMBR,
              fillColor = ~pal(PRODUCCION_1),
              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_1,
    title = "Produccion de Cafe en Nariño [Ton] (2018)",
    opacity = 1
  )
mapa

Si se desea hacer una comparacion con el cultivo de Papa el cual es el segundo mas importante de acuerdo a agronet (2017). Solo se debe repetir los pasos del punto 4 y 5 con el fin de construir otro mapa.

datos %>% 
  filter(ANO==2018 & GRUPO=="TUBERCULOS Y PLATANOS") %>% 
  group_by(CULTIVO, MUN) %>%
  summarize(max_rend = max(RENDIMIENTO, na.rm = TRUE)) %>%
    slice(which.max(max_rend)) -> max_rend2
## `summarise()` regrouping output by 'CULTIVO' (override with `.groups` argument)
max_rend2
datos6<-datos
datos6$TEMP <-  as.character(datos6$COD_MUN)
datos6$MPIO_CCDGO<- as.factor(datos6$TEMP)
head(datos6)
datos6 <- datos2 %>% filter(CULTIVO== "PAPA", ANO==2018)
datos6
datos7 <- datos6 %>% dplyr::select(MUN, MPIO_CCDGO, ANO, PRODUCCION, RENDIMIENTO)
datos7
datos7 %>% 
  gather("ANO", "PRODUCCION", "RENDIMIENTO" , key = variable, value = number)
datos7 %>% 
  group_by(MPIO_CCDGO) %>% 
  mutate(Visit = 1:n()) %>% 
  gather("ANO", "PRODUCCION", "RENDIMIENTO", key = variable, value = number) %>% 
  unite(combi, variable, Visit) %>%
  spread(combi, number) -> datos8

head(datos8)
tail(datos8)
nar_munic3 <- nar_munic
nar_munic_stat2 = left_join(nar_munic3, datos8, by="MPIO_CCDGO")
summary(nar_munic_stat2)
##   DPTO_CCDGO         MPIO_CCDGO         MPIO_CNMBR         MPIO_CRSLC       
##  Length:64          Length:64          Length:64          Length:64         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##    MPIO_NAREA        MPIO_NANO     DPTO_CNMBR          Shape_Leng    
##  Min.   :  25.31   Min.   :2017   Length:64          Min.   :0.2409  
##  1st Qu.:  96.64   1st Qu.:2017   Class :character   1st Qu.:0.4636  
##  Median : 236.75   Median :2017   Mode  :character   Median :0.7825  
##  Mean   : 492.15   Mean   :2017                      Mean   :1.1660  
##  3rd Qu.: 547.73   3rd Qu.:2017                      3rd Qu.:1.4100  
##  Max.   :3611.35   Max.   :2017                      Max.   :6.4534  
##                                                                      
##    Shape_Area           MUN                ANO_1       PRODUCCION_1   
##  Min.   :0.002050   Length:64          Min.   :2018   Min.   :   6.0  
##  1st Qu.:0.007825   Class :character   1st Qu.:2018   1st Qu.:  39.0  
##  Median :0.019186   Mode  :character   Median :2018   Median : 101.5  
##  Mean   :0.039820                      Mean   :2018   Mean   : 666.7  
##  3rd Qu.:0.044323                      3rd Qu.:2018   3rd Qu.: 540.0  
##  Max.   :0.291654                      Max.   :2018   Max.   :6552.0  
##                                        NA's   :36     NA's   :36      
##  RENDIMIENTO_1            geometry 
##  Min.   : 1.00   MULTIPOLYGON :64  
##  1st Qu.:11.75   epsg:4326    : 0  
##  Median :18.00   +proj=long...: 0  
##  Mean   :17.36                     
##  3rd Qu.:20.50                     
##  Max.   :45.00                     
##  NA's   :36
bins <- c(0, 250, 500, 1000, 2000, 5000, 10000, 15000)
pal <- colorBin("YlOrRd", domain = nar_munic_stat2$PRODUCCION_1, bins = bins)

  mapa2 <- leaflet(data = nar_munic_stat2) %>%
  addTiles() %>%
  addPolygons(label = ~PRODUCCION_1,
              popup = ~MPIO_CNMBR,
              fillColor = ~pal(PRODUCCION_1),
              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_1,
    title = "Producción de Papa en Nariño [Ton/ha] (2018)",
    opacity = 1
  )
mapa2

A modo de conclusion, cabe recalcar que el departamento de Nariño en los ultimos años viene presentando un cambio en su producción agricola, en donde se observa una disminucion en los cultivos de cereales como el Trigo o de tuberculos como el cultivo de papa y por el contrario se presenta un aumento de la presencia de cultivos mas “tecnificados”. Esto se debe a que el sistema de produccion del departemento se basa principalmente en economia campesina minifundista. Sin embargo, en los ultimos años sectores mas agroindustriales estan tomando una meyor relevancion en el departamento como el sectro cafetero el cual el costo total de producción en el año se calcula en 200 mil millones de pesos de acuerdo a lo establecido por el ministerio de agricultura (2019).

6. Bibliografia

Ministerio de Agricultura. (2019). PLAN DEPARTAMENTAL DE EXTENSIÓN AGROPECUARIA DEL DEPARTAMENTO DE NARIÑO PDEA - NARIÑO (pp. 12-36). San Juan de Pasto, Colombia.

Agronet. (2018). Evaluaciones Agropecuarias Municipales, Nariño.