Datos de la estación de monitoreo en el Centro Cultural R. Mijares en Torreón Coahuila Los datos se obtuvieron desde el portal nacional de transparencia. Folio 051260800035525Exp.355-2025 Folio_051260800026325Exp.263-2025. Folio_051260800035525Exp.355-2025

Recuperación de datos

Conversión de datos alfanuméricos (caracter) a variables de tiempo La función myPOSIX cambia la hora 00:00 del día a las 24:00 del día anterior para efectos de visualización y cálculo Uso del paquete lubridate que manipula datos de fecha y tiempo para agruparlos de acuerdo a un propósito por una variable de tiempo

ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/data/"
infile <- "2calairetorreon.CSV"
nomarchi<-paste0(ruta,infile)
caliairehora<-read.csv(nomarchi,stringsAsFactors = FALSE)
caliairehora$tc<-as.POSIXct(caliairehora$Fecha,format = "%d/%m/%Y %H:%M",tz=Sys.timezone())

class(caliairehora$tc) <- c("myPOSIX", class(caliairehora$tc))

caliairehora$time2<- format.myPOSIX(caliairehora$tc)

caliairehora$date <- ymd_hms(caliairehora$time2)
#head(caliairehora$date)
caliairehora$pm10<-as.numeric(caliairehora$PM10)
caliairehora$pm25<-as.numeric(caliairehora$PM2.5)
caliairehora$o3<-as.numeric(caliairehora$O3)
caliairehora$no2<-as.numeric(caliairehora$NO2)
caliairehora$so2<-as.numeric(caliairehora$SO2)
caliairehora$co<-as.numeric(caliairehora$CO)
caliairehora$wd<-as.numeric(caliairehora$DV)
caliairehora$ws<-as.numeric(caliairehora$VV)

caliairehora$year<-year(caliairehora$date)
caliairehora$hour<-hour(caliairehora$date)


caliairehora %>% dplyr::select(date,year,hour) 
head(caliairehora)
aq<-caliairehora %>% dplyr::select(date,ws,wd,no2,o3,pm10,so2,co,pm25,year,hour) 
dim(aq)

[1] 15732 11

#str(aq)
aq %>% select(date, year, hour) %>% head()

OPENAIR Software libre Casos de uso libre

https://rpubs.com/mrahman/airdatavis

Quick summary of the data

summaryPlot(aq)

# Very easy approach for doing time average 
aq_daily <- timeAverage(aq, avg.time = "day")
head(aq_daily)
aq_monthly <- timeAverage(aq, avg.time = "month")
head(aq_monthly)
#summary(aq)

Windrose, frecuencia del conteo de la dirección del viento

windRose(aq,ws.int=.4,breaks=4)

windRose(aq,ws.int=.4,breaks=4,type = "year", layout = c(4, 2))

Windrose, Contaminantes

pollutionRose(aq, pollutant = "pm25",main="Frecuencia de conteo por dirección del viento")

pollutionRose(aq, pollutant = "pm10",cols = "jet")

pollutionRose(aq, pollutant = "no2",cols = c("yellow", "green"))

pollutionRose(aq, pollutant = "so2",cols = "hue")

pollutionRose(aq, pollutant = "o3",cols = c("red", "white", "blue"))

pollutionRose(aq, pollutant = "co")

pollutionRose(aq, pollutant = "pm25", type = "pm10", layout = c(4, 1))

pollutionRose(aq, pollutant = "pm25", statistic = "prop.mean")

pollutionRose(aq, pollutant = "pm25", normalise = TRUE, seg = 1)

pollutionRose(aq, type = "season",pollutant = "pm25")

pollutionRose(aq, type = "season",pollutant = "pm10")

pollutionRose(aq, type = "season",pollutant = "no2")

pollutionRose(aq, type = "season",pollutant = "so2")

pollutionRose(aq, type = "season",pollutant = "o3")

pollutionRose(aq, type = "season",pollutant = "co")

pollutionRose(aq, type = "monthyear",pollutant = "pm25")

pollutionRose(aq, type = "monthyear",pollutant = "pm10")

pollutionRose(aq, type = "monthyear",pollutant = "no2")

pollutionRose(aq, type = "monthyear",pollutant = "so2")

pollutionRose(aq, type = "monthyear",pollutant = "o3")

pollutionRose(aq, type = "monthyear",pollutant = "co")

pollutionRose(aq, type = "monthyear",pollutant = "pm25",grid.line=5)

pollutionRose(aq, type = "monthyear",pollutant = "pm10",grid.line=5)

pollutionRose(aq, type = "monthyear",pollutant = "no2",grid.line=5)

pollutionRose(aq, type = "monthyear",pollutant = "so2",grid.line=5)

pollutionRose(aq, type = "monthyear",pollutant = "o3",grid.line=5)

pollutionRose(aq, type = "monthyear",pollutant = "co",grid.line=5)

a<-pollutionRose(aq, type = "pm10",pollutant = "pm25",grid.line=5)

b<-pollutionRose(aq, type = "pm25",pollutant = "pm10",grid.line=5)

pollutionRose(aq, type = "pm25",pollutant = "no2",grid.line=5)

pollutionRose(aq, type = "pm25",pollutant = "so2",grid.line=5)

pollutionRose(aq, type = "pm25",pollutant = "o3",grid.line=5)

pollutionRose(aq, type = "pm25",pollutant = "co",grid.line=5)

print(a, split = c(1, 1, 2, 1))
print(b, split = c(2, 1, 2, 1), newpage = FALSE)

PercentilRose, Contaminantes

percentileRose(aq, pollutant = "pm25",
               percentile = c(25, 50, 75, 90, 95, 99, 99.9),
               col = "brewer1", key.position = "right", smooth = TRUE)

percentileRose(aq, pollutant = "pm10",
               percentile = c(25, 50, 75, 90, 95, 99, 99.9),
               col = "brewer1", key.position = "right", smooth = TRUE)

percentileRose(aq, pollutant = "no2",
               percentile = c(25, 50, 75, 90, 95, 99, 99.9),
               col = "brewer1", key.position = "right", smooth = TRUE)

percentileRose(aq, pollutant = "so2",
               percentile = c(25, 50, 75, 90, 95, 99, 99.9),
               col = "brewer1", key.position = "right", smooth = TRUE)

percentileRose(aq, pollutant = "o3",
               percentile = c(25, 50, 75, 90, 95, 99, 99.9),
               col = "brewer1", key.position = "right", smooth = TRUE)

percentileRose(aq, poll="pm25", percentile = 95, method = "cpf",
               col = "darkorange", smooth = TRUE)

percentileRose(aq, poll="pm25", type = c("season", "daylight"),
col = "Set3", mean.col = "black")

timeProp(aq, pollutant = "pm25", proportion = "wd", avg.time = "week")

PolarPlot

polarPlot(aq, pollutant = "pm25", col = "jet", key.position = "bottom",
          key.header = "mean pm25 (ug/m3)", key.footer = NULL)

polarFreq(aq)

polarFreq(aq, type = "year")

polarAnnulus(aq, poll = "pm10", period = "trend", main = "Trend")

polarAnnulus(aq, poll = "pm10",period = "hour", main = "Hour")

polarAnnulus(aq, poll = "pm25", period = "trend", main = "Trend")

TimePlot

timePlot(selectByDate(aq, year = 2025, month = "abr"),
pollutant = c("so2","no2", "o3", "pm25", "pm10", "ws"),
y.relation = "free")

timePlot(aq, pollutant = "pm25", ylab = "pm25 (ug/m3)")

timePlot(selectByDate(aq, year = 2024, month = "abr"), pollutant = "pm25", avg.time = "day")

TimeVariation

timeVariation(subset(aq),
              pollutant = "pm25", ylab = "pm25 (ug/m3)")

timeVariation(aq, pollutant = c("pm25", "co", "no2", "o3"), normalise = TRUE)

timeVariation(aq, pollutant = "pm10", statistic = "median",
col = "firebrick")

timeVariation(aq, pollutant = "pm25", statistic = "median",
col = "firebrick")

myOutput <- timeVariation(aq, pollutant = "pm25")

head(myOutput$data$day.hour)
plot(myOutput, subset = "hour")

trendLevel(aq, x = "season", y = "daylight", pollutant = "pm25")

trendLevel(aq, pollutant = "pm25", y = "wd", border = "white",
cols = "jet")

Regresión, modelo lineal

lm_aq <- aq %>% filter(year == "2024") 
lm_aq <-  lm (data = aq, formula = pm25 ~ pm10 + co + hour)
summary(lm_aq)

Call: lm(formula = pm25 ~ pm10 + co + hour, data = aq)

Residuals: Min 1Q Median 3Q Max -195.784 -3.746 -0.529 2.960 296.674

Coefficients: Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.71380 0.18715 19.84 <2e-16 pm10 0.17400 0.00115 151.33 <2e-16 co 5.90596 0.12391 47.66 <2e-16 hour -0.28301 0.01294 -21.87 <2e-16 — Signif. codes: 0 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘.’ 0.1 ’ ’ 1

Residual standard error: 10.69 on 15138 degrees of freedom (590 observations deleted due to missingness) Multiple R-squared: 0.7138, Adjusted R-squared: 0.7137 F-statistic: 1.258e+04 on 3 and 15138 DF, p-value: < 2.2e-16

Calendario

calendarPlot(aq, pollutant = "pm25", cols = "viridis")

calendarPlot(aq, pollutant = "pm25", year =2025)

calendarPlot(aq, pollutant = "pm25", year = 2024,
breaks = c(0, 15, 41, 79, 130,400),
labels = c("buena", "acepta", "mala", "muy mala","extmala"),
cols = c("green", "yellow", "orange", "red","purple"),
statistic = "max")

calendarPlot(aq, pollutant = "pm25", 
breaks = c(0, 15, 41, 79, 130,400),
labels = c("buena", "acepta", "mala", "muy mala","extmala"),
cols = c("green", "yellow", "orange", "red","purple"),
statistic = "max")

calendarPlot(aq, pollutant = "pm25", year = 2024, annotate = "value",
lim =41, cols = "Purples", col.lim = c("black", "orange"),
layout = c(4, 3))

#### FUNCIÓN Para calcular promedios móviles #Se crea un nuevo campo de rollingPM25

aq <- rollingMean(aq, pollutant = "pm25", hours = 8, new.name = "rollingPM25", data.thresh = 75)

tail(aq)
calendarPlot(aq, pollutant = "rollingPM25", year = 2025,
breaks = c(0, 15, 41, 79, 130,400),
labels = c("buena", "acepta", "mala", "muy mala","extmala"),
cols = c("green", "yellow", "orange", "red","purple"),
statistic = "max")

## Cálculo de dos promedios pero que sólo se calcula si hay al menos el 75% de datos disponibles


twoweek <- timeAverage(aq, avg.time = "2 week", data.thresh = 75)
head(twoweek)
### correlación de datos

corPlot(subset(select(aq,-c(year,hour))), dendrogram = TRUE)

aq <- rollingMean(aq, pollutant = "pm25", hours = 8, new.name = "rollingPM25", data.thresh = 75)

tail(aq)

Relación entre variables

scatterPlot(aq, x = "pm25", y = "pm10")

scatterPlot(aq, x = "pm25", y = "pm10",smooth=TRUE)

scatterPlot(aq, x = "pm25", y = "pm10", method = "hexbin", cols = "inferno")

Estimaciones de tendencia Theilsend

TheilSen(aq, pollutant = "pm25")

Estimaciones tendencias suaves no paramétricas

smoothTrend(aq, pollutant = "pm25")

smoothTrend(aq, pollutant = "pm25", 
            type = "wd",
            date.breaks = 4)

smoothTrend(aq, pollutant = c("pm25", "pm10"), type = c("wd", "year"),
date.breaks = 4, lty = 0)