Skip to content

Commit

Permalink
Merge pull request #15 from bcgov/update-2019
Browse files Browse the repository at this point in the history 
Update 2021
  • Loading branch information
@KarHarker
KarHarker committed Apr 13, 2021
2 parents 562d582 + c6d5af1 commit 1f4791f
Show file tree
Hide file tree
Showing 8 changed files with 300 additions and 188 deletions.
90 changes: 85 additions & 5 deletions 02_clean.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,15 @@ stn_names <- read_csv("data/stn_names_reporting.csv") %>%
mutate(ems_id = str_pad(ems_id, 7, "left", "0")) %>%
rename(orig_stn_name = station_name)

max_year <- 2017
max_year <- 2018

## Set stations to exclude from analyis (those in indsutrial settings):
excluded_stations <- stations$EMS_ID[grepl("industr", stations$STATION_ENVIRONMENT, ignore.case = TRUE)]

# remove bc_hydro stations (fort st James) as not reported
bchydro_stations <- unique(stations$EMS_ID[grepl("INDUSTRY-BCH", stations$STATION_OWNER, ignore.case = TRUE)])
excluded_stations <- unique(c(excluded_stations, bchydro_stations))

## Exclude Kitimat Smeltersite which is industrial but not labelled as such in
## the stations metadata
excluded_stations <- unique(c(excluded_stations, "E290529"))
Expand All @@ -47,6 +51,7 @@ combo_squamish_id <- paste(squamish_ems_ids, collapse = "-")
## Met == meteorological stns using Campbell loggers;
## BAM == Beta Attenuation Monitoring for PM measurement.
select_pattern <- "_60$|Met$|OLD$|BAM$|Squamish Gov't Bldg"

stations_clean <- rename_all(stations, tolower) %>%
mutate(ems_id = case_when(ems_id %in% squamish_ems_ids ~ combo_squamish_id,
TRUE ~ gsub("_.+$", "", ems_id))) %>%
Expand All @@ -67,7 +72,7 @@ stations_clean <- rename_all(stations, tolower) %>%
station_id = "ems_id",
coords = c("longitude", "latitude"))

## Format dates, extract 2014-2016, set variable names.
## Format dates, extract time period , set variable names.

pm25 <- pm25_all %>%
filter(!EMS_ID %in% excluded_stations) %>%
Expand Down Expand Up @@ -109,10 +114,58 @@ pm25 <- pm25_all %>%
) %>%
distinct()


# temp : check stations with NA's for instrument type
na_stations <- pm25 %>%
filter(is.na(instrument_type)) %>%
group_by(station_name, year)%>%
summarise(count = n())
na_stations

## Plot deployments of different instruments at each station

plot_station_instruments(pm25)
plot_station_instruments(pm25, instrument = "instrument_type")

library(ggplot2)

# Temporary check outputs -------------------------------------------------
# follow up with checks on instrument type.
st_nms <- unique(pm25$station_name)
st_nms <- st_nms[1:25]
pm25a <- pm25 %>%
filter(station_name %in% st_nms)
p1_25 <- plot_station_instruments(pm25a)
p1_25
ggsave( "tmp/pm25_p1_25.jpg", plot = last_plot())

st_nms <- unique(pm25$station_name)
st_nms <- st_nms[26:50]
pm25a <- pm25 %>%
filter(station_name %in% st_nms)
p26_50 <- plot_station_instruments(pm25a)
p26_50
ggsave( "tmp/pm25_p26_50.jpg", plot = last_plot())
#
st_nms <- unique(pm25$station_name)
st_nms <- st_nms[51:75]
pm25a <- pm25 %>%
filter(station_name %in% st_nms)
p51_75 <- plot_station_instruments(pm25a)
p51_75
ggsave( "tmp/pm25_p51_75.jpg", plot = last_plot())
#
st_nms <- unique(pm25$station_name)
st_nms <- st_nms[76:100]
pm25a <- pm25 %>%
filter(station_name %in% st_nms)
p76_100 <- plot_station_instruments(pm25a)
p76_100
ggsave( "tmp/pm25_p76_100.jpg", plot = last_plot())
#

##------------------------------------------------------------------------

## Summarise the dates during the most recent three years (caaqs timeframe)
## that different PM2.5 monitoring instrument types were deployed at each
## station so we can get the most data
Expand All @@ -126,14 +179,14 @@ instrument_deployments <- mutate(pm25, date = as.Date(date_time)) %>%
n_days = n()) %>%
ungroup()

## Select the monitor at each station that hast the most days
## Select the monitor at each station that has the most days
max_deployment_by_station <- group_by(instrument_deployments, ems_id) %>%
summarise(which_instrument = instrument_type[which.max(n_days)])

squamish <- filter(
pm25, instrument_type == "FEM",
(ems_id == squamish_ems_ids[1] & year == 2015) |
(ems_id == squamish_ems_ids[2] & year %in% 2016:2017)
(ems_id == squamish_ems_ids[2] & year %in% 2016:2018)
) %>%
# If BAM and SHARP were running in concert - take BAM
group_by(date_time) %>%
Expand All @@ -153,7 +206,7 @@ pm25_clean <- pm25 %>%
inner_join(select(stations_clean, ems_id, station_name),
by = "ems_id") %>%
# tsibble time series package to make sure hourly data
as_tsibble(key = id(ems_id, station_name, instrument, instrument_type),
as_tsibble(key = c(ems_id, station_name, instrument, instrument_type),
regular = FALSE) %>%
group_by(ems_id, station_name, instrument, instrument_type) %>%
index_by(date_hr = ceiling_date(date_time, "hour") - 1) %>%
Expand All @@ -168,3 +221,30 @@ plot_station_instruments(pm25_clean, instrument = "instrument_type")

save(pm25_clean, stations_clean, max_year, file = "tmp/pm25_clean.rda")


# Temporary check outputs -------------------------------------------------
# follow up with checks on instrument type.
st_nms <- unique(pm25_clean$station_name)
st_nms <- st_nms[1:25]
pm25a <- pm25_clean %>%
filter(station_name %in% st_nms)
p1_25 <- plot_station_instruments(pm25a)
p1_25
ggsave( "tmp/pm25c_p1_25.jpg", plot = last_plot())

st_nms <- unique(pm25_clean$station_name)
st_nms <- st_nms[26:50]
pm25a <- pm25_clean %>%
filter(station_name %in% st_nms)
p26_50 <- plot_station_instruments(pm25a)
p26_50
ggsave( "tmp/pm25c_p26_50.jpg", plot = last_plot())

st_nms <- unique(pm25_clean$station_name)
st_nms <- st_nms[51:75]
pm25a <- pm25_clean %>%
filter(station_name %in% st_nms)
p51_75 <- plot_station_instruments(pm25a)
p51_75
ggsave( "tmp/pm25c_p51_75.jpg", plot = last_plot())

30 changes: 25 additions & 5 deletions 03_analysis.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,18 @@ excludes <- filter(get_daily(pm25_caaqs_24h), avg_24h > 28,
between(month(date), 5, 9)) %>%
select(site_group_vars, date)


# output for Rmd file # 61 stations with
no.site.with.extreme.fire <- filter(get_daily(pm25_caaqs_24h), avg_24h > 28,
between(month(date), 5, 9),
between(year(date),2016,2018)) %>%
select(site_group_vars, date) %>%
ungroup() %>%
select(ems_id) %>%
distinct %>%
pull


# PM25 24 Hour - Management -----------------------------------------------

pm25_24h_caaqs_mgmt <- caaqs_management(pm25_caaqs_24h, exclude_df = excludes,
Expand Down Expand Up @@ -72,7 +84,6 @@ pm_annual_caaqs_results <- station_results_pipeline(pm25_annual_caaqs_mgmt)
# Combo reporting df of annual and 24h caaqs results
pm_caaqs_combined_results <- bind_rows(pm_annual_caaqs_results, pm_24h_caaqs_results)


# Airzone results ---------------------------------------------------------

airzone_caaqs_pm_annual <- filter(pm_caaqs_combined_results, metric == "pm2.5_annual") %>%
Expand All @@ -89,15 +100,24 @@ az_ambient <- bind_rows(airzone_caaqs_pm24h, airzone_caaqs_pm_annual) %>%
az_mgmt <- az_ambient %>%
group_by(airzone) %>%
slice(which.max(mgmt_level)) %>%
mutate(caaqs_year = 2017L) %>%
mutate(caaqs_year = max_year) %>%
select(caaqs_year, airzone, mgmt_level, rep_metric = metric,
metric_value = metric_value_mgmt, rep_stn_id = rep_stn_id_mgmt,
rep_stn_name = station_name_mgmt)


stations.with.exlusion <- pm_24h_caaqs_results %>%
filter(ems_id %in% no.site.with.extreme.fire) %>%
select(ems_id) %>%
distinct()

no.stations.with.exclusion = length(stations.with.exlusion$ems_id)


save(list = ls(), file = "tmp/analysed.RData")

dir.create("out", showWarnings = FALSE)
write_csv(az_ambient, "out/pm2.5_airzone_results_2017.csv", na = "")
write_csv(pm_caaqs_combined_results, "out/pm2.5_caaqs_combined_results_2017.csv", na = "")
write_csv(az_mgmt, "out/pm2.5_airzone_management_levels_2017.csv", na = "")
write_csv(az_ambient, "out/pm2.5_airzone_results.csv" , na = "")
write_csv(pm_caaqs_combined_results, "out/pm2.5_caaqs_combined_results.csv", na = "")
write_csv(az_mgmt, "out/pm2.5_airzone_management_levels.csv", na = "")

15 changes: 11 additions & 4 deletions 04_output.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@ if (!exists("az_final")) load("tmp/analysed.RData")
# Create output directory:
dir.create("out", showWarnings = FALSE)

rep_year <- max_year

## @knitr pre

az <- st_intersection(airzones(), st_geometry(bc_bound())) %>%
Expand Down Expand Up @@ -129,12 +131,12 @@ names(stn_plots) <- emsids
for (emsid in emsids) {

# Create plots
p_24 <- plot_ts(pm25_caaqs_24h, id = emsid, , id_col = "ems_id",
rep_yr = 2017, plot_caaqs = TRUE, plot_exceedances = FALSE,
p_24 <- plot_ts(pm25_caaqs_24h, id = emsid, id_col = "ems_id",
rep_yr = max_year, plot_caaqs = TRUE, plot_exceedances = FALSE,
base_size = 14)

p_annual <- plot_ts(pm25_caaqs_annual, id = emsid, , id_col = "ems_id",
rep_yr = 2017, plot_caaqs = TRUE, plot_exceedances = FALSE,
p_annual <- plot_ts(pm25_caaqs_annual, id = emsid, id_col = "ems_id",
rep_yr = max_year, plot_caaqs = TRUE, plot_exceedances = FALSE,
base_size = 14)

p_annual <- p_annual + coord_cartesian(ylim = c(0, 80))
Expand Down Expand Up @@ -244,6 +246,7 @@ width <- 778
height <- 254

for (i in seq_along(stn_plots)) {
#i = 1
emsid <- names(stn_plots[i])
daily_plot <- stn_plots[[i]]$daily
annual_plot <- stn_plots[[i]]$annual
Expand All @@ -270,3 +273,7 @@ save(
pm_mgmt_chart,
file = "tmp/plots.RData"
)


save(list = ls(), file = "tmp/analysed.RData")

116 changes: 45 additions & 71 deletions 05_databc_output.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,83 +20,57 @@ library(readr)
library(dplyr)
library(bcmaps)
library(tidyr)
library(bcdata)
library(tidyverse)

if (!exists("az_final")) load("tmp/analysed.RData")

az <- st_intersection(airzones(), st_geometry(bc_bound())) %>%
group_by(airzone = Airzone) %>%
summarize()
if (!exists("az_final")) load("tmp/analysed.RData")
dir.create("out/databc", showWarnings = FALSE)

az_mgmt_sf <- az %>%
left_join(az_mgmt) %>%
mutate_at("mgmt_level", ~ replace_na(.x, "Insufficient Data"))
# station results

dir.create("out/databc", showWarnings = FALSE)
stations_old <- bcdc_get_data('699be99e-a9ba-403e-b0fe-3d13f84f45ab',
resource = 'bfa3fdd8-2950-4d3a-b190-52fb39a5ffd4')

stations_summary <- pm_caaqs_stations_all %>%
select(-c(flag_yearly_incomplete, flag_two_of_three_years,flag_daily_incomplete)) %>%
rename(latitude = lat, longitude = lon) %>%
mutate(caaqs_ambient = as.character(caaqs_ambient),
mgmt_level = as.character(mgmt_level))

# Station-level results
pm_2013 <- read_csv(soe_path("Operations ORCS/Indicators/air/fine_pm/2015/pm25_site_summary.csv")) %>%
select(-regional_district) %>%
rename_all(tolower) %>%
rename(station_name = display_name, instrument_type = monitor,
caaqs_year = caaq_year, metric_value_ambient = metric_value,
caaqs_ambient = caaqs, mgmt_level = mgmt) %>%
left_join(select(stations_clean, ems_id, city))

pm_2016 <- read_csv(soe_path("Operations ORCS/Indicators/air/fine_pm/2017/pm25_site_summary.csv")) %>%
rename_all(tolower) %>%
mutate(caaqs_year = 2016L) %>%
rename(metric_value_ambient = metric_value, caaqs_ambient = caaqs)

station_caaqs_combined <- pm_caaqs_combined_results %>%
select(-starts_with("flag")) %>%
rename(latitude = lat, longitude = lon) %>%
bind_rows(pm_2016, pm_2013) %>%
bind_rows(stations_old, stations_summary) %>%
arrange(caaqs_year) %>%
select(ems_id, station_name, city, longitude, latitude, instrument_type,
airzone, metric, caaqs_year, min_year, max_year, n_years,
metric_value_ambient, caaqs_ambient, excluded, metric_value_mgmt,
mgmt_level)

# Ambient airzone caaqs
az_ambient_2016 <- read_csv(soe_path("Operations ORCS/Indicators/air/fine_pm/2017/pm25_ambient_airzone_caaqs_summary.csv")) %>%
rename_all(tolower)

az_ambient_2016_tidy <- select(az_ambient_2016, airzone, contains("annual")) %>%
mutate(metric = "pm2.5_annual") %>%
rename_all(function(x) gsub("_annual", "", x)) %>%
bind_rows(select(az_ambient_2016, airzone, contains("24")) %>%
mutate(metric = "pm2.5_24h") %>%
rename_all(function(x) gsub("_24h", "", x))) %>%
mutate(caaqs_year = 2016L) %>%
select(airzone, metric, caaqs_year, metric_value = pm2.5_metric,
caaqs, rep_stn_id = rep_id, rep_stn_name = rep_stn, n_years)

az_ambient_combined <- az_ambient %>%
select(airzone, metric, contains("ambient")) %>%
rename_all(function(x) gsub("_ambient", "", x)) %>%
rename(rep_stn_name = station_name) %>%
# join to fill out airzones with missing values (e.g., NW)
right_join(select(az_ambient_2016_tidy, airzone, metric)) %>%
mutate(caaqs_year = 2017L) %>%
bind_rows(az_ambient_2016_tidy) %>%
replace_na(list(caaqs = "Insufficient Data")) %>%
select(airzone, metric, caaqs_year, everything()) %>%
arrange(caaqs_year, airzone, metric)

# Airzone management levels
az_mgmt_2016 <- read_csv(soe_path("Operations ORCS/Indicators/air/fine_pm/2017/pm25_mgmt_airzone_caaqs_summary.csv")) %>%
rename_all(tolower) %>%
rename_all(function(x) gsub("^caaq_mgmt_", "", x)) %>%
rename(mgmt_level = caaq_mgmt, rep_metric = metric) %>%
mutate(caaqs_year = 2016L)

az_mgmt_combined <- st_set_geometry(az_mgmt_sf, NULL) %>%
mutate(caaqs_year = 2017L) %>%
bind_rows(az_mgmt_2016) %>%
write_csv("out/databc/pm25sitesummary.csv", na = "")


# fine particulate matter - ambient caaqs by airzone

az_ambient_old <- bcdc_get_data('699be99e-a9ba-403e-b0fe-3d13f84f45ab',
resource = '5dd4fcf9-f7c6-49cf-90a0-0a5c9bc00334')

az_ambient <- az_ambient %>%
select(c(airzone, metric, n_years_ambient, metric_value_ambient,
caaqs_ambient, rep_stn_id_ambient, station_name_ambient)) %>%
rename_at(vars(ends_with("_ambient")), ~ gsub("_ambient", "",.)) %>%
mutate(caaqs = as.character(caaqs),
caaqs_year = max_year,
rep_stn_name = station_name) %>%
select(-station_name)

bind_rows(az_ambient_old, az_ambient) %>%
arrange(caaqs_year) %>%
replace_na(list(mgmt_level = "Insufficient Data"))
write_csv("out/databc/pm25-airzone-caaqs.csv", na = "")


# air management zone

az_mgmt_old <- bcdc_get_data('699be99e-a9ba-403e-b0fe-3d13f84f45ab',
resource = '700a7155-0b68-4e5a-bbe0-11d4b844ec57')

az_mgmt <- az_mgmt %>%
mutate(mgmt_level = as.character(mgmt_level))

write_csv(station_caaqs_combined, "out/databc/pm25sitesummary.csv", na = "")
write_csv(az_ambient_combined, "out/databc/pm25-airzone-caaqs.csv", na = "")
write_csv(az_mgmt_combined, "out/databc/pm25-airzone-management-levels.csv", na = "")
bind_rows(az_mgmt_old, az_mgmt) %>%
arrange(caaqs_year)%>%
write_csv("out/databc/pm25-airzone-management-levels.csv.csv", na = "")

19 changes: 0 additions & 19 deletions combine_prev_year.R
View file

This file was deleted.

Loading

0 comments on commit 1f4791f

Please sign in to comment.