#Librerias con las que vamos a trabajar
library(tidyverse)
library(sf)
library(sp)
library(readxl)
library(grid)
library(ggplot2)
library(plotly)
library(leaflet)
library(osmdata)
library(httr)
library(hrbrthemes)
library(ggmap)
options(scipen = 999)
options(OutDec= ".")
options(digits = 2)
formato <- function(x){
  format(x, digits = 2, big.mark = ".", decimal.mark = ",")
}
#Bases con las que vamos a trabajar
estimaciones <- read_xlsx("TP_3.xlsx", sheet = "Hoja2")
retenciones <- read_xlsx("TP_3.xlsx", sheet = "Retenciones")
argentina_departamentos <- read_sf(dsn = "departamento.json")
argentina_provincias <- read_sf(dsn ="provincia.geojson")
#Modificaciones a las bases con las que vamos a trabajar
##Renombrar los shapes para luego poder unirlos de manera más sencilla
argentina_departamentos <- rename(argentina_departamentos, cde = in1)
argentina_provincias <- rename(argentina_provincias, cpr = IN1 )

##Renombrar las variables para facilitar el trabajo
estimaciones <- rename(estimaciones, Sup_Sembrada = "Sup. Sembrada (Ha)")
estimaciones <- rename(estimaciones, Sup_Cosechada = "Sup. Cosechada (Ha)")
estimaciones <- rename(estimaciones, Produccion = "Producción (Tn)")
estimaciones <- rename(estimaciones, Rendimiento = "Rendimiento (Kg/Ha)")

##Introducción

A partir del decreto 37/2019, el poder ejecutivo dejó sin efecto el límite de $4 por cada dólar para los derechos de exportación, establecido en el artículo 2º del Decreto 793/2018, y demás modificaciones. En la práctica significó el aumento de los derechos de exportación para los cultivos tradicionales de aproximadamente un 5%. En este trabajo nos concentraremos en la distribución territorial del cultivo de cebada, girasol, maíz, soja, sorgo y trigo por su relevancia en las cuentas nacionales de nuestro país.

El mencionado decreto puede ser encontrado a continuación: https://www.boletinoficial.gob.ar/detalleAviso/primera/223859/20191214

Contamos con información detallada por cultivo, campaña, provincia y departamento con los siguientes datos: superficie sembrada, superficie cosechada, producción y rendimiento. Las superficies estan expresadas en hectáreas, la producción en toneladas y el rendimiento en kg/ha. Estos datos son provenientes de la Dirección de Análisis Económico Agroindustrial - Dirección de Estimaciones agrícolas http://datosestimaciones.magyp.gob.ar/reportes.php?reporte=Estimaciones

##Una primera aproximación

Realizamos una primera aproximación a este dataset:

str(estimaciones)
## Classes 'tbl_df', 'tbl' and 'data.frame':    63936 obs. of  13 variables:
##  $ ID Provincia   : chr  "06" "06" "06" "06" ...
##  $ Provincia      : chr  "BUENOS AIRES" "BUENOS AIRES" "BUENOS AIRES" "BUENOS AIRES" ...
##  $ ID Departamento: chr  "854" "854" "854" "588" ...
##  $ Departamento   : chr  "25 DE MAYO" "25 DE MAYO" "25 DE MAYO" "9 DE JULIO" ...
##  $ Id Cultivo     : num  33 33 33 33 33 33 33 33 33 33 ...
##  $ Cultivo        : chr  "Cebada" "Cebada" "Cebada" "Cebada" ...
##  $ ID Campaña     : num  48 49 50 48 49 50 48 49 50 48 ...
##  $ Campana        : chr  "2016/17" "2017/18" "2018/19" "2016/17" ...
##  $ Sup_Sembrada   : num  4742 5660 13507 10087 3750 ...
##  $ Sup_Cosechada  : num  4742 4940 12697 8871 3750 ...
##  $ Produccion     : num  18020 15808 50788 39920 16875 ...
##  $ Rendimiento    : num  3800 3200 4000 4500 4500 5000 3600 3000 3200 3300 ...
##  $ cde            : chr  "06854" "06854" "06854" "06588" ...
summary(estimaciones)
##  ID Provincia        Provincia         ID Departamento   
##  Length:63936       Length:63936       Length:63936      
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
##                                                          
##  Departamento         Id Cultivo   Cultivo            ID Campaña
##  Length:63936       Min.   :12   Length:63936       Min.   : 1  
##  Class :character   1st Qu.:25   Class :character   1st Qu.:11  
##  Mode  :character   Median :26   Mode  :character   Median :23  
##                     Mean   :25                      Mean   :24  
##                     3rd Qu.:32                      3rd Qu.:38  
##                     Max.   :33                      Max.   :50  
##    Campana           Sup_Sembrada     Sup_Cosechada       Produccion     
##  Length:63936       Min.   :      1   Min.   :      0   Min.   :      0  
##  Class :character   1st Qu.:    700   1st Qu.:    450   1st Qu.:    750  
##  Mode  :character   Median :   3620   Median :   2700   Median :   5400  
##                     Mean   :  17593   Mean   :  15645   Mean   :  43846  
##                     3rd Qu.:  16000   3rd Qu.:  13000   3rd Qu.:  30985  
##                     Max.   :1096100   Max.   :1066100   Max.   :3399048  
##   Rendimiento        cde           
##  Min.   :    0   Length:63936      
##  1st Qu.: 1200   Class :character  
##  Median : 2000   Mode  :character  
##  Mean   : 2364                     
##  3rd Qu.: 3000                     
##  Max.   :33333

Donde vemos que tenemos variables numéricas, Sup. Sembrada (Ha) por ejemplo, y categóricas, Cultivo es una de ellas.

Primero filtraremos nuestra base solo para la última campaña para la que tenemos información (2018/19) y graficaremos la cantidad de hectáreas sembradas y la producción en toneladas por cultivo.

campana_2018_19 <- estimaciones %>% 
  filter(Campana == "2018/19") %>% 
  group_by(Cultivo) %>% 
  summarise(Produccion = sum(Produccion),
            Sup_Cosechada = sum(Sup_Cosechada),
            Sup_Sembrada = sum(Sup_Sembrada),
            Rendimiento = sum(Rendimiento))
#Calculamos diferentes datos que vamos a utilizar en nuestro informe
soja_sembrada <- campana_2018_19 %>% 
  filter(Cultivo == "Soja") %>%
  select(Sup_Sembrada)
soja_sembrada <- soja_sembrada$Sup_Sembrada
#falta coma de decimales
soja_sembrada_porcentaje <- campana_2018_19 %>%
  summarise(total = sum(Sup_Sembrada)) %>% 
  summarise(porcentaje = soja_sembrada/total*100)
soja_sembrada_porcentaje <- round(soja_sembrada_porcentaje$porcentaje, 2)
campana_2018_19 %>%
  arrange(desc(Sup_Sembrada)) %>%
  mutate(Cultivo = factor(Cultivo, levels = Cultivo)) %>% 
ggplot(aes(x = Cultivo, y = Sup_Sembrada, fill= Cultivo)) +
  geom_col()+
  scale_fill_hue(l=40) +
  theme_ipsum() +
      theme(axis.text.x = element_text(angle=45, vjust=0.6),
            axis.title.x = element_text(size = 12),
            axis.title.y = element_text(size = 12)) +
  
      scale_y_continuous(labels=function(x) format(x, big.mark = ".", decimal.mark = ",", scientific = FALSE)) +
      labs(title = "Superficie sembrada por cultivo tradicional",
       subtitle = "Campaña 2018/19",
       caption = "Fuente: elaboración propia en base a datos de Ministerio de Agricultura, Ganadería y Pesca",
        y = "Superficie sembrada (Ha)",
           x = "Cultivo") +
      guides(fill = FALSE)

Si nos concentramos en la superficie sembrada por cultivo para la campaña 2018/19, la soja ocupa el primer lugar con 17010277 hectáreas y representa el 46.61% del total sembrado.

#Calculamos diferentes datos que vamos a utilizar en nuestro informe
maiz_producido <- campana_2018_19 %>% 
  filter(Cultivo == "Maíz") %>%
  select(Produccion)
maiz_producido <- maiz_producido$Produccion
#falta coma de decimales
maiz_producido_porcentaje <- campana_2018_19 %>%
  summarise(total = sum(Produccion)) %>% 
  summarise(porcentaje = maiz_producido/total*100)
maiz_producido_porcentaje <- round(maiz_producido_porcentaje$porcentaje, 2)
campana_2018_19 %>%
  arrange(desc(Produccion)) %>%
  mutate(Cultivo = factor(Cultivo, levels = Cultivo)) %>% 
ggplot(aes(x = Cultivo, y = Produccion, fill= Cultivo)) +
  geom_col()+
  scale_fill_hue(l=40) +
  theme_ipsum() +
      theme(axis.text.x = element_text(angle=45, vjust=0.6),
            axis.title.x = element_text(size = 12),
            axis.title.y = element_text(size = 12)) +
      scale_y_continuous(labels=function(x) format(x, big.mark = ".", decimal.mark = ",", scientific = FALSE)) +
      labs(title = "Producción por cultivo tradicional",
       subtitle = "Campaña 2018/19",
       caption = "Fuente: elaboración propia en base a datos de Ministerio de Agricultura, Ganadería y Pesca",
        y = "Produccion (Tn)",
           x = "Cultivo") +
      guides(fill = FALSE)

Sin embargo, si nos enfocamos en la producción el maíz supera a la soja, con 56860704 toneladas en la campaña 2018/19 (40% del total). Esto se debe a los mayores rendimientos que se obtienen del maíz por hectarea, como se verifica en el gráfico nº3. Incluso el sorgo tiene un mayor rendimiento que la soja, pero como se cultivan una cantidad de hectáreas mucho menor, no se ubica entre los principales cuando se ve solo la cantidad de toneladas producidas.

campana_2018_19 %>%
  arrange(desc(Rendimiento)) %>%
  mutate(Cultivo = factor(Cultivo, levels = Cultivo)) %>% 
ggplot(aes(x = Cultivo, y = Rendimiento, fill= Cultivo)) +
  geom_col()+
  scale_fill_hue(l=40) +
  theme_ipsum() +
      theme(axis.text.x = element_text(angle=45, vjust=0.6),
            axis.title.x = element_text(size = 12),
            axis.title.y = element_text(size = 12)) +
      scale_y_continuous(labels=function(x) format(x, big.mark = ".", decimal.mark = ",", scientific = FALSE)) +
      labs(title = "Rendimiento por cultivo tradicional",
       subtitle = "Campaña 2018/19",
       caption = "Fuente: elaboración propia en base a datos de Ministerio de Agricultura, Ganadería y Pesca",
        y = "Rendimiento (Kg/Ha)",
           x = "Cultivo") +
      guides(fill = FALSE)

##Distribución territorial

Conocemos por estudios previos que el cultivo de estos productos agrícolas tradicionales se realiza en mayor medida en la región centro de nuestro país, pero nos gustaría con mayor detalle su distribución territorial. De esta manera, podremos preveer en que provincias y departamentos se verán afectados los productores.

Lo primero que tenemos que hacer es unir los datos de hectáreas sembradas y de producción con nuestro shapefile por departamentos.

#Filtramos el dataset estimaciones para quedarnos con la última campaña y agrupamos por departamento para calcular la superficie sembrada
campana_2018_19_sembrada <- estimaciones %>% 
  filter(Campana == "2018/19") %>% 
  group_by(cde) %>% 
  summarise(Sup_Sembrada = sum(Sup_Sembrada))

#Realizamos el mismo proceso para la producción por departamento
campana_2018_19_producida <- estimaciones %>% 
  filter(Campana == "2018/19") %>% 
  group_by(cde) %>% 
  summarise(Produccion = sum(Produccion))

#Creamos dos shapefiles nuevos: el primero con la superficie sembrada por departamento y el segundo con la producción por departamento
argentina_sembrada <- left_join(argentina_departamentos, campana_2018_19_sembrada)
## Joining, by = "cde"
argentina_producida <- left_join(argentina_departamentos, campana_2018_19_producida)
## Joining, by = "cde"
#Eliminamos el dataset que no vamos a utilizar más
remove(campana_2018_19_sembrada, campana_2018_19_producida)

Una vez que realizamos esto, podemos plotearlos. El primer mapa nos muestra la superficie sembrada por departamento. Las jurisdicciones que más se destacan se encuentran en Córdoba y Santa Fe.

ggplot() +
  geom_sf(data = argentina_sembrada, aes(fill = Sup_Sembrada), color = NA) +
  geom_sf(data = argentina_provincias, fill = NA, color = "gray60") +
  scale_fill_viridis_c(option = "magma") +
  theme_void() +
  theme(panel.grid.major = element_line(colour = "white")) +
  coord_sf(xlim = c(-52,-74), ylim = c(-20,-56), expand = FALSE) +
  labs(title = "Superficie sembrada de cultivos tradicionales por departamento",
       subtitle = "Campaña 2018/19",
       caption = "Fuente: elaboración propia en base a datos de Ministerio de Agricultura, Ganadería y Pesca",
       fill = "Hectáreas")

El segundo mapa nos muestra la producción de los mismos cultivos por departamento, con mínimos cambios. La región centro concentra gran parte de los cultivos tradicionales, así como el norte del país.

ggplot() +
  geom_sf(data = argentina_producida, aes(fill = Produccion), color = NA) +
  geom_sf(data = argentina_provincias, fill = NA, color = "gray60") +
  scale_fill_viridis_c(option = "magma")+
  theme_void() +
  theme(panel.grid.major = element_line(colour = "white")) +
  coord_sf(xlim = c(-52,-74), ylim = c(-20,-56), expand = FALSE) +
  labs(title = "Producción de cultivos tradicionales por departamento",
       subtitle = "Campaña 2018/19",
       caption = "Fuente: elaboración propia en base a datos de Ministerio de Agricultura, Ganadería y Pesca",
       fill = "Toneladas")

##Infraestructura vial relacionada

Un análisis que nos permiten los datos de Open Street Map es el de obtener información acerca de la infraestructura vial que conecta a los productores de estos cultivos con los puertos a través de los cuales se exporta.

Primero obtenemos un bounding box con los limites políticos de nuestro país.

#Debemos ser específicos con el nombre y el tipo de jurisdicción que queremos encontrar
#
bbox <- getbb("República Argentina", featuretype = "country")

Luego realizamos una consulta solicitando solo las rutas principales de la Argentina.

#rutas_principales <- opq(bbox, memsize = 891289600, timeout = 900)%>% #Utilizamos memsize porque necesitabamos mayor memoria para descargar los datos
 # add_osm_feature(key = "highway", value = "trunk") %>%
  #osmdata_sf()

Lo que también queremos obtener son los puertos.

#puertos_argentinos <- opq(bbox) %>%
 # add_osm_feature(key = "harbour") %>%
  #osmdata_sf()

Una vez que tenemos la infromación nos quedamos con lo que vamos a necesitar para plotear en el mapa.

#rutas <- rutas_principales$osm_lines
#puertos <- puertos_argentinos$osm_points

#rutas <- rutas %>% 
#  select(osm_id, name, FIXME, geometry)

#st_write(rutas, "rutas.shp")

#puertos <- puertos %>% 
 # select(osm_id, name, geometry)

#st_write(puertos, "puertos.shp")

rutas <- st_read(dsn = "rutas.shp")
## Reading layer `rutas' from data source `D:\RESTO\Curso Big Data\MODULO 3\TP\rutas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 21537 features and 3 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: -74 ymin: -55 xmax: -54 ymax: -21
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
puertos <- st_read(dsn = "puertos.shp")
## Reading layer `puertos' from data source `D:\RESTO\Curso Big Data\MODULO 3\TP\puertos.shp' using driver `ESRI Shapefile'
## Simple feature collection with 165 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -73 ymin: -55 xmax: -54 ymax: -22
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
ggplot() +
  geom_sf(data = rutas) +
 geom_sf(data = puertos, fill = "blue", color = "blue")

Si queremos limitarnos a tener las rutas y los puertos que se encuentran ubicados en nuestro país, podemos realizar un join espacial.

argentina_puertos <- st_join(puertos, argentina_provincias)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
argentina_rutas <- st_join(rutas, argentina_provincias)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

Luego filtramos aquellas rutas y puertos que no tengan nombres de provincias argentinas asociadas.

argentina_puertos <- argentina_puertos %>% 
    group_by(NAM) %>% 
    filter(NAM != is.na(NAM))

argentina_rutas <- argentina_rutas %>% 
  group_by(NAM) %>% 
  filter(NAM != is.na(NAM))

Para finalmente graficar las rutas y puertos en conjunto con la superficie sembrada de los cultivos tradicionales.

ggplot() +
  geom_sf(data = argentina_sembrada, aes(fill = Sup_Sembrada), color = NA) +
  geom_sf(data = argentina_provincias, fill = NA, color = "gray60") +
  geom_sf(data = argentina_puertos, size = 2.5, color = "blue") +
  geom_sf(data = argentina_rutas, color = "white") +
  scale_fill_viridis_c(option = "magma") +
  theme_void() +
  theme(panel.grid.major = element_line(colour = "white")) +
  coord_sf(xlim = c(-52,-74), ylim = c(-20,-56), expand = FALSE) +
  labs(title = "Superficie sembrada + infraestructura vial y portuaria",
       subtitle = "Campaña 2018/19",
       caption = "Fuente: elaboración propia en base a datos de Ministerio de Agricultura, Ganadería y Pesca y Open Street Map",
       fill = "Hectáreas")

Podemos visualizar la buena conectividad existente entre las áreas de mayor siembra y los puertos.

##Tweets

library(rtweet)
## 
## Attaching package: 'rtweet'
## The following object is masked from 'package:purrr':
## 
##     flatten

Tuvimos la intencion de utilizar el paquete rtweet para capturar tweets que hubiesen utilizado la palabra “retenciones”. Sin embargo, solo obtuvimos 2587 tweets, a pesar de pasar diferentes parámetros, de los cuales solo 6 tuvieron coordenadas.

#tweets_retenciones <- search_tweets(q = "retenciones",
 #               geocode = "-32.6284584,-62.7389599,500mi",
  #            include_rts = FALSE,
   #           n = 10000000,
    #          retryonratelimit = TRUE)
#coordenadas <- function(campo_coordenadas) {
 #   extraer_coordenadas <- function(lista_coords) {
  #      data_frame(lon = lista_coords[1],
   #                lat = lista_coords[2])
    #}
    #
  #  map_df(campo_coordenadas, extraer_coordenadas)
#}
#
#tweets_retenciones <- tweets_retenciones %>% 
 #   cbind(coordenadas(tweets_retenciones$coords_coords)) %>% 
  #  select(-geo_coords, -coords_coords, -bbox_coords) 

#tweets_retenciones_geo <- tweets_retenciones %>% 
 #   filter(!is.na(lat), !is.na(lon))