library(fields)
library(geoR)
library(akima)
library(dplyr)
library(sqldf)
library(tidyr)
library(sp)
library(sf)
library(maps)
library(ggplot2)

Planteamiento del Problema

La calidad del aire es un factor crítico en la salud pública y el bienestar ambiental.
En particular, el ozono troposférico (O\(_3\)) y el dióxido de nitrógeno (NO\(_2\)) son contaminantes atmosféricos regulados bajo los Estándares Nacionales de Calidad del Aire Ambiental (NAAQS, por sus siglas en inglés).

En el presente informe se analizan datos oficiales de la Agencia de Protección Ambiental de Estados Unidos (EPA) correspondientes al estado de Texas en el año 2024.
El problema que se plantea es determinar si los niveles de O\(_3\) y NO\(_2\) cumplen con los límites normativos establecidos, y en qué medida podrían representar un riesgo para la salud en un momento específico del tiempo.

Descripción de las Variables

Se utilizarán dos variables provenientes de la red de monitoreo de calidad del aire en Texas:

  • Ozono (O\(_3\)):
    Variable: Daily Max 8-hour Ozone Concentration, que corresponde al máximo promedio móvil de 8 horas de ozono reportado diariamente por cada estación de monitoreo.
    • Unidad: partes por millón (ppm).
    • Escala: Razón (el valor cero indica ausencia total de ozono; se dice que una estación excede los niveles de O\(_3\) cuando la concentración máxima registrada en el promedio móvil de 8 horas supera los 0.070 ppm).
  • Dióxido de Nitrógeno (NO\(_2\)):
    Variable: Daily Max 1-hour NO2 Concentration, que corresponde a la concentración máxima de NO\(_2\) registrada en una hora dentro de cada día.
    • Unidad: partes por billón (ppb).
    • Escala: Razón (el valor cero indica ausencia total de NO\(_2\); se dice que una estación excede los niveles de NO\(_2\) cuando la concentración máxima registrada en una hora del día supera los 100 ppb).

Escala Temporal

Ambas variables se registran diariamente durante el año 2024.
Para efectos de este análisis preliminar, se toma como referencia un momento específico; es decir, desde el 01 de enero de 2024 hasta el 31 de diciembre de 2024.

Mapa de Texas

load("OZ_coor.RData")

load("NO_coor.RData")

texas <- st_read("Texas_State/State_Boundary.shp")
## Reading layer `State_Boundary' from data source 
##   `G:\Mi unidad\Camilo UNAL\Estadística Espacial\Proyecto\Primera Entrega\Texas_State\State_Boundary.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -11871800 ymin: 2978934 xmax: -10409240 ymax: 4369694
## Projected CRS: WGS 84 / Pseudo-Mercator
oz_sf <- st_as_sf(OZ_coor, coords = c("Site.Longitude", "Site.Latitude"), crs = 4326)
no_sf <- st_as_sf(NO_coor, coords = c("Site.Longitude", "Site.Latitude"), crs = 4326)

Se presenta el mapa de Texas con las localizaciones de las estaciones desde las cuales se tomaron los datos.

ggplot() +
  geom_sf(data = texas, fill = "gray95", color = "black") +
  geom_sf(data = oz_sf, color = "blue", shape = 16, size = 3) +
  geom_sf(data = no_sf, color = "red", shape = 17, size = 3) +
  labs(title = "Estaciones de Monitoreo de O3 y NO2 en Texas (2024)",
       caption = "Fuente: EPA") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Modelamiento con el Ozono

Análisis de Estacionariedad

Primero, se importan los datos correspondientes al Ozono en el estado de Texas, Estados Unidos.

OZ <- read.csv("OZ.csv")[,-1]

OZg <- as.geodata(OZ)


head(OZ)
##    Easting Northing Ozone
## 1 -1663915  3669343 0.034
## 2 -1689317  3670665 0.037
## 3 -1822422  3509765 0.037
## 4 -1813926  3522080 0.039
## 5 -1797108  3475404 0.038
## 6 -1497689  3445674 0.030

par(mfrow = c(2, 1), mar = c(2.75, 3, 0.5, 0.5), mgp = c(1.7, 0.7, 0))

plot(x = OZ$Ozone, y = OZ$Northing,
     xlab = "Ozono", ylab = "Norte")

plot(x = OZ$Easting, y = OZ$Ozone,
     xlab = "Este", ylab = "Ozono")

De estos gráficos es posible ver que la variable presenta estacionariedad con respecto a la media y, además, no se evidencian problemas con respecto a heteroscedasticidad. Por lo tanto, no es necesario ajustar un modelo lineal para intentar alcanzar estacionariedad.

Semivariograma

A pesar de haber ajustado el modelo lineal, el Ozono en Texas parece ser estacionario desde un comienzo.

De esta forma, se halla el semivariograma empírico para estos dos casos.

Datos originales

vargrm1 <- variog(OZg, trend = "cte")
## variog: computing omnidirectional variogram
plot(vargrm1, pch = 19, col = "springgreen3",
     xlab = "Distancia",
     ylab = "Semivarianza",
     main = "Variograma empírico de Ozono (O3) en Texas 2024")

lines(vargrm1$u, vargrm1$v, col = "springgreen3")

Como se espera, a puntos cercanos hay menor variabilidad en cuanto a la semivarianza, como es de esperar. Se puede ver que se presenta un efecto pepita pequeño, pues la semivarianza no pasa por el orígen. Además, se puede ver que el rango ocurre en aproximadamente \(400.000\) unidades de distancia.

Resistente a datos atípicos

vg1_at <- variog(OZg, trend = "cte",
                 estimator.type = "modulus")
## variog: computing omnidirectional variogram
plot(vg1_at, pch = 19, col = "springgreen4",
     xlab = "Distancia",
     ylab = "Semivarianza",
     main = "Variograma empírico de Ozono (O3) en Texas 2024")

lines(vg1_at$u, vg1_at$v, col = "springgreen4")

Es posible ver que es similar al semivariograma basado en el promedio. Esto se puede estar dando debido a que no hay presencia de datos atípicos que estén afectando el nivel medio de la variable.

Modelamiento con el Dióxido de Carbono

Análisis de Estacionariedad

Se procede a importar los datos relacionados al dióxido de nitrógeno en el estado de Texas.

NO <- read.csv("NO.csv")[,-1]

NOg <- as.geodata(NO)


head(NO)
##    Easting Northing   NO
## 1 -1689317  3670665  5.9
## 2 -1822422  3509765 10.4
## 3 -1813926  3522080  1.0
## 4 -1797108  3475404  8.3
## 5 -1798947  3506464 19.9
## 6 -1497689  3445674 20.1

A forma de exploración, analizamos la relación que hay entre el dióxido de carbono y las coordenadas.

par(mfrow = c(2, 1), mar = c(2.75, 3, 0.5, 0.5), mgp = c(1.7, 0.7, 0))

plot(x = NO$NO, y = NO$Northing,
     xlab = "Dióxido de Nitrógeno", ylab = "Norte")
plot(x = NO$Easting, y = NO$NO,
     xlab = "Este", ylab = "Dióxido de Nitrógeno")

Parece ser estacionario con respecto a la media pues no tiene una tendencia que dependa del espacio.

Semivariograma

Al igual que con la variable de ozono, el dióxido de nitrógeno parece tener un comportamiento estacionario a través del espacio desde un principio.

Datos originales

vgOZ_1 <- variog(NOg, trend = "cte")
## variog: computing omnidirectional variogram
plot(vgOZ_1, pch = 19, col = "#CD4F39",
     xlab = "Distancia",
     ylab = "Semivarianza",
     main = "Variograma empírico de NO2 en Texas")

lines(vgOZ_1$u, vgOZ_1$v, col = "#CD4F39")

Resistente a datos atípicos

vgOZ_at <- variog(NOg, trend = "cte",
                  estimator.type = "modulus")
## variog: computing omnidirectional variogram
plot(vgOZ_at, pch = 19, col = "#912CEE",
     xlab = "Distancia",
     ylab = "Semivarianza",
     main = "Variograma empírico de NO2 en Texas")

lines(vgOZ_at$u, vgOZ_at$v, col = "#912CEE")