Influencia de la sequía en la producción agrícola en Sonora

La prducción agricola en Sonora

Introducción

Sonora es un estado en el cual la práctica agrícola es de mucha importancia ya que destaca por sus altos volúmenes de producción, los cuales se han logrado posicionar en el mercado internacional. las 10 hortalizas que más se producen en Sonora son: Sandía, papa, pepino, chile verde, calabacita, espárrago, calabaza, tomate, melón y cebolla.

Cabe aclarar que en la agricultura las condiciones climatológicas juegan un gran papel y el cambio clímatico puede afectar las condiciones en las que se producen los cultivos en la región, es por eso que en este proyecto analizaremos si la falta de agua por sequías ha afectado la producción de hortalizas en el estado de Sonora.

Objetivos

Hacer uso de datos de producción agrícola para saber si el cambio climático ha afectado la producción de hortalizas por la falta de agua que se ha visto en los últimos años mediante esta herramienta tecnológica la cual nos permitirá concluir si la producción se ha visto afectada durante estos últimos años y de esta manera sugerir distintas técnicas de riego u optar por otros tipos de siembre que requieran de menos cantidad de agua.

Antecedentes

Varios modelos de predicción de cambio climatico sugieren que en ciertas regiones ña productividad se podría reducir hasta en un %50, especialmente en zonas secas, de la cual, Sonora no está excento debido a que es una zona desértica.

Método

Se utilizará un método de análisis estadístico en el cual se analizará la producción agrícola, así como es que han variado las sequías a lo largo de los años.

Datos

library(readxl)
library(DT)
setwd("~/EAMJ1130")
datosAgricolasSonora2019 <- read_excel("SonoraAgricultureData.xlsx")
datosAgricolasSonora2010 <- read_excel("Sonora2010AgricultureData.xlsx")
datosAgricolasSonora2000 <- read_excel("Sonora2000AgricultureData.xlsx")
# Depuración de datos
datosRendimiento2019 <- (datosAgricolasSonora2019$Rendimiento)
datosRendimiento2010 <- (datosAgricolasSonora2010$Rendimiento)
datosRendimiento2000 <- (datosAgricolasSonora2000$Rendimiento)

aDatosRendimiento2019 <- cumsum(datosAgricolasSonora2019$Rendimiento)
aDatosRendimiento2010 <- cumsum(datosAgricolasSonora2010$Rendimiento)
aDatosRendimiento2000 <- cumsum(datosAgricolasSonora2000$Rendimiento)
datatable(datosAgricolasSonora2019)
library(ggplot2)
datosAgricolasSonora2019 <- read_excel("datos_AgricolasSonora2019.xlsx")
datosAgricolasSonora2010 <- read_excel("datos_AgricolasSonora2010.xlsx")
datosAgricolasSonora2000 <- read_excel("datos_AgricolasSonora2000.xlsx")


ggplot(datosAgricolasSonora2000, aes(fill=Nomcultivo, y=Sembrada, x= Nommodalidad)) + geom_bar(position="dodge", stat="identity") + xlab("Tipo de cultivo") + ylab("Hectareas sembradas") + ggtitle("Principales cultivos en Sonora en primavera-verano durante el 2000")

ggplot(datosAgricolasSonora2010, aes(fill=Nomcultivo, y=Sembrada, x= Nommodalidad)) + geom_bar(position="dodge", stat="identity") + xlab("Tipo de cultivo") + ylab("Hectareas sembradas") + ggtitle("Principales cultivos en Sonora durante el año 2010")

ggplot(datosAgricolasSonora2019, aes(fill=Nomcultivo, y=Sembrada, x= Nommodalidad)) + geom_bar(position="dodge", stat="identity") + xlab("Tipo de cultivo") + ylab("Hectareas sembradas") + ggtitle("Principales cultivos en Sonora durante el año 2019")

(Gráfica 3)

Sequia en Sonora

En Sonora, la falta de lluvias y el ambiente seco ha repercutido en la producción agrícola y ganadera, en este último ha reducido un 42% el hato ganadero; es decir, se han perdido 300 mil reses y se prevé que incremente a 456 mil entre mayo y junio. De acuerdo al Monitor de Sequía de México, dependiente del Servicio Meteorológico Nacional (SMN) y la Comisión Nacional del Agua (Conagua), en el último reporte del 30 de abril se contabilizaron 15 municipios con la condición más crítica de sequedad y en el informe del 31 de marzo había cinco en este nivel.

Además, el monitoreo más reciente señaló que 100% del territorio estatal se encuentra con sequía, la mayor parte (41%) enfrenta una situación severa, 38. 4% se encuentra en nivel extremo, mientras que 8.5% de la entidad en excepcional y el resto entre anormalmente seco y moderada.

A pesar de que toda la entidad es azotada por las inclemencias del clima y la escasez de lluvias hay municipios que atraviesan por una etapa más crítica, entre las que destacan Agua Prieta, Nogales, Naco, Sáric, Cananea, Altar, Tubutama, Magdalena, Ímuris, Bacoachi, Fronteras, Átil, Santa Cruz y Caborca, entre otros.

Los niveles que ha exhibido la Conagua están por alcanzar los obtenidos entre 2011 y 2012, período en qué se registró una de las peores sequías en México cuando alrededor del 95% de la nación se encontraba en esta situación, actualmente se tiene contabilizado hasta el 80% de territorio federal.

Visualizacion

Sequia anual en Sonora

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(readxl)
SeqAnual <- read_excel("~/EAMJ1130/SeqAnual.xlsx")
G <- ggplot(data = SeqAnual, aes(x = Fecha, y = D3, geom_area(stat = "bin", color = "blue", fill = "#00AFBB")))

  G

G1 <- plot_ly(data = SeqAnual, x = ~Fecha, y = ~D3, color = ~Fecha)
          
          G1
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

En estas dos graficas se muestra como ha ido avanzando la sequia extrema en Sonora desde 2002, a principios de la grafica se puede ver que hasta el 2008 se comporta de una manera considerable, pero en diciembre de 2011 la sequia sube de una manera anorma, en ese ano hubo miles de muertes ganadera en el estado, de igual manera afecto demasiado la agricultura. Si observamos desde 2019 ha empezado a subir la escases de agua en Sonora y siguen subiendo, hasta ahora son incontables el numero de perdidas de ganado y las perdidas agricolas que se han tenido desde 2019 y que seguira subiendo en el transcurso del ano.

Sequia extrema por municipio

library(readr)
SequiasSon <- read_csv("~/EAMJ1130/SequiasSon.csv")
## Warning: Missing column names filled in: 'X301' [301]
## 
## -- Column specification --------------------------------------------------------
## cols(
##   .default = col_character(),
##   `Clave del Municipio` = col_double(),
##   `30-abr-05` = col_logical(),
##   X301 = col_logical(),
##   `Anormalmente Seco (D0)` = col_double(),
##   Moderada_D1 = col_double(),
##   `Severa(D2)` = col_double(),
##   Extrema_D3 = col_double(),
##   `Excepcional(D4)` = col_double(),
##   `Total (D0 a D4)` = col_double()
## )
## i Use `spec()` for the full column specifications.
G2 <- plot_ly(data = SequiasSon, x = ~Municipio, y = ~Extrema_D3, color = ~Municipio)
          
          G2
## No trace type specified:
##   Based on info supplied, a 'bar' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#bar
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

Si bien tenemos en cuenta los cultivos que cosecha cada municipio por su localizacion vemos como la seqia extrema va aumentando en cada municipio afectando directamente los cultivos. Esto afecta en la produccion y en la calidad dde producto que se obtiene.

Analisis de datos

Correlacion lineal

 library(readr)
 SequiasSon <- read_csv("SequiasSon.csv")
## Warning: Missing column names filled in: 'X301' [301]
## 
## -- Column specification --------------------------------------------------------
## cols(
##   .default = col_character(),
##   `Clave del Municipio` = col_double(),
##   `30-abr-05` = col_logical(),
##   X301 = col_logical(),
##   `Anormalmente Seco (D0)` = col_double(),
##   Moderada_D1 = col_double(),
##   `Severa(D2)` = col_double(),
##   Extrema_D3 = col_double(),
##   `Excepcional(D4)` = col_double(),
##   `Total (D0 a D4)` = col_double()
## )
## i Use `spec()` for the full column specifications.
library(pacman)
p_load("MASS", "ggplot2")

Diagrama de dispersion para municipios con sequia extrema y moderada

ggplot(data = SequiasSon, aes(x= Municipio , y= Extrema_D3)) +
  geom_point(colour="red4") +
  ggtitle("Diagrama de dispersion de sequia extrema por municipio") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

Analisis de normalidad

par (mfrow = c(1,2))
hist(SequiasSon$'Moderada_D1', breaks = 10, main = "", xlab = "Sequia moderada", border="darkred")
hist(SequiasSon$'Extrema_D3', breaks = 10, main = "", xlab = "Sequia extrema", border="blue")

Grafico cuantilico

qqnorm(SequiasSon$'Moderada_D1', main = "Sequia moderada", col= "darkred")
qqline(SequiasSon$'Moderada_D1')

qqnorm(SequiasSon$'Extrema_D3', main = "Sequia extrema", col= "darkred")
qqline(SequiasSon$'Extrema_D3')

Prueba (test) de hipotesis para el analisis de normalidad

#Test de Shapiro Wilk para sequia moderada
shapiro.test(SequiasSon$`Moderada_D1`)
## 
##  Shapiro-Wilk normality test
## 
## data:  SequiasSon$Moderada_D1
## W = 0.93455, p-value = 0.001018

Como se puede observar el valor de P es 0.001018 lo que quiere decir que no alcanza el minimo (0.05 o 5%) para que sean considerados datos normales

#Test de Shapiro Wilk para Sequia extrema
shapiro.test(SequiasSon$`Extrema_D3`)
## 
##  Shapiro-Wilk normality test
## 
## data:  SequiasSon$Extrema_D3
## W = 0.85216, p-value = 5.88e-07

En este test el valor de p es sumamente pequeno y estos datos tampoco son normales

par (mfrow = c(1,2))
hist(SequiasSon$`Extrema_D3`, breaks = 10, main = "", xlab = " Log10 (Sequia extrema)", border="blue")
qqnorm(log10(SequiasSon$`Extrema_D3`), main = "", col = "blue")
qqline(log10(SequiasSon$`Extrema_D3`))

Ahora, con los datos ajustados a escala logaritmica de base, tenemos el siguiente analisis de normalidad

shapiro.test(log10(SequiasSon$`Extrema_D3`))
## 
##  Shapiro-Wilk normality test
## 
## data:  log10(SequiasSon$Extrema_D3)
## W = 0.95154, p-value = 0.007548

Ahora con los datos ajustados vemos que aun no alcanza el test de normalidad, lo que es probablemente un problema con lo que esta pasando en la vida diaria.

REGRESION LINEAL MULTIPLE

La regresión lineal múltiple permite generar un modelo lineal en el que el valor de la variable dependiente o respuesta (Y ) se determina a partir de un conjunto de variables independientes llamadas predictores (X1, X2, X3…). Es una extensión de la regresión lineal simple, por lo que es fundamental comprender esta última. Los modelos de regresión múltiple pueden emplearse para predecir el valor de la variable dependiente o para evaluar la influencia que tienen los predictores sobre ella (esto último se debe que analizar con cautela para no malinterpretar causa-efecto).

Los modelos lineales múltiples siguen la siguiente ecuación:

\[ Y_{i}=(\beta_{0}+\beta_{1}X_{1i}+\beta_{2}X_{2i}+\cdots+\beta_{n}X_{ni})+e_{i} \]

  • β0 es la ordenada en el origen, el valor de la variable dependiente Y cuando todos los predictores son cero.

  • βi : es el efecto promedio que tiene el incremento en una unidad de la variable predictora Xi sobre la variable dependiente Y , manteniéndose constantes el resto de variables. Se conocen como coeficientes parciales de regresión.

  • ei: es el residuo o error, la diferencia entre el valor observado y el estimado por el modelo.

Es importante tener en cuenta que la magnitud de cada coeficiente parcial de regresión depende de las unidades en las que se mida la variable predictora a la que corresponde, por lo que su magnitud no está asociada con la importancia de cada predictor. Para poder determinar qué impacto tienen en el modelo cada una de las variables, se emplean los coeficientes parciales estandarizados, que se obtienen al estandarizar (sustraer la media y dividir entre la desviación estándar) las variables predictoras previo ajuste del modelo.

Condiciones para la regresion lineal multiple

En si, mismas condiciones que los modelos lineales simples, mas algunas adicionales

Analisis de cosecha de trigo en el periodo primavera-verano 2020

#Librerias a utilizar
library(gplots)
## Warning: package 'gplots' was built under R version 4.0.5
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
library(dplyr)
library(psych)
## Warning: package 'psych' was built under R version 4.0.5
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(GGally)
## Warning: package 'GGally' was built under R version 4.0.5
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(readr)
library(readxl)
library(prettydoc) 
library(dplyr)
library(ggplot2)
library(data.table)
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
library(plotly)
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v tibble  3.0.6     v stringr 1.4.0
## v tidyr   1.1.2     v forcats 0.5.1
## v purrr   0.3.4
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x psych::%+%()          masks ggplot2::%+%()
## x psych::alpha()        masks ggplot2::alpha()
## x data.table::between() masks dplyr::between()
## x dplyr::filter()       masks plotly::filter(), stats::filter()
## x data.table::first()   masks dplyr::first()
## x dplyr::lag()          masks stats::lag()
## x data.table::last()    masks dplyr::last()
## x MASS::select()        masks dplyr::select(), plotly::select()
## x purrr::transpose()    masks data.table::transpose()
library(modelr)
DA2 <- read_excel("TrigoPrim-Ver2020.xlsx")

Tabla de predictores

En esta parte solo tomaremos en cuenta la valla del yaqui

Superficie_sembrada_Ha <- DA2$'SuperficieSembrada'
Superficie_cosechada_Ha <- DA2$'Superficiecosechada'
Superficie_siniestrada_Ha <- DA2$'Superficiesiniestrada'
Produccion <- DA2$'Produccion_ton'
Valor_de_la_produccion <- DA2$'Rendimiento(ton/ha)'

DA2 <- data.table(Superficie_sembrada_Ha, Superficie_cosechada_Ha, Superficie_siniestrada_Ha, Valor_de_la_produccion)

round( cor( x = DA2, method = "pearson"), 3)
##                           Superficie_sembrada_Ha Superficie_cosechada_Ha
## Superficie_sembrada_Ha                         1                      NA
## Superficie_cosechada_Ha                       NA                       1
## Superficie_siniestrada_Ha                     NA                      NA
## Valor_de_la_produccion                        NA                      NA
##                           Superficie_siniestrada_Ha Valor_de_la_produccion
## Superficie_sembrada_Ha                           NA                     NA
## Superficie_cosechada_Ha                          NA                     NA
## Superficie_siniestrada_Ha                         1                     NA
## Valor_de_la_produccion                           NA                      1

Distribucion de datos

library(psych)
multi.hist(x = DA2, dcol = c("blue", "red"), dlty = c("dotted", "solid"),main = "")

  • Ahora exploremos este mismo analisis usando GGally
library(GGally)
ggpairs(DA2, lower = list(continuous = "smooth"),
        diag = list(continuous = "barDiag"), axisLabels = "none")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (stat_bin).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 6 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (stat_bin).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 6 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value
## Warning: Removed 6 rows containing non-finite values (stat_smooth).
## Warning: Removed 6 rows containing missing values (geom_point).
## Warning: Removed 6 rows containing non-finite values (stat_smooth).
## Warning: Removed 6 rows containing missing values (geom_point).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 6 rows containing non-finite values (stat_bin).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 6 rows containing missing values
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 6 rows containing non-finite values (stat_smooth).
## Warning: Removed 6 rows containing missing values (geom_point).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (stat_bin).

Genereando el modelo

Ahora una vez que entendemos la forma en la cual se relacionan las vatiables, podemos empezar a experimentar con la generacion del modelo

modelo <- lm(Produccion ~ Superficie_sembrada_Ha + Superficie_cosechada_Ha + Superficie_siniestrada_Ha + Valor_de_la_produccion, data = DA2 )

summary(modelo)
## 
## Call:
## lm(formula = Produccion ~ Superficie_sembrada_Ha + Superficie_cosechada_Ha + 
##     Superficie_siniestrada_Ha + Valor_de_la_produccion, data = DA2)
## 
## Residuals:
## ALL 3 residuals are 0: no residual degrees of freedom!
## 
## Coefficients: (2 not defined because of singularities)
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)               -3.883e+01         NA      NA       NA
## Superficie_sembrada_Ha     8.167e+00         NA      NA       NA
## Superficie_cosechada_Ha           NA         NA      NA       NA
## Superficie_siniestrada_Ha         NA         NA      NA       NA
## Valor_de_la_produccion     5.145e-12         NA      NA       NA
## 
## Residual standard error: NaN on 0 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:      1,  Adjusted R-squared:    NaN 
## F-statistic:   NaN on 2 and 0 DF,  p-value: NA

Conclusión

La sequía es una factor climatológico que afecta a la producción de hortalizas en el estado de Sonora, ya que al realizar los análisis de producción se observo que la producción fue decreciendo durante los últimos años, mientras que la sequía fue incrementando, observandose en los anális de esta manera que esta es una tendencia que parece seguir continuando por los siguientes años, por lo que se sugeriría, que en la región se opte por utilizar cultivos que requieran de menos agua o cambiar técnicas de riego para que se use menos gasto de agua.

##Conclusión personal Con el análisis de datos realizado, observamos como el cambio de la producción desde el año 2000 al 2019 fue disminuyendo en gran cantidad, ya que en estos últimos años los datos de la sequía indican que se encuentra una tendencia de incremento, lo que significa que hay menos agua disponible para los cultivos, lo cual es causa del cambio climatico durante os últimos años, por lo que en este trabajo a realizar una regresión lineal en la producción del trigo se observa que existe una correlación etre la superficie cosechada y la siniestrada de 0.507, por lo que podemos concluir que la falta de agua es un factor que no permite que la producción agrícola crezca.

Referencias

Altieri, M. A., & Nicholls, C. (2008). Los impactos del cambio climático sobre las comunidades campesinas y de agricultores tradicionales y sus respuestas adaptativas. Agroecología, 3, 7-24. Recuperado a partir de https://revistas.um.es/agroecologia/article/view/95471

Sánchez, R. (2016, 5 septiembre). El monitor de la sequía en México. Scielo. http://www.scielo.org.mx/scielo.php?script=sci_arttext&pid=S2007-24222016000500197

Las 10 principales hortalizas producidas en sonora (2017). Recuperado del 10 de Diciembre del 2020, de http://oiapes.sagarhpa.sonora.gob.mx/notas/econo/10hortalizas.pdf