Informe tecnico Santander

Este cuaderno pretende mostrar un acercamiento hacia la economia agricola del departamento, basado en datos espaciales y no espaciales, especificamente un archivo de evaluaciones agropecuarias del departamento y datos del marco geoestadistico nacional del año 2017

##esta primera linea borra cualquier objeto que pudiera tener almacenado el programa

rm(list=ls())

estas librerias seràn de gran utilidad para el desarrollo de la pràctica, con esta linea nos aseguramos de que se instalen solo una vez en el programa

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)

here() starts at C:/Users/juane/Desktop

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.0     v purrr   0.3.3
v tibble  2.1.3     v dplyr   0.8.4
v tidyr   1.0.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-2, (SVN revision 621)
 GEOS runtime version: 3.6.1-CAPI-1.10.1 
 Linking to sp version: 1.4-1 
 Polygon checking: TRUE 

library(maptools)

Checking rgeos availability: TRUE

library(raster)

Attaching package: 㤼㸱raster㤼㸲

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

    select

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

    extract

library(sf)

Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

library(viridis)

Loading required package: viridisLite

library(rnaturalearth)
library(GSODR)

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

library(ggrepel)
library(cowplot)

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

En este objeto vamos a asignar los datos no espaciales asociados a las evaluaciones agropecuarias del departamento, donde se hace un resumen de los cultivos producidos en cada municipio, ademas de datos estadisticos como el area cultivada por cada elemento ademas de valores de produccion anual de los mismos

datos <- read_csv2("D:/unal/geomatica/santander/evaluacion agropecuaria_Santander.csv")

Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
Parsed with column specification:
cols(
  COD_DPTO = col_double(),
  DPTO = col_character(),
  COD_MUN = col_double(),
  MUN = col_character(),
  GRU_CUL = col_character(),
  SGRU_CUL = col_character(),
  CULT = col_character(),
  YEAR = col_double(),
  PERIODO = col_character(),
  AREA_SEM = col_double(),
  AREA_COS = col_double(),
  PROD = col_double(),
  REND = col_double(),
  ESTADO = col_character(),
  NOM_CIEN = col_character(),
  CICLO_CULTIVO = col_character()
)

datos

las siguientes dos lineas nos dejan ver unas muestras de los datos anteriores, la priemra nos deja ver las primeras seis filas y la segunda nos deja ver las ultimas seis del archivo

head(datos)

tail(datos)

vamos a iniciar el trabajo mostrando el rendimiento promedio de la produccion de los diferentes cultivos en cada uno de los muncipios, con esto podemos ver que grupo de cultivos predominan en tema de rendimiento anual en cada municipio

datos %>%
  group_by(MUN, GRU_CUL) %>%
  summarise(rend_prom = mean(REND, na.rm = TRUE)) -> REND_resumen
REND_resumen

##En esta nueva exploracion, vemos otra agrupacion de los datos, en donde se muestra el rendimiento promedio de cada grupo de cultivo en Ton/año en el departamento, donde observamos que se destacan las hortalizas y las oleaginosas, mas adelante veremos que importancia pueden tener estos datos

datos %>% group_by(GRU_CUL) %>% summarise(PROD_DPTO=mean(PROD,na.rm= TRUE))->PROD_Santander
PROD_Santander

aqui vamos a mostrar que municipio tuvo un mayor rendimiento para cada grupo de cultivo especificamente en el año 2018

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

REND_max_2018

para esta nueva visualizacion, realizaremos un filtro de los datos donde que municipio tuvo una mayor area cosechada para cada grupo de cultivo, observando que hay municipios que destinan una gran superficie para el tema agricola, como lo son Puerto Wilches

datos %>%filter(YEAR==2018) %>% group_by(GRU_CUL, MUN) %>%
  summarize(max_AREA_COS = max(AREA_COS, na.rm = TRUE)) %>%
    slice(which.max(max_AREA_COS)) -> area_cosecha_max

area_cosecha_max

aqui realizamos un filtro de los datos para explorar el grupo de cultivo Otros permanentes para ver que cual fue la maxima produccion de estos cultivos dentro del tiempo que cubren los datos

datos %>% filter(GRU_CUL=="Otros Permanentes") %>% group_by(MUN,SGRU_CUL="Cacao") %>% summarise (max_PROD=max(PROD, na.rm = TRUE)) %>% slice(which.max(max_PROD))->PROD_max
PROD_max

he decicido realizar una observacion de los datos en el municipio de San Benito, debido a que en la consulta anterior, se observo que este municipio presenta una alta produccion de cacao, llegando llegando hasta un maximo de 36000 hectareas del mismo

datos %>% 
  filter(MUN=="San Benito" & SGRU_CUL=="Cacao") %>% 
  group_by(YEAR, CULT) ->  VICEN_CACAO
VICEN_CACAO

para ver estos datos de una forma grafica, nos aprovecharemos de la funcion ggplot para vrlos de forma grafica y analizarlos de mejor manera, donde observamos que el año de mayor produccion en el municipio fue el 2013 llegando a las 200 toneladas en ese año

g <- ggplot(aes(x=YEAR, y=PROD/100), data = VICEN_CACAO) + geom_bar(stat='identity') + labs(y='Cacao Production [Ton x 1000]')

g + ggtitle("Evolution of cacao production in San Benito from 2013 to 2018") + labs(caption= "Based on EMA data (DANE, 2018)")

ahora, se realizara otra consulta, esta vez para observar la superficie cosechada total por cada grupo de cultivo en el año 2018, donde se observa que los cultivos Otros Permanentes detacan ampliamente sobre el resto

datos %>% 
  filter(YEAR==2018) %>% 
  group_by(GRU_CUL) %>%
  summarize(sum_AREA_COS = sum(AREA_COS, na.rm = TRUE)) %>%
     arrange(desc(sum_AREA_COS)) -> total_AREA_COS
total_AREA_COS

basados en lo anterior, vamos a realizar un filtro sobre el mismo año, pero tomando solo los cultivos otros permanentes, esto para ver cual cultivo en especifico es el que predomina respecto al area cosechada en el departamento. Como se observa, el cacao junto al cafe, presntan un gran porcentaje del total, donde destaca el cacao

datos %>% filter(GRU_CUL=="Otros Permanentes"& YEAR==2018) %>% group_by(CULT) %>% summarize(sum_cosecha = sum(AREA_COS, na.rm = TRUE)) %>%
     arrange(desc(sum_cosecha)) -> TOT_COSECHA

TOT_COSECHA

ahora que ya sabemos que el cacao es el cultivo otro permanente de mayor area cosechada, con este nuevo filtro veremos que municipio presenta una mayor area cosechada de cacao en el año 2018, donde vemos que el mayor exponente es el municipio San Vicente De Chucuri

datos %>% 
  filter(YEAR==2018 & GRU_CUL=="Otros Permanentes") %>% 
  group_by(CULT, MUN) %>%
  summarize(AREA_COS_MUN = max(AREA_COS, na.rm = TRUE)) %>%
    slice(which.max(AREA_COS_MUN)) ->max_area2 

max_area2

antes de realizar el ploteo, a cada una de las clases, que seran los grupos de cultivo, vamos a realizarles una abreviatura, para que se pueda leer con claridad en el ploteo

total_AREA_COS$CROP <- abbreviate(total_AREA_COS$GRU_CUL, 3)

este ploteo nos muestra el area cultivada por cada grupo de cultivo en el año 2018, veremos la misma informacion de hace varias lineas, pero con este grafico, sera mucho mas facil de sacar conclusiones sobre el mismo

g <- ggplot(aes(x=CROP, y=sum_AREA_COS), data = total_AREA_COS) + geom_bar(stat='identity') + labs(y='Total Harvested Area [Ha]')

g+ ggtitle("Total harvested area by crop groups in 2018 for Santander") + theme(plot.title = element_text(hjust = 0.5)) +
   labs(caption= "Based on EMA data (DANE, 2018)")

la segunda parte del informe consistira en realizar una union entre los datos no espacilaes mostrados en la primeraparte y un archivo de datos espaciales proveniente del marco geoestadistico del pais provenineto del DANE

Sant_MUN <-sf::st_read("D:/unal/geomatica/santander/marco geoestadisitico/68_SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Reading layer `MGN_MPIO_POLITICO' from data source `' using driver `ESRI Shapefile'
Simple feature collection with 87 features and 9 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -74.52895 ymin: 5.707536 xmax: -72.47706 ymax: 8.14501
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

veamos que tiene este archivo .shp

Simple feature collection with 87 features and 9 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -74.52895 ymin: 5.707536 xmax: -72.47706 ymax: 8.14501
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
   DPTO_CCDGO MPIO_CCDGO      MPIO_CNMBR
1          68      68001     BUCARAMANGA
2          68      68013          AGUADA
3          68      68020         ALBANIA
4          68      68051         ARATOCA
5          68      68077         BARBOSA
6          68      68079       BARICHARA
7          68      68081 BARRANCABERMEJA
8          68      68092         BETULIA
9          68      68101         BOLIVAR
10         68      68121         CABRERA
                             MPIO_CRSLC MPIO_NAREA
3                  Ordenanza 33 de 1919  166.21697      2017  SANTANDER
4                                  1750  169.79155      2017  SANTANDER
5  Ordenanza 30 del 25 de Abril de 1936   46.66489      2017  SANTANDER
6                                  1799  137.27581      2017  SANTANDER
7  Ordenanza 13 del 17 de Abril de 1922 1326.83512      2017  SANTANDER
8                                  1874  431.24871      2017  SANTANDER
9                                  1844 1010.11035      2017  SANTANDER
10                                 1869   65.57431      2017  SANTANDER
   Shape_Leng  Shape_Area                       geometry
1   0.6922752 0.012513526 POLYGON ((-73.08418 7.23063...
2   0.4758098 0.006146093 POLYGON ((-73.56261 6.24032...
3   0.8761299 0.013570466 POLYGON ((-73.73616 5.87092...
4   0.6746922 0.013882031 POLYGON ((-72.98158 6.76065...
5   0.2703415 0.003810882 POLYGON ((-73.58988 5.99809...
6   0.5610888 0.011223345 POLYGON ((-73.22126 6.73288...
7   2.7351901 0.108590745 POLYGON ((-73.6939 7.254447...
8   1.2718180 0.035286963 POLYGON ((-73.53993 7.15392...
9   4.3603864 0.082514928 POLYGON ((-74.50132 6.27574...
10  0.3515706 0.005360404 POLYGON ((-73.25696 6.6213,...

para realizar el join necesitamos un campo en comun entre ambos paquetes de datos, por lo que usaremos el codigo de municipio para ello, estas dos lineas nos serviran para crear dos campos nuevos a los datos no espaciales, conel objetivo de convertir la columna COD_MUN que es un objeto de tipo numerico a un objeto de tipo factor, cosa fundamental para realizar el join

datos_2 <- datos
datos_2$TEMP <-  as.character(datos_2$COD_MUN)

datos_2

ahora al objeto datos_2 le realizaremos un filtro que consistira en que solo tomara los datos en el que se analize el cacao en cada uno de los municipios del departamento. El resultado de esto se almacenara en un nuevo objeto llamado datos_3

datos_2 %>% filter(CULT == "Cacao")  -> datos_3
datos_3

Aqui vemos la clase del objeto recien creado, vemos que se trata de un dataframe, cosa que es buena para realizar el join

class(datos_3)

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

para este nuevo objeto tomaremos solo los campos que nos sean de real interes del objeto datos_3, entre ellos el campo creado previamente que es un factor

datos_4 <- datos_3 %>% dplyr::select(MUN, MPIO_CCDGO, YEAR, PROD, REND) 

ahora al resultado de la linea anterior vamos a organizarlo de manera cronologica, es decir, veremos los datos de cada municipio en un año, luego los del siguiente asi hasta cubrir todos los datos

datos_4 %>% 
  gather("YEAR", "PROD", "REND" , key = variable, value = number)

datos_4

ahora vamos a pasar de una organizacion de los datos en forma larga, es decir con muchas filas y pocas columnas a una orgnizacion larga con muchas columnas y pocas filas, de manera que en cada columna veamos los datos de produccion y rendimiento de cada añ en columnas independientes, el resultado de este proceso se alamacenara en un nuevo objeto de nombre datos_5

datos_4 %>% 
  group_by(MPIO_CCDGO) %>% 
  mutate(Visit = 1:n()) %>% 
  gather("YEAR", "PROD", "REND", key = variable, value = number) %>% 
  unite(combi, variable, Visit) %>%
  spread(combi, number) -> datos_5
head(datos_5)

ya que tenemos todos los elementos necesarios, procedemos a realizar el join entre los datos no espaciales y los datos espaciales, tomando como punto comun el campo MPIO_CCDGO

Sant_MUN_stat =left_join(Sant_MUN, datos_5, by="MPIO_CCDGO")
summary(Sant_MUN_stat)

 DPTO_CCDGO   MPIO_CCDGO           MPIO_CNMBR   MPIO_CRSLC
 68:87      68001  : 1   AGUADA         : 1   1899   :11  
            68013  : 1   ALBANIA        : 1   1887   : 9  
            68020  : 1   ARATOCA        : 1   1799   : 8  
            68051  : 1   BARBOSA        : 1   1743   : 2  
            68077  : 1   BARICHARA      : 1   1750   : 2  
            68079  : 1   BARRANCABERMEJA: 1   1762   : 2  
            (Other):81   (Other)        :81   (Other):53  
   MPIO_NAREA        MPIO_NANO        DPTO_CNMBR   Shape_Leng    
 Min.   :  19.69   Min.   :2017   SANTANDER:87   Min.   :0.2700  
 1st Qu.:  91.01   1st Qu.:2017                  1st Qu.:0.4765  
 Median : 201.73   Median :2017                  Median :0.7392  
 Mean   : 351.28   Mean   :2017                  Mean   :1.0298  
 3rd Qu.: 426.17   3rd Qu.:2017                  3rd Qu.:1.2240  
 Max.   :3174.28   Max.   :2017                  Max.   :4.3604                                                                   
   Shape_Area          MUN               PROD_10      
 Min.   :0.00161   Length:87          Min.   :   0.0   Min.   :   6.0  
 1st Qu.:0.00744   Class :character   1st Qu.:   4.0   1st Qu.:  61.5  
 Median :0.01649   Mode  :character   Median :  30.0   Median : 140.0  
 Mean   :0.02873                      Mean   : 468.1   Mean   : 661.3  
 3rd Qu.:0.03486                      3rd Qu.: 100.0   3rd Qu.: 375.0  
 Max.   :0.25946                      Max.   :7200.0   Max.   :6625.0  
                                      NA's   :40       NA's   :52      
    PROD_11          PROD_12           PROD_2           PROD_3       Min.   :  20.0   Min.   :  28.0   Min.   :   0.0   Min.   :   0.0  
 1st Qu.:  63.0   1st Qu.:  88.0   1st Qu.:   7.0   1st Qu.:  14.0  
 Median : 137.5   Median : 166.0   Median :  39.0   Median :  43.0  
 Mean   : 702.9   Mean   : 837.4   Mean   : 460.5   Mean   : 432.1   3rd Qu.: 423.2   3rd Qu.: 578.0   3rd Qu.: 100.0   3rd Qu.: 110.0  
 Max.   :6540.0   Max.   :6540.0   Max.   :7490.0   Max.   :5640.0  
 NA's   :55       NA's   :58       NA's   :40       NA's   :42      
     PROD_4           PROD_5           PROD_6           PROD_7       
 Min.   :   2.0   Min.   :   2.0   Min.   :   2.0   Min.   :   2.00  
 1st Qu.:  20.5   1st Qu.:  28.0   1st Qu.:  27.0   1st Qu.:  24.25  
 Median :  53.0   Median :  62.0   Median :  94.0   Median : 100.50  
 Mean   : 478.8   Mean   : 502.1   Mean   : 499.9   Mean   : 538.35  
 3rd Qu.: 135.5   3rd Qu.: 144.0   3rd Qu.: 237.0   3rd Qu.: 285.50  
 Max.   :5350.0   Max.   :5014.0   Max.   :4993.0   Max.   :5145.00  
 NA's   :43       NA's   :44       NA's   :46       NA's   :47       
     PROD_8           PROD_9           REND_1         REND_10    
 Min.   :   2.0   Min.   :   3.0   Min.   : 1.00   Min.   : 1.0  
 1st Qu.:  36.0   1st Qu.:  54.0   1st Qu.: 5.00   1st Qu.: 5.0  
 Median : 116.0   Median : 140.5   Median : 7.50   Median : 6.0  
 Mean   : 505.9   Mean   : 602.7   Mean   :25.12   Mean   :10.4  
 3rd Qu.: 270.0   3rd Qu.: 325.5   3rd Qu.:51.50   3rd Qu.: 7.0  
 Max.   :5000.0   Max.   :5400.0   Max.   :76.00   Max.   :55.0  
 NA's   :50       NA's   :51       NA's   :47      NA's   :52    
    REND_11          REND_12           REND_2          REND_3     
 Min.   : 1.000   Min.   : 1.000   Min.   : 1.00   Min.   : 1.00  
 1st Qu.: 5.000   1st Qu.: 5.000   1st Qu.: 5.00   1st Qu.: 5.00  
 Median : 6.000   Median : 6.000   Median : 8.00   Median : 7.00  
 Mean   : 9.812   Mean   : 8.862   Mean   :30.05   Mean   :20.82  
 3rd Qu.: 8.000   3rd Qu.: 8.000   3rd Qu.:55.25   3rd Qu.:47.50  
 Max.   :55.000   Max.   :55.000   Max.   :75.00   Max.   :83.00  
 NA's   :55       NA's   :58       NA's   :45      NA's   :43     
     REND_4          REND_5           REND_6          REND_7     
 Min.   : 1.00   Min.   :  1.00   Min.   : 1.00   Min.   : 4.00  
 1st Qu.: 5.00   1st Qu.:  5.00   1st Qu.: 5.00   1st Qu.: 5.00  
 Median : 6.00   Median :  6.00   Median : 5.00   Median : 6.00  
 Mean   :19.09   Mean   : 14.26   Mean   :12.53  
 3rd Qu.: 8.00   3rd Qu.:  8.00   3rd Qu.: 8.00   3rd Qu.: 7.25  
 Max.   :94.00   Max.   :116.00   Max.   :94.00   Max.   :75.00  
 NA's   :43      NA's   :44       NA's   :46      NA's   :47     
     REND_8          YEAR_1        YEAR_10    
 Min.   : 1.00   Min.   : 1.000   Min.   :2007   Min.   :2016   1st Qu.: 5.00   1st Qu.: 5.000   1st Qu.:2007   1st Qu.:2016  
 Median : 6.00   Median : 5.000   Median :2007   Median :2016  
 Mean   :12.22   Mean   :2009   Mean   :2016  
 3rd Qu.: 8.00   3rd Qu.: 6.000   3rd Qu.:2009   3rd Qu.:2016  
 Max.   :67.00   Max.   :55.000   Max.   :2017   Max.   :2018  
 NA's   :50      NA's   :51       NA's   :40     NA's   :52    
    YEAR_11        YEAR_12         YEAR_2         YEAR_3    
 Min.   :2017   Min.   :2018   Min.   :2008   Min.   :2009  
 1st Qu.:2017   1st Qu.:2018   1st Qu.:2008   1st Qu.:2009  
 Median :2017   Median :2018   Median :2008   Median :2009  
 Mean   :2017   Mean   :2018   Mean   :2010   Mean   :2010  
 3rd Qu.:2017   3rd Qu.:2010   3rd Qu.:2011  
 Max.   :2018   Max.   :2018   Max.   :2018   Max.   :2018  

     YEAR_4         YEAR_5         YEAR_6         YEAR_7    
 Min.   :2010   Min.   :2011   Min.   :2012   Min.   :2013  
 1st Qu.:2010   1st Qu.:2011   1st Qu.:2012   1st Qu.:2013   Median :2010   Median :2011   Median :2012   Median :2013  
 Mean   :2011   Mean   :2012   Mean   :2013   Mean   :2014  
 3rd Qu.:2011   3rd Qu.:2012   3rd Qu.:2013   3rd Qu.:2014  
 Max.   :2018   Max.   :2018   Max.   :2018   Max.   :2018  
 NA's   :43     NA's   :44     NA's   :46     NA's   :47    
     YEAR_8         YEAR_9              geometry 
 Min.   :2014   Min.   :2015   POLYGON      :87  
 1st Qu.:2014   1st Qu.:2015   epsg:4326    : 0  
 Median :2014   Median :2015   +proj=long...: 0  
 Mean   :2014   Mean   :2015                     
 3rd Qu.:2014   3rd Qu.:2015                     
 Max.   :2017   Max.   :2018                     
 NA's   :50     NA's   :51                       

para el mapa final estas librerias seran de gran ayuda

library(leaflet)

esta instruccion nos permite ver en un mapa del departamento los valores de produccion de cada departamente en un año especifico

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

  mapa <- leaflet(data = Sant_MUN_stat) %>%
  addTiles() %>%
  addPolygons(label = ~PROD_12,
              popup = ~MPIO_CNMBR,
              fillColor = ~pal(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 = ~PROD_12,
    title = "Cacao production in Santander [Ton] (2018)",
    opacity = 1
  )
  

este es el resultado

mapa