Ejercicio 2 de modelación ambiental

Este es un ejercicio enfocado en modelar la concentración de O3 en la tropósfera utilizando datos de la estación de monitoreo de la ERNO UNAM de Hermosillo y datos de movilidad urbana de “Google Mobility reports”

Librerías

Para instalar librerías usamos por ejemplo: install.packages(“ggplot2”) y ejecutando en la consola

library(ggplot2)
library(scales)
library(prettydoc)
library(randomForest)

## randomForest 4.7-1.1

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:ggplot2':
## 
##     margin

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

Importando los datos

En este apartado importaremos los datos de calidad del aire de la estación ERNO

#Leer y dar formato basico a la primer tabla 
datos <- read.csv("2020-01-ERNO_hora_L1_caire.csv", skip=7)
datos <- datos[-1,]
#Seleccionar nuestras columnas de interes 
datos <- datos[, c("Time", "O3", "PM10")]
#cambiamos el nombre de Time a tiempo
colnames(datos)[1] <- "tiempo"
#convertimos la columna time de tipo chr a formato tiempo
datos$tiempo <- as.POSIXct(datos$tiempo, format = "%Y-%m-%d %H:%M:%S"  )
#Convertimos las columnas O3 y PM10 
datos$O3 <- as.numeric(datos$O3)
datos$PM10 <- as.numeric(datos$PM10)

Visualizando la primer serie de tiempo de contaminantes atmosféricos

Primer gráfico (datos horarios )

Para esto utilizaremos la librería ggplot2 de tidyverse

ggplot(data= datos, aes(x = tiempo)) +
  geom_line(aes(y = O3, color ="O3")) +
  geom_line(aes(y = PM10, color ="PM10")) +
  scale_color_manual(values = c("blue","red"  )) +
  labs(x ="Tiempo", y = "Valor", color = "Contaminante") +
  scale_x_datetime(date_labels= "%d %m %H:%M", date_breaks = "1 hour"  ) +
  theme_bw()

### Segundo gráfico (datos diarios y semanales) Para esto utilizaremos la librería ggplot2 de tidyverse

ggplot(data= datos, aes(x = tiempo)) +
  geom_line(aes(y = O3, color ="O3")) +
  geom_line(aes(y = PM10, color ="PM10")) +
  scale_color_manual(values = c("blue","red" ), name = "Contaminante" ) +
  labs(x ="Tiempo", y = "Valor", title ="Valores de contaminantes atmosfericos de Hermosillo") +
  ggtitle("Valores de contaminantes atmosfericos de Hermosillo") +
  scale_x_datetime(date_labels= "%d %m %H:%M", date_breaks = "1 week"  ) +
  theme_bw() +
  theme(plot.title = element_text(hjust= 0.5  ),
        legend.position = "bottom",
        legend.box = "horizontal")

## Uniendo los datos de enero y febrero

#Leer y dar formato basico a la primer tabla 
datos2 <- read.csv("2020-02-ERNO_hora_L1_caire.csv", skip=7)
datos2 <- datos2[-1,]
#Seleccionar nuestras columnas de interes 
datos2 <- datos2[, c("Time", "O3", "PM10")]
#cambiamos el nombre de Time a tiempo
colnames(datos2)[1] <- "tiempo"
#convertimos la columna time de tipo chr a formato tiempo
datos2$tiempo <- as.POSIXct(datos2$tiempo, format = "%Y-%m-%d %H:%M:%S"  )
#Convertimos las columnas O3 y PM10 
datos2$O3 <- as.numeric(datos2$O3)
datos2$PM10 <- as.numeric(datos2$PM10)

unir (bind)

enefeb <- rbind(datos,datos2)

incorporando el mes de marzo

#Leer y dar formato basico a la primer tabla 
datos3 <- read.csv("2020-03-ERNO_hora_L1_caire.csv", skip=7)
datos3 <- datos3[-1,]
#Seleccionar nuestras columnas de interes 
datos3 <- datos3[, c("Time", "O3", "PM10")]
#cambiamos el nombre de Time a tiempo
colnames(datos3)[1] <- "tiempo"
#convertimos la columna time de tipo chr a formato tiempo
datos3$tiempo <- as.POSIXct(datos3$tiempo, format = "%Y-%m-%d %H:%M:%S"  )
#Convertimos las columnas O3 y PM10 
datos3$O3 <- as.numeric(datos3$O3)
datos3$PM10 <- as.numeric(datos3$PM10)

unir (bind)

enefebmar<- rbind(datos,datos2, datos3)

Verificando si existen datos nulos

colSums(is.na(enefebmar))

## tiempo     O3   PM10 
##      0      0      0

Graficando el ozono para los 3 primeros meses de 2020

plot(enefebmar$tiempo, enefebmar$O3)

Promedios diarios

Para poder unir las tabla de los contaminantes con la de movilidad urbana y de esta forma poder crear modelos, necesitamos hacer promedios diarios de los contaminantes dado que los datos de movilidad son datos diarios

promedios_diarios <- aggregate(enefebmar[, 2:3 ], by =
list(as.Date(enefebmar$tiempo)), mean )
colnames(promedios_diarios)[1] <- "tiempo"
promedios_diarios$tiempo <- as.POSIXct(promedios_diarios$tiempo)
str(promedios_diarios)

## 'data.frame':    92 obs. of  3 variables:
##  $ tiempo: POSIXct, format: "2020-01-01" "2020-01-02" ...
##  $ O3    : num  25.4 20.2 20.9 16.7 17.6 ...
##  $ PM10  : num  16.8 23.9 25.3 60.9 48.4 ...

Grafico de los primeros 3 meses de 2020

ggplot(data= promedios_diarios, aes(x = tiempo)) +
  geom_line(aes(y = O3, color ="O3")) +
  geom_line(aes(y = PM10, color ="PM10")) +
  scale_color_manual(values = c("blue","red" ), name = "Contaminante" ) +
  labs(x ="Tiempo", y = "Valor", title ="Valores de contaminantes atmosfericos de Hermosillo") +
  ggtitle("Valores de contaminantes atmosfericos de Hermosillo") +
  scale_x_datetime(date_labels= "%d %m %H:%M", date_breaks = "2 week"  ) +
  theme_bw() +
  theme(plot.title = element_text(hjust= 0.5  ),
        legend.position = "bottom",
        legend.box = "horizontal")

## Primer modelo de regresion lineal

la premisa es que mayor concentracion de PM10 indica que hay mas vehiculos circulando, entonces por efecto de smog fotoquimico se crea mas ozono en la troposferca

modelo <- lm(O3 ~ PM10, data = enefebmar)
summary(modelo)

## 
## Call:
## lm(formula = O3 ~ PM10, data = enefebmar)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.773 -11.598  -1.583  11.181  49.089 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.405915   0.411642   61.72   <2e-16 ***
## PM10        -0.105892   0.007719  -13.72   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.88 on 2182 degrees of freedom
## Multiple R-squared:  0.07939,    Adjusted R-squared:  0.07897 
## F-statistic: 188.2 on 1 and 2182 DF,  p-value: < 2.2e-16

Grafico de la recta de minimos cuadrados

ggplot(enefebmar, aes(x=PM10, y=O3)) +
  geom_point() +
  geom_smooth(method="lm", color = "red") +
  labs ( title = "Regresion lineal PM10 vs O3")

## `geom_smooth()` using formula = 'y ~ x'

Datos de movilidad urbana de google

movilidad <- read.csv("2020_MX_Region_Mobility_Report.csv")
movilidad <- subset(movilidad, iso_3166_2_code =="MX-SON")
movilidad <- subset(movilidad, select = -(1:8))
colnames(movilidad)[1] <- "tiempo"
colnames(movilidad)[2] <- "recreacion"
colnames(movilidad)[3] <- "abarrotesfarmarcias"
colnames(movilidad)[4] <- "parques"
colnames(movilidad)[5] <- "estaciones"
colnames(movilidad)[6] <- "trabajo"
colnames(movilidad)[7] <- "residencial"

formateando el tiempo

movilidad$tiempo <- as.POSIXct(movilidad$tiempo, format = "%Y-%m-%d")
str(movilidad)

## 'data.frame':    321 obs. of  7 variables:
##  $ tiempo             : POSIXct, format: "2020-02-15" "2020-02-16" ...
##  $ recreacion         : int  7 6 4 3 1 1 0 1 3 -3 ...
##  $ abarrotesfarmarcias: int  6 6 4 -1 -1 0 -1 -1 0 -1 ...
##  $ parques            : int  0 -1 2 5 3 12 3 5 6 8 ...
##  $ estaciones         : int  -7 -3 3 1 1 3 8 -1 0 0 ...
##  $ trabajo            : int  3 2 9 8 5 4 10 3 2 0 ...
##  $ residencial        : int  -1 0 -1 -1 -1 -1 -1 0 0 1 ...

Ajustando ventanas de tiempo

Datos de movilidad

movilidad_filtrado <- subset(movilidad, tiempo >= as.POSIXct("2020-02-15") & tiempo <= as.POSIXct("2020-03-31")       )

Datos de contaminantes

contaminantes_filtrado <- subset(promedios_diarios, tiempo >= as.POSIXct("2020-02-14") & tiempo <= as.POSIXct("2020-03-31")       )

Uniendo tablas de movilidad y de contaminantes

movcont <- merge(contaminantes_filtrado,movilidad_filtrado, by = "tiempo" )
str(movcont)

## 'data.frame':    46 obs. of  9 variables:
##  $ tiempo             : POSIXct, format: "2020-02-15" "2020-02-16" ...
##  $ O3                 : num  23.7 23.3 17.2 13.9 15.5 ...
##  $ PM10               : num  47.6 32.6 44 73.9 62.7 ...
##  $ recreacion         : int  7 6 4 3 1 1 0 1 3 -3 ...
##  $ abarrotesfarmarcias: int  6 6 4 -1 -1 0 -1 -1 0 -1 ...
##  $ parques            : int  0 -1 2 5 3 12 3 5 6 8 ...
##  $ estaciones         : int  -7 -3 3 1 1 3 8 -1 0 0 ...
##  $ trabajo            : int  3 2 9 8 5 4 10 3 2 0 ...
##  $ residencial        : int  -1 0 -1 -1 -1 -1 -1 0 0 1 ...

Modelando la concentracion de ozono troposferico usando el algoritmo de bosques aleatorios

movcont <- movcont[, -which(names(movcont) == "tiempo")]
model <- randomForest(O3 ~ ., data = movcont)

Evaluacion del modelo

set.seed(123)
train_index <- sample(nrow(movcont), 0.8 * nrow(movcont))
train_data <- movcont[train_index, ]
test_data <- movcont[-train_index, ]

#entrenando el modelo para su testeo 
model <- randomForest(O3 ~ ., data = movcont)
prediction <- predict(model, newdata = test_data)
cor(prediction, test_data$O3)^2

## [1] 0.8735303

Importancia relativa de las variables

importance <- importance(model)

plot_ly(x = rownames(importance), y = importance[, 1]) %>%
  add_bars()%>%
  layout(title = "Importancia relativa de las variables")

modelacion2

Badouin

2023-04-25

Ejercicio 2 de modelación ambiental

Librerías

Importando los datos

Visualizando la primer serie de tiempo de contaminantes atmosféricos

Primer gráfico (datos horarios )

unir (bind)

incorporando el mes de marzo

unir (bind)

Verificando si existen datos nulos

Graficando el ozono para los 3 primeros meses de 2020

Promedios diarios

Grafico de los primeros 3 meses de 2020

Grafico de la recta de minimos cuadrados

Datos de movilidad urbana de google

formateando el tiempo

Ajustando ventanas de tiempo

Uniendo tablas de movilidad y de contaminantes

Modelando la concentracion de ozono troposferico usando el algoritmo de bosques aleatorios

Evaluacion del modelo

Importancia relativa de las variables