Workflowcaliaire
Claudia Guerrero S
2025-08-07
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 051260800019725 Exp.197-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$t<-as.POSIXct(caliairehora$Fecha,format = "%d/%m/%Y")
caliairehora$tc<-as.POSIXct(caliairehora$Fecha,format = "%d/%m/%Y %H:%M",tz=Sys.timezone())
caliairehora$tcm<-as.POSIXct(caliairehora$Fecha,format = "%d/%m/%Y %H:%M",tz=Sys.timezone())
class(caliairehora$tcm) <- c("myPOSIX", class(caliairehora$tcm))
caliairehora$time2<- format.myPOSIX(caliairehora$tcm)
caliairehora<-caliairehora %>%
mutate(Date = lubridate::ymd(caliairehora$t),
anual= format(Date, "%Y"),
date_dia = format(Date, "%Y-%m-%d"),
tiempo =lubridate::ymd_hms(caliairehora$tc,quiet=FALSE),
time=format(tiempo,"%H:%M"),
horadia=substr(time2, start=12,stop=16))
# Conversión a datos numéricos
caliairehora$PM10<-as.numeric(caliairehora$PM10)
caliairehora$PM2.5<-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)Limpieza de datos
Conteo de datos NA y su frecuencia Cálculo de promedios en datos NA Elección de atributos y su orden primera graficación outliers y relación de partículas PM2.5 y PM10
calidadairehor <-caliairehora %>%
dplyr::select(contains("tc"),contains("tcm"),contains("Fecha"),contains("date_dia"),contains("horadia"),contains("PM10"),contains("PM2.5"),contains("O3"),contains("NO2"),contains("SO2"),contains("CO"))
nu_Total_NA <-sum(is.na(calidadairehor)) # count the number of NA
nu_Total_non_NA <-sum(!is.na(calidadairehor)) # number of non- NA values
nu_Total <- (nu_Total_NA + nu_Total_non_NA)
(fract_NA <-nu_Total_NA/nu_Total) [1] 0.06125326
#** SUMMARY** About 5.04% of total data from merged file is NA
dat1<-calidadairehor
is.na(dat1) <- do.call(cbind,lapply(dat1, is.infinite))
#Conteo de infinitos
sum(do.call(cbind,lapply(dat1, is.infinite)))[1] 0
#Conteo de NAs
sum(do.call(cbind,lapply(dat1, is.na)))[1] 10600
#conteo de no infinitos y no faltantes
sum(do.call(cbind,lapply(dat1, is.nan)))[1] 0
dat2 <- dat1 %>%
mutate_all(funs(replace(., .<0, NA)))
(nu_NA_in_dat1 <-sum(do.call(cbind,lapply(dat1, is.na)))) # number of NA in dat1[1] 10600
(nu_NA_in_dat2 <-sum(do.call(cbind,lapply(dat2, is.na)))) # number of NA in dat2[1] 12397
(nu_Neg_Values_in_dat1 <- (nu_NA_in_dat2-nu_NA_in_dat1))[1] 1797
#Reemplazo por la media en valores NA
dat3 <- data.frame(lapply(dat2,function(x) {
if(is.numeric(x)) ifelse(is.na(x),median(x,na.rm=T),x) else x}))
# Count of na in dat3
(sum(do.call(cbind,lapply(dat3, is.na))))[1] 0
#Nuevas variables
# library(lubridate)
dat3$Year <-year(dat3$tcm) # create new variable year
dat3$month <-month(dat3$tcm) # create new variable month
dat3$day <-day(dat3$tcm) # create new variable day
dat3$wday <-wday(dat3$tcm, label = TRUE) # create new variable week day (Factor Type)
dat3$hour <-hour(dat3$tcm) # create new variable hour
dat3$minute <-minute(dat3$tcm) # create new variable minute
#Elección de atributos
dat5 <- subset(dat3, select = -c(tcm,Fecha,date_dia,horadia)) # Drop column (variable) Date and Time from dat4 dataframe
#str(dat5)
# Modificación del lugar de los atributos en la tabla
dat5 <-dat5[,c(1,8,9,10,11,12,13,3,2,7,5,4,6)]
par(mfrow=c(1,3))
boxplot(dat5$PM2.5)
boxplot(dat5$PM10)
plot(dat5$PM10,dat5$PM2.5)
Cálculo de valores extremos con valores Z estandarizados
Resumen de resultados para cada una de las variables (acomodadas y de la 8 al 13) sumarios
Gráfica de valores normalizados (z) y sus histogramas
library(outliers)
z.scores <- dat5$PM2.5 %>% scores(type = "z") #Checking outliers for PM2.5 pollutant
z.scores %>% summary()Min. 1st Qu. Median Mean 3rd Qu. Max. -1.0472 -0.5889 -0.2833 0.0000 0.2259 21.5633
dat5$PM2.5zscore<-dat5$PM2.5 %>% scores(type = "z")
dat5$PM10zscore<-dat5$PM10 %>% scores(type = "z")
dat5$COzscore<-dat5$CO %>% scores(type = "z")
dat5$NO2zscore<-dat5$NO2 %>% scores
dat5$O3zscore<-dat5$O3 %>% scores(type = "z")
dat5$SO2zscore<-dat5$SO2 %>% scores(type = "z")
for(i in 8:13)
{
print(" ")
print("--------------",quote = FALSE)
colnames(dat5[i])
z.scores <- dat5[i] %>% scores(type = "z") #Checking outliers for PM2.5 pollutant
print(z.scores %>% summary())
a<-sum(abs(z.scores)>3)
fract_outl <-a/dim(dat5)[1]
#print(a/dim(dat5)[1])
cat("Fraction of outliers:",fract_outl,"\n")
} [1] ” ” [1] ————– PM2.5
Min. :-1.0472
1st Qu.:-0.5889
Median :-0.2833
Mean : 0.0000
3rd Qu.: 0.2259
Max. :21.5633
Fraction of outliers: 0.01392067 [1] ” ” [1] ————– PM10
Min. :-1.0965
1st Qu.:-0.5940
Median :-0.2754
Mean : 0.0000
3rd Qu.: 0.3127
Max. :42.8824
Fraction of outliers: 0.01036105 [1] ” ” [1] ————– CO
Min. :-0.9495
1st Qu.:-0.4727
Median :-0.2536
Mean : 0.0000
3rd Qu.: 0.1072
Max. :34.4109
Fraction of outliers: 0.0151284 [1] ” ” [1] ————– NO2
Min. :-1.1880
1st Qu.:-0.7309
Median :-0.3499
Mean : 0.0000
3rd Qu.: 0.4881
Max. : 5.5162
Fraction of outliers: 0.01150521 [1] ” ” [1] ————– O3
Min. :-1.46987
1st Qu.:-0.89930
Median :-0.01752
Mean : 0.00000
3rd Qu.: 0.70865
Max. : 8.69656
Fraction of outliers: 0.003051106 [1] ” ” [1] ————– SO2
Min. :-0.3464
1st Qu.:-0.3085
Median :-0.3085
Mean : 0.0000
3rd Qu.:-0.3085
Max. : 5.5994
Fraction of outliers: 0.04487668
#** SUMMARY** The results indicates each calidadaire variables has a small fraction of outliers, less than 2%“.
#*
#*
# Normalize the data
normalized_data <- data.frame(
PM2.5_normalized = scale(dat5$PM2.5),
PM10_normalized = scale(dat5$PM10),
CO_normalized = scale(dat5$CO),
NO2_normalized =scale(dat5$NO2),
O3_normalized =scale(dat5$O3),
SO2_normalized =scale(dat5$SO2))
# Visualize the original and normalized data (Optional)
par(mfrow = c(2, 3)) # Arrange plots in a 2x3 grid
hist(dat5$PM2.5, main = "PM2.5", xlab = "PM2.5")
hist(normalized_data$PM2.5_normalized, main = "Normalized PM2.5", xlab = "PM2.5 (Normalized)")
par(mfrow = c(2, 3)) # Arrange plots in a 2x3 grid
hist(dat5$PM10, main = "PM10", xlab = "PM2.5")
hist(normalized_data$PM10_normalized, main = "Normalized PM10", xlab = "PM10 (Normalized)")
#Excluding Outliers
#An attempt was done to clean the dataset dat5 by excluding data with zscore>3; however due time constrains this was not completed.
par(mfrow=c(1,3))
boxplot(dat5$PM2.5)
boxplot(dat5$PM2.5zscore)
dat5_norm_clean<- dat5 %>% filter(abs(dat5$PM2.5zscore) <=3 )
boxplot(dat5_norm_clean$PM2.5zscore)
Identificación del índice de la calidad de aire, de acuerdo a el monitoreo de las partículas menores o iguales a los 2.5 micras, 10 y contaminaciónd CO, NO2, O3, y SO2
Cálculo en https://www.qld.gov.au/environment/management/monitoring/air/air-monitoring/air-quality-categories
Basado en: https://rpubs.com/Jamriska/394004
#https://www.qld.gov.au/environment/management/monitoring/air/air-monitoring/air-quality-categories
pol1_standard = 25 # [ug/m3] PM2.5 (24 hour avg) National Air NEPM Standard
pol2_standard = 50 # [ug/m3] PM10 (24 hour avg) National Air NEPM Standard
pol3_standard = 9 # [ppm] CO ( 8 hour avg) National Air NEPM Standard
pol4_standard = 0.08 # [ppm] NO2 ( 1 hour avg) National Air NEPM Standard
pol5_standard = 0.065 # [ppm] O3 ( 1 hour avg) National Air NEPM Standard
pol6_standard = .1 # [ppm] SO2 ( 1 hour avg) National Air NEPM Standard
#Calculation of AQI_pol
#The AQI was calculated using Formula: AQI_pol = (pol_reading/pol_standard)*100. Below is a script used to calculate AQI for each of the selected pollutants:
#Create a dataframe with AQI_po1, …, AQI_pol6
#For further processing conc values for each pollutant were assigned to new variables (later merged into a df)
pol1_conc<-dat5[,8] # pol1 PM2.5
pol2_conc<-dat5[,9] # pol2 PM10
pol3_conc<-dat5[,10] # pol3 CO
pol4_conc<-dat5[,11] # pol4 NO2
pol5_conc<-dat5[,12] # pol5 O3
pol6_conc<-dat5[,13] # pol6 SO2
plot(pol1_conc)
summary(pol1_conc)Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 9.00 15.00 20.56 25.00 444.00
# Esl cálculo para cada contaminante se obtiene el leído entre el estándar por 100
AQI_pol1<-(pol1_conc/pol1_standard)*100
plot(AQI_pol1)
AQI_pol2<-(pol2_conc/pol2_standard)*100
AQI_pol3<-(pol3_conc/pol3_standard)*100
AQI_pol4<-(pol4_conc/pol4_standard)*100
AQI_pol5<-(pol5_conc/pol5_standard)*100
AQI_pol6<-(pol6_conc/pol6_standard)*100
par(mfrow=c(1, 2))
boxplot(pol1_conc,pol2_conc,pol3_conc,pol4_conc,pol5_conc,pol6_conc)
boxplot(AQI_pol1,AQI_pol2,AQI_pol3,AQI_pol4,AQI_pol5,AQI_pol6)
#Min-Max Normalisation method
#Defined as: X - min(X))/(max(X) - min(X)) Define a normalization function
normalize <- function(x) {
return ((x - min(x)) / (max(x) - min(x)))
}
dat5_Norm <- as.data.frame(lapply(dat5[8:13], normalize))
#CONC_pol1_6 <-data.frame(cbind(pol1_conc,pol2_conc,pol3_conc,pol4_conc,pol5_conc,pol6_conc))
#colnames(CONC_pol1_6)[1] <- "PM2.5_(ug/m3)" # change the name of column
#colnames(CONC_pol1_6)[2] <- "PM10_(ug/m3))" # change the name of column
#colnames(CONC_pol1_6)[3] <- "CO_(ppm)" # change the name of column
#colnames(CONC_pol1_6)[4] <- "NO2_(ppm)" # change the name of column
#colnames(CONC_pol1_6)[5] <- "O3_(ppm)" # change the name of column
#colnames(CONC_pol1_6)[6] <- "Vis_(1/Mm)" # change the name of column
#head(CONC_pol1_6)
dat6 <-data.frame(cbind(dat5,AQI_pol1,AQI_pol2,AQI_pol3,AQI_pol4,AQI_pol5,AQI_pol6))
#colnames(dat6)
AQI_pol1_6_names <-c("AQI_PM2.5","AQI_PM10","AQI_CO","AQI_NO2","AQI_O3","AQI_SO2")
for(i in 1:6){
i_pol=i+19
colnames(dat6)[i_pol] <- AQI_pol1_6_names[i]
}
#colnames(dat6)
# Cálculo del índice de la calidad del aire
# Se identifica el valor máximo de cada uno de los índices de concentración calculados
AQI_pol1_6 <-data.frame(cbind(dat5,AQI_pol1,AQI_pol2,AQI_pol3,AQI_pol4,AQI_pol5,AQI_pol6))
#AQI_pol1_6$max <- apply(dat6[, 20:25], 1, max)
AQI_Overall <-AQI_pol1_6[, "max"] <- apply(dat6[, 20:25], 1, max) # a single column - vector representing max of 6 observation for all rows for
#class(AQI_Overall) # AIQ_Overall is numeric
#summary(AQI_Overall)
hist(AQI_Overall)
qqnorm(AQI_Overall);qqline(AQI_Overall, col = 2) # checking data normality using q-q plot
ln_AQI_Overall <- log(AQI_Overall) # natural log of AQI_Overall
hist(ln_AQI_Overall)
qqnorm(ln_AQI_Overall);qqline(ln_AQI_Overall, col = 2) # improved result-> closer to normal distribution
AQI_Overall_Levels <- cut(AQI_Overall, breaks=c(0,33,66,99,149,Inf),
labels=c("Muy bueno","Bueno","Razonable","Pobre","Malo"))
par(mfrow=c(1, 1))
hist(AQI_Overall)
dat7<-cbind(dat6,AQI_Overall, AQI_Overall_Levels)
#colnames(dat7)
plot(dat7$AQI_Overall,
main = "Indice de calidad del aire por hora para Torreón 2023-2025",
xlab = "Índice horario",
ylab = "AQI [-]",
xlim=c(0, 9000),
ylim=c(0, 600)
)
#minor.tick(nx = 4, ny=5, tick.ratio=0.3)
legend("topright", inset=.05, title="AQI Range (upper bounds)",
c("Muy bueno","Bueno","Razonable","Pobre","Malo"), col=c("green", "blue", "yellow", "red"),
lty=c(2,2,2,2), lwd=c(3, 3,3,3),
cex = 0.77,
horiz=TRUE)
abline(h=c(33,66,99,149), col=c("green", "blue", "yellow", "red"), lty=c(2,2,2,2), lwd=c(3, 3,3,3))
Gráfica de las PM2.5 de acuerdo a los valores de la NORMA Oficial
Mexicana NOM-172-SEMARNAT-2023, Lineamientos para la obtención y
comunicación del índice de calidad del airey riesgos a la salud. https://www.dof.gob.mx/nota_detalle.php?codigo=5715154&fecha=25/01/2024#gsc.tab=0
dat7$buena<-15
dat7$acepta<-41
dat7$mala<-79
dat7$muymala<-130
dat7$extmala<-400
custom <- colorRampPalette(c("green","yellow", "orange", "red","purple"))
dat7<-cbind(dat6,AQI_Overall, AQI_Overall_Levels)
dat7$buena<-15
dat7$acepta<-41
dat7$mala<-79
dat7$muymala<-130
dat7$extmala<-400
#dat7<-dat7 %>% filter(Year==2024 & month == "3")
dat7<-dat7 %>% filter(Year==2024 & month == "1" & day==1)
a<-ggplot(dat7,
# aes(x=caliairehoradia$Fecha,
aes(x=dat7$tc, y=dat7$PM2.5))+
geom_line()+
geom_line(aes(x=dat7$tc, y=dat7$buena),color="green")+
geom_line(aes(x=dat7$tc, y=dat7$acepta),color="yellow")+
geom_line(aes(x=dat7$tc, y=dat7$mala),color="orange")+
geom_line(aes(x=dat7$tc, y=dat7$muymala),color="red")+
geom_line(aes(x=dat7$tc, y=dat7$extmala),color="purple")+
#geom_point()+ geom_line() +
theme_void()+
theme(text = element_text(size=10),
axis.text.x = element_text(angle=90, hjust=1)) +
ggtitle(paste0("PM2.5 por hora del día ","1 de enero 2024")) +
xlab("hora") + ylab("PM2.5 u/m3")
#https://stackoverflow.com/questions/62983435/how-do-i-using-scale-y-datetime-breaks-time-1-hours
b<- a +scale_x_datetime("Hora", breaks = scales::date_breaks('1 hour'),
labels=scales::date_format('%H:%M'))
b+theme_aire()
#https://www.rpubs.com/tomm1116/543579
#https://rpubs.com/HetalS/871233
#https://rpubs.com/HetalS/871233
#https://davidcarslaw.com/files/openairmanual.pdf
#https://archive.ics.uci.edu/datasets
#https://www.pranaair.com/blog/what-is-air-quality-index-aqi-and-its-calculation/
lincrit<-c("buena","acepta","mala","muy mala","ext mala")
lincolor<-c("green","yellow","orange","red","purple")
linpm2.5<-c(25,45,79,147,400)
datlimi<-data.frame(lincrit,lincolor,linpm2.5)
dat7<-cbind(dat6,AQI_Overall, AQI_Overall_Levels)
dat7$buena<-15
dat7$acepta<-41
dat7$mala<-79
dat7$muymala<-130
dat7$extmala<-400
#dat7<-dat7 %>% filter(Year==2024 & month == "1" & day==1)
#dat7<-dat7 %>% filter(Year==2024 & month == "1")
dat7<-dat7 %>% filter(Year==2024)
dfs<-dat7 %>%
mutate(date = lubridate::ymd_hms(tc),
date_dia = format(date, "%Y-%m-%d")) %>%
group_by(date_dia)%>%
summarize(promedio = mean(PM2.5,tr=0.1,na.rm=TRUE))%>%
mutate(linpmbuena=25,linpmacepta=45,linpmmala=79,linpmmuymala=147,linpmextmala=200)
dfs$date_dia<-as.POSIXct(dfs$date_dia,format = "%Y-%m-%d")
dfs$id<-as.numeric(row.names(dfs))
a<-dfs %>% ggplot(aes(x=date_dia,y=promedio))+
#geom_point()
geom_line()+
geom_line(aes(x=dfs$date_dia, y=dfs$linpmbuena),color="green")+
geom_area(aes(x=dfs$date_dia, y=dfs$linpmbuena),fill="green",alpha=.1)+
geom_ribbon(aes(x=dfs$date_dia, ymax=dfs$linpmacepta, ymin=dfs$linpmbuena), fill="yellow", alpha=.1) +
geom_line(aes(x=dfs$date_dia, y=dfs$linpmacepta),color="yellow")+
geom_ribbon(aes(x=dfs$date_dia, ymax=dfs$linpmmala, ymin=dfs$linpmacepta), fill="orange", alpha=.1) +
geom_line(aes(x=dfs$date_dia, y=dfs$linpmmala),color="orange")+
geom_ribbon(aes(x=dfs$date_dia, ymax=dfs$linpmmuymala, ymin=dfs$linpmmala), fill="red", alpha=.1) +
geom_line(aes(x=dfs$date_dia, y=dfs$linpmmuymala),color="red")+
#geom_line(aes(x=as.Date(dfs$date_dia), y=dfs$linpmextmala),color="purple")+
#theme_void()+
theme(text = element_text(size=10),
axis.text.x = element_text(angle=90, hjust=1)) +
ggtitle(paste0("PM2.5 promedio por día ","año 2024")) +
labs(x=id)+
xlab("Día ") + ylab("PM2.5 u/m3")
b<- a +scale_x_datetime("Día", breaks = scales:: breaks_pretty(12))
b<- a +scale_x_datetime("Día", breaks = scales:: breaks_pretty(24),labels=label_date_short())
b+theme_aire()
dat7<-cbind(dat6,AQI_Overall, AQI_Overall_Levels)
dat7$buena<-15
dat7$acepta<-41
dat7$mala<-79
dat7$muymala<-130
dat7$extmala<-400
#dat7<-dat7 %>% filter(Year==2024 & month == "1" & day==1)
#dat7<-dat7 %>% filter(Year==2024 & month == "1")
dat7<-dat7 %>% filter(Year==2024)
dfs<-dat7 %>%
mutate(date = lubridate::ymd_hms(tc),
date_mes = substring(date, 1, 7))%>%
group_by(date_mes)%>%
summarize(promedio = mean(PM2.5,tr=0.1,na.rm=TRUE))%>%
mutate(linpmbuena=15,linpmacepta=41,linpmmala=79,linpmmuymala=130,linpmextmala=400)
dfs$id<-as.numeric(row.names(dfs))
dfs$date_mes<-as.Date(paste0(dfs$date_mes,"-15"))
dfs$date_mes<-as.POSIXct(dfs$date_mes,format = "%Y-%m")
a<-dfs %>% ggplot(aes(x=date_mes,y=promedio))+
#geom_point()
geom_line()+
geom_line(aes(x=date_mes, y=dfs$linpmbuena),color="green")+
geom_area(aes(x=date_mes, y=dfs$linpmbuena),fill="green",alpha=.1) +
geom_ribbon(aes(x=date_mes, ymax=dfs$linpmacepta, ymin=dfs$linpmbuena), fill="yellow", alpha=.1) +
geom_line(aes(x=date_mes, y=dfs$linpmacepta),color="yellow")+
geom_ribbon(aes(x=date_mes, ymax=dfs$linpmmala, ymin=dfs$linpmacepta), fill="orange", alpha=.1) +
geom_line(aes(x=date_mes, y=dfs$linpmmala),color="orange")+
geom_ribbon(aes(x=date_mes, ymax=dfs$linpmmuymala, ymin=dfs$linpmmala), fill="red", alpha=.1) +
geom_line(aes(x=date_mes, y=dfs$linpmmuymala),color="red")+
#geom_line(aes(x=as.Date(dfs$date_dia), y=dfs$linpmextmala),color="purple")+
#theme_void()+
theme(text = element_text(size=10),
axis.text.x = element_text(angle=90, hjust=1)) +
ggtitle(paste0("PM2.5 promedio por mes ","año 2024")) +
labs(x=id)+
xlab("mes") + ylab("PM2.5 u/m3")
#a
#b<- a +scale_x_continuous("Mes", breaks = scales::breaks_width(5))
#b<- a +scale_x_datetime("Mes", breaks = scales:: breaks_pretty(12))
#b+theme_aire()
b<- a +scale_x_datetime("Mes", breaks = scales:: breaks_pretty(12),labels=label_date_short())
b+theme_aire()
No hay datos suficientes para generar series de tiempo Se solicitó información de la estación desde el año 2016, y respondieron que sólo contaban desde octubre 2023
#https://www.r-bloggers.com/2016/12/air-quality-in-indian-cities/
#análisi tiempo no hay suficientes periodos
dat7<-cbind(dat6,AQI_Overall, AQI_Overall_Levels)
dat7$buena<-15
dat7$acepta<-41
dat7$mala<-79
dat7$muymala<-130
dat7$extmala<-400
#dat7<-dat7 %>% filter(Year==2024 & month == "1" & day==1)
#dat7<-dat7 %>% filter(Year==2024 & month == "1")
#dat7<-dat7 %>% filter(Year==2024)
dfs<-dat7 %>%
mutate(date = lubridate::ymd_hms(tc),
date_mes = substring(date, 1, 7))%>%
group_by(date_mes)%>%
summarize(promedio = mean(PM2.5,tr=0.1,na.rm=TRUE))%>%
ungroup()
pm25_monthly_ts <- dfs %>%
select(-date_mes) %>%
map(function(x){ts(x, start = c(2023, 10), frequency = 12)}) %>%
do.call(cbind, .)
monthplot(pm25_monthly_ts)
library(ggseas)
pm<-pm25_monthly_ts %>%
as.data.frame() %>%
mutate(time = time(pm25_monthly_ts),
month = time - floor(time))
colnames(pm)[1] <- "valor"
acf(pm)
ap_df <- tsdf(pm25_monthly_ts)
head(ap_df) x y1 2023.750 13.83817 2 2023.833 22.92935 3 2023.917 26.75482 4 2024.000 27.28371 5 2024.083 14.85047 6 2024.167 11.64974
ggplot(ap_df, aes(x = x, y = y)) +
geom_line(colour = "grey75") 
ggplot(ap_df, aes(x = x, y = y)) +
stat_index(index.ref = 1:10)
ggplot(ap_df, aes(x = x, y = y)) +
stat_index(index.ref = 10, index.basis = 200)
ggplot(ap_df, aes(x = x, y = y)) +
geom_line(colour = "grey75") +
stat_rollapplyr(width = 12, align = "center") +
labs(x = "", y = "PM2.5)")
ggplot(ap_df, aes(x = x, y = y)) +
geom_line(colour = "grey75") +
stat_rollapplyr(width = 12, align = "center", geom = "point",
size = 0.5, colour = "steelblue") +
labs(x = "", y = "pm2.5)")
ggplot(ap_df, aes(x = x, y = y)) +
geom_line(colour = "grey50") +
stat_seas(colour = "blue") +
stat_stl(s.window = 7, colour = "red")
#ggsdc(ap_df, aes(x = x, y = y), method = "seas") + geom_line()