calidadaire

Contexto geográfico

Incorporación de la zona municipal de Torreón en el polígono urbano Uso de consultas para la lectura filtrada del Shapefile del contexto estatal filtrado municipal

Fuentes de datos

https://inegi.org.mx/contenidos/productos/prod_serv/contenidos/espanol/bvinegi/productos/geografia/marcogeo/794551132173_s.zip

Calidad de aire https://www.kaggle.com/datasets/elianaj/mexico-air-quality-dataset

ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/Marcogeoestadisticonacional/mg_2023_integrado/conjunto_de_datos/"
infile <- "00ent.shp"
nomarchi<-paste0(ruta,infile)
datageom<- read_sf(nomarchi,options = "ENCODING=latin1")
datageom= st_transform(datageom, crs = 4326)


ruta<-"C:/Users/cguer/Documents/Claudia/Midropbox/Investigacion y escritos/karamanis/data/"
infile <- "stations_rsinaica.CSV"
nomarchi<-paste0(ruta,infile)

#*stations_hourly.CSV
#*stations_daily.CSV
estaciones<-read.csv(nomarchi)
estaciones <- estaciones[which(!is.na(estaciones$lat) & !is.na(estaciones$lon)),]
estaciones_sf <- st_as_sf(estaciones, crs = 4326,
  coords = c("lon", "lat"))

ggplot() + geom_sf(data = datageom) + 
    geom_sf(data = estaciones_sf, color="red")

Información por hora

Se selecciona NO2, se filtra a los datos posteriores al 01 de enero de 2021, se eliminan las estaciones sin datos o con valores de cero Se grafica el promedio por estación

infile <- "stations_hourly.CSV"
nomarchi<-paste0(ruta,infile)
medidashor<-read.csv(nomarchi)

medidashor2 <- medidashor[, c("datetime", "station_id", "NO2")]
class(medidashor2)

[1] “data.frame”

Sys.setenv(TZ = "UTC") # don't use local time zone
medidashor2$t<-as.POSIXct(medidashor2$datetime)
medidashor2 <-medidashor2[order(medidashor2$t),]

medidashor2 <-medidashor2%>% 
  filter(t > ymd_hms("2021-01-01 01:00:00"))  |> dplyr::filter(!is.na(NO2) & NO2 != 0) 

medidashor2st<-right_join(medidashor2,estaciones_sf)
class(medidashor2st)

[1] “data.frame”

medidashor2st <-medidashor2st%>% 
   dplyr::filter(!is.na(NO2)) 

medidashor_sf <- st_as_sf(medidashor2st, crs = st_crs(4326))
class(medidashor_sf)

[1] “sf” “data.frame”

#ggplot() + geom_sf(data = datageom) + 
#    geom_sf(data = medidashor_sf, aes(color = NO2))



s <- split(medidashor_sf$NO2, medidashor_sf$station_id)
sf<-lapply(s, function(x) mean(x, na.rm = T))
library(dplyr)


stackdf <- do.call(rbind, sf)
class(stackdf)

[1] “matrix” “array”

library(tidyverse)

dfs <- as.data.frame(stackdf) %>% 
  rownames_to_column() 

colnames(dfs) <- c("station_id", "NO2a")
dfs$station_id<-as.numeric(dfs$station_id)
mediasNO2<-left_join(dfs,estaciones_sf)
class(mediasNO2)

[1] “data.frame”

mediasNO2_sf <- st_as_sf(mediasNO2, crs = st_crs(4326))
class(mediasNO2_sf)

[1] “sf” “data.frame”

ggplot() + geom_sf(data = datageom) + 
    geom_sf(data = mediasNO2_sf, aes(color = NO2a))

Se prepara la grid para la interpolación y se utiliza la interpolación del peso de la distancia inversa

crs <- st_crs("EPSG:32632")
datageom<-st_transform(datageom,crs)
library(stars) |> suppressPackageStartupMessages()
st_bbox(datageom) |>
  st_as_stars(dx = 10000) |>
  st_crop(datageom) -> grd
grd

stars object with 2 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s values 0 0 0 0 0 0 234895 dimension(s): from to offset delta refsys x/y x 1 633 -11120618 10000 WGS 84 / UTM zone 32N [x] y 1 500 16570602 -10000 WGS 84 / UTM zone 32N [y]

mediasNO2_sf<-st_transform(mediasNO2_sf,crs)
library(gstat)
i <- idw(NO2a~1, mediasNO2_sf, grd)

[inverse distance weighted interpolation]

# [inverse distance weighted interpolation]


ggplot() + geom_stars(data = i, 
                      aes(fill = var1.pred, x = x, y = y)) + 
    xlab(NULL) + ylab(NULL) +
    geom_sf(data = st_cast(datageom, "MULTILINESTRING")) + 
    geom_sf(data = mediasNO2_sf)

crs = st_crs(4326)
datageom<-st_transform(datageom,crs)
mediasNO2_sf<-st_transform(mediasNO2_sf,crs)
i<-st_transform(i,crs)


ggplot() + geom_stars(data = i, 
                      aes(fill = var1.pred)) + 
    xlab(NULL) + ylab(NULL) +
    geom_sf(data = st_cast(datageom, "MULTILINESTRING")) + 
    geom_sf(data = mediasNO2_sf)

Se trabaja el archivo de las mediciones diarias con el parámetro de PM2.5 Valores posteriores al 01 de enero del 2014

Se contabilizaron otros parámetros (sólo se grafica el PM2.5) Rango PM2.5 mayores a 30 y menores de 250

infile <- "stations_daily.CSV"
nomarchi<-paste0(ruta,infile)
medidasdia<-read.csv(nomarchi)
class(medidasdia)

[1] “data.frame”

Sys.setenv(TZ = "UTC") # don't use local time zone
medidasdia$t<-as.POSIXct(medidasdia$datetime)
medidasdia <-medidasdia[order(medidasdia$t),]

medidasdia  <-medidasdia%>% 
  filter(t > ymd_hms("2014-01-01 01:00:00"))
library(dplyr)


b<-medidasdia %>% dplyr::filter(!is.na(PM10) & PM10 != 0) %>% 
  group_by(PM10) %>% 
  summarise(n = n())


c<-medidasdia %>%  dplyr::filter(!is.na(NOx) & NOx != 0) %>% 
  group_by(NOx) %>% 
  summarise(n = n())


d<-medidasdia %>% dplyr::filter(!is.na(O3) & O3 != 0) %>% 
  group_by(O3) %>% 
  summarise(n = n())


e<-medidasdia %>% dplyr::filter(!is.na(CO) & CO!= 0) %>% 
  group_by(CO) %>% 
  summarise(n = n())


f<-medidasdia %>%  dplyr::filter(!is.na(HR) & HR != 0) %>% 
  group_by(HR) %>% 
  summarise(n = n())

g<-medidasdia %>% dplyr::filter(!is.na(NO) & NO != 0) %>% 
  group_by(NO) %>% 
  summarise(n = n())


h<-medidasdia %>% dplyr::filter(!is.na(NO2) & NO2!= 0) %>% 
  group_by(NO2) %>% 
  summarise(n = n())


i<-medidasdia %>%  dplyr::filter(!is.na(TMP) & TMP != 0) %>% 
  group_by(TMP) %>% 
  summarise(n = n())

medidasdia<-medidasdia %>% dplyr::filter(!is.na(PM2.5) & (PM2.5 > 50 & PM2.5 <= 9.9e+02)) 
a<-medidasdia %>% 
  group_by(PM2.5) %>% 
  summarise(n = n())

s <- split(medidasdia$PM2.5, medidasdia$station_id)
sf<-lapply(s, function(x) mean(x, na.rm = T))
library(dplyr)


stackdf <- do.call(rbind, sf)
class(stackdf)

[1] “matrix” “array”

library(tidyverse)

dfs <- as.data.frame(stackdf) %>% 
  rownames_to_column() 

colnames(dfs) <- c("station_id", "PM2.5a")
dfs$station_id<-as.numeric(dfs$station_id)
mediasNO2<-left_join(dfs,estaciones_sf)
class(mediasNO2)

[1] “data.frame”

mediasNO2<-mediasNO2 |> dplyr::filter(!is.na(PM2.5a)) 
mediasNO2_sf <- st_as_sf(mediasNO2, crs = st_crs(4326))
class(mediasNO2_sf)

[1] “sf” “data.frame”

ggplot() + geom_sf(data = datageom) + 
    geom_sf(data = mediasNO2_sf, aes(color = PM2.5a))

crs <- st_crs("EPSG:32632")
datageom<-st_transform(datageom,crs)
library(stars) |> suppressPackageStartupMessages()
st_bbox(datageom) |>
  st_as_stars(dx = 20000) |>
  st_crop(datageom) -> grd
grd

stars object with 2 dimensions and 1 attribute attribute(s): Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s values 0 0 0 0 0 0 58824 dimension(s): from to offset delta refsys x/y x 1 317 -11120618 20000 WGS 84 / UTM zone 32N [x] y 1 250 16570602 -20000 WGS 84 / UTM zone 32N [y]

mediasNO2_sf<-st_transform(mediasNO2_sf,crs)
library(gstat)
i <- idw(PM2.5a~1, mediasNO2_sf, grd)

[inverse distance weighted interpolation]

# [inverse distance weighted interpolation]


ggplot() + geom_stars(data = i, 
                      aes(fill = var1.pred, x = x, y = y)) + 
    xlab(NULL) + ylab(NULL) +
    geom_sf(data = st_cast(datageom, "MULTILINESTRING")) + 
    geom_sf(data = mediasNO2_sf)

crs = st_crs(4326)
datageom<-st_transform(datageom,crs)
mediasNO2_sf<-st_transform(mediasNO2_sf,crs)
i<-st_transform(i,crs)


ggplot() + geom_stars(data = i, 
                      aes(fill = var1.pred)) + 
    xlab(NULL) + ylab(NULL) +
    geom_sf(data = st_cast(datageom, "MULTILINESTRING")) + 
    geom_sf(data = mediasNO2_sf)