project_workflows/IndonesiaIESR/plot_video_IESR_using_functions.R

In energyandcleanair/creapuff: CREAPUFF: R package to work with CALMET and CALPUFF

require(terra)
require(sf)
require(tidyverse)
require(magrittr)
require(lubridate)
library(readxl)
library(creapuff) 
#list.files(path='R', full.names=T) %>% sapply(source)
require(rcrea)
require(creahelpers)


# Parameters ###################################################################
# ============================= Project specific ===============================
project_dir="G:/IndonesiaIESR"       # calpuff_external_data-2 persistent disk (project data)
gis_dir <- "H:/gis"                    # The folder where we store general GIS data

#load(file.path(project_dir, 'CALPUFF_setup.RData'))

input_dir <- file.path(project_dir,"calpuff_suite") # Where to read all CALPUFF generated files
output_dir <- file.path(project_dir, 'animation'); dir.create(output_dir)
emissions_dir <- file.path(project_dir,"emissions")

select <- dplyr::select
source('project_workflows/read_IESR_emissions.R')
emis %>% distinct(cluster, .keep_all=T) %>% to_spdf %>% st_as_sf -> point_sources_to_plot


#load grid parameters
calmet_result <- readRDS(file.path(input_dir,"calmet_result.RDS" ))
UTMZ <- calmet_result$params[[01]]$IUTMZN
UTMH <- calmet_result$params[[01]]$UTMHEM

sf_use_s2(F)

calpuff_files <- get_calpuff_files(ext=".csv", gasunit = 'ug', dir=input_dir, hg_scaling=1e-3)
grids = get_grids_calpuff(calpuff_files, UTMZ, UTMH, map_res=5)

plot_bb = point_sources_to_plot %>% st_transform(crs=crs(grids$gridR)) %>% extent() %>% add(800)
grids$gridR %<>% crop(plot_bb)

scenarios_to_process=calpuff_files$scenario %>% 
  (function(x) paste0('v', substr(x,2,8))) %>% unique

emissions_data %<>% mutate(video_name = paste0('v', substr(tolower(emission_names),2,8)))


#list and select files to process
list.files(path=input_dir, pattern='tser.*1hr.*\\.dat', full.names = T) %>% 
  subset(grepl('pm25_1hr', .)) -> infiles

infiles %<>% grep(paste(paste0('_', scenarios_to_process, '\\.'), collapse='|'), ., value=T)

#convert tseries.dat files to RDS for faster reading and processing
tseries_to_rds(infiles) -> infiles

#select time steps to process
times <- read_times_from_rds(infiles) %>% subset(month(.)==4 & hour(.) %in% ((0:7)*3))

infiles <- infiles[order(as.numeric(str_match(infiles, "[0-9]+\\.")))]

#interpolate concentration time series files with receptor data into gridded rasters
tseries_rds_to_raster(infiles, grids, times=times)

#sum concentrations for multiple simulations for each timestep
sum_tseries_rasters(infiles, case_name='allsources', times=times, overwrite=F)



# Plots ########################################################################
target_crs <- crs(grids$gridR)

#get basemap
require(ggspatial); require(ggmap); require(ggrepel)
register_google(readLines("~/google_api_key.txt"))

plot_bb = point_sources_to_plot %>% st_transform(crs=4326) %>% extent() %>% add(2)
basemap <- get_basemap(plot_bb, zoom=4)

plot_bb %>% as('SpatialPolygons') %>% st_as_sf %>% st_set_crs(4326) %>% st_transform(3857) -> plot_bb_polygon
ggmap(basemap) + layer_spatial(plot_bb_polygon, fill=NA, linewidth=2, color='red')

#query the files to plot
calpuff_files_all <- get_calpuff_files(filetype='tseries', dir=output_dir) %>% 
  filter(scenario=='allsources') %>% 
  mutate(title=paste(speciesName, 'concentration from coal power plants in Indonesia'),
         subtitle=format(datetime, "%b %d, %H:%M"))

#output plots
calpuff_files_all %>% group_by(datetime) %>% 
  plot_video_frames(plot_bb=plot_bb, 
                    point_sources=point_sources_to_plot,
                    basemap=basemap,
                    output_dir=output_dir,
                    quicksave_options=list(width=10, height=6, scale=1.33, logo=T),
                    ggplot_theme=theme(legend.position='bottom', legend.direction = "horizontal"),
                    contour_guide=guide_legend(nrow=1))