## Original code from Kristen Foley, February 2022
##
## Modified by David S. Dixon on 3 April 2024 to remove dependencies on 
##   M3 package (removed by CRAN) and raster package (superceded by terra package).
##   Also, hires moved from the map package to the mapdata package.
##

# Download CMAQ data from the CMAS data warehouse:
#https://drive.google.com/drive/folders/1VnSxXL3MTovGuBAN__8T4nO3LsAkoB0T?usp=sharing
# Download file: HR2DAY_LST_ACONC_EQUATES_v532_12US1_2017.nc
#Location of file on atmos: /work/MOD3EVAL/fib/data_warehouse/ 
# DOI for EQUATES 2017 CMAQ output: https://doi.org/10.15139/S3/F2KJSK


library(ncdf4)
library(terra)
library(fields)
library(maps)
library(mapdata)
library(leaflet)
library(htmlwidgets)

#> Assumption buit into M3
EARTH_RADIUS <- 6370000

#>#################################################
#> Convenience function to decode the date integer 
#>#################################################
format_date <- function(int_date){
  return(as.Date(paste0(as.integer(int_date / 1000), "-01-01"), format="%Y-%m-%d"))
}

#>#################################################
#> Convenience function to decode the time integer
#>#################################################
format_time <- function(int_time){
  return(sub("(\\d\\d)(\\d\\d)(\\d\\d)", "\\1:\\2:\\3", sprintf("%06.0f",int_time)))
}

#> Set location where you have downloaded the CMAQ netcdf file.  
# setwd("/work/MOD3EVAL/fib/EQUATES/data_sharing/daily_files")
setwd("/media/ddixon/ext0/Data_A_E/EPA/EQUATES")

#> Name of CMAQ netcdf file. 
#cctm.file <- "HR2DAY_LST_ACONC_EQUATES_v532_12US1_2017.nc"
cctm.file <- "HR2DAY_LST_ACONC_EQUATES_v532_12US1_2008-001.nc"

#> Open the CCTM file.
cctm.in <- nc_open(cctm.file)

#> Print information about a netCDF file, including the variables and dimensions it contains.
print(cctm.in)

#> Create a list of all model variables in cctm.file. 
#> CMAQ netcdf files are formatted following I/O API conventions (https://www.cmascenter.org/ioapi/).  I/O API is a data storage library designed specifically for CMAQ data. 
#> A variable called “TFLAG” will always be the first variable listed in a CMAQ I/O API file, so remove the first element of the list.
all.mod.names <- unlist(lapply(cctm.in$var, function(var)var$name))[-1]

#> Create a list units for all the model variables in the cctm.file. 
#> A variable called “TFLAG” will always be the first variable listed in an I/O API file, so remove the first element of the list.
#> Use gsub to strip out extra spaces.
all.mod.units <- gsub(" ","",unlist(lapply(cctm.in$var, function(var)var$units))[-1])

#> Pull out the time steps and the grid coordinates associated with the data.
sdate <- ncatt_get(cctm.in, varid=0, attname="SDATE")$value
stime <- ncatt_get(cctm.in, varid=0, attname="STIME")$value
tstep <- ncatt_get(cctm.in, varid=0, attname="TSTEP")$value
time.series.length <- cctm.in$dim$TSTEP$len
# julian day minus one gives days *after* january first
start.date <- as.Date((sdate %% 1000) - 1, origin = format_date(sdate)) 
# start datetime
start.date.time <- strptime(paste0(as.character(start.date), " ", 
                                   format_time(stime)),
                            format="%Y-%m-%d %H:%M:%S", tz="GMT")
# timestep in seconds
timestep <- eval(parse(text = sub("(\\d\\d)(\\d\\d)(\\d\\d)", 
                                  "\\1*3600+\\2*60+\\3", 
                                  sprintf("%06.0f",tstep))))
# the time dimension in date format
date.seq <- start.date.time + timestep * 0:(time.series.length - 1)
# the time dimension in string format
format.date.seq <- format.Date(date.seq,"%m/%d/%Y")

#> Lambert projected coordinates of the grid cell CENTERS (unit=km).
#> These are the unique x, y coordinates of the grid cell centers -- NOT the coordinates for every grid cell, since the data are on a regular grid.
#> You can also use the get.coord.for.dimension() function to extract the grid cell edges by changing the “position” argument.
ncols <- ncatt_get(cctm.in, varid=0, attname="NCOLS")$value
nrows <- ncatt_get(cctm.in, varid=0, attname="NROWS")$value
x.origin <- ncatt_get(cctm.in, varid=0, attname="XORIG")$value
y.origin <- ncatt_get(cctm.in, varid=0, attname="YORIG")$value
x.cell.width <- ncatt_get(cctm.in, varid=0, attname="XCELL")$value
y.cell.width <- ncatt_get(cctm.in, varid=0, attname="YCELL")$value
#> for the center of each cell
grid.offset <- 0.5
#> x-coordinate locations of grid centers
x.proj.coord <- seq(from = x.origin + grid.offset * x.cell.width, 
                    by = x.cell.width, 
                    length = ncols)
length(x.proj.coord)
#[1] 459
#> y-coordinate locations of grid centers
y.proj.coord <- seq(from = y.origin + grid.offset * y.cell.width, 
                    by = y.cell.width, 
                    length = nrows)
length(y.proj.coord)
#[1] 299

#> Also get the grid cell centers of all of the grid cell with units=meters.  We will use this later when we convert the data to an object raster.
xy.proj.coord.meters <- as.matrix(expand.grid(x.proj.coord, y.proj.coord))
colnames(xy.proj.coord.meters) <- c("x", "y")
dim(xy.proj.coord.meters)
#[1] 137241      2

#> Convert lambert coordinates for grid cell centers to long/lat
p.alpha <-  ncatt_get(cctm.in, varid=0, attname="P_ALP")$value
p.beta <-  ncatt_get(cctm.in, varid=0, attname="P_BET")$value
p.gamma <-  ncatt_get(cctm.in, varid=0, attname="P_GAM")$value
y.center <- ncatt_get(cctm.in, varid=0, attname="YCENT")$value
#> Get the projection string from the file
lambert.proj.string <- paste0("+proj=lcc +lat_1=", p.alpha, 
                              " +lat_2=", p.beta,
                              " +lat_0=", y.center,
                              " +lon_0=", p.gamma,
                              " +a=", EARTH_RADIUS, 
                              " +b=", EARTH_RADIUS)
lambert.proj.string
#[1] "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000"
lonlat.proj.string <- "+proj=longlat"
#> Project from Lambert to longitude/latitude (using terra::project)
long.lat.coord <- as.data.frame(project(xy.proj.coord.meters, 
                                from = lambert.proj.string,
                                to = lonlat.proj.string))

#> US, Canada, Mexico map lines in Lambert projected coordinates.
canusamex.natl <- map("worldHires", regions=c("Canada", "USA", "Mexico"),
                      exact=FALSE, resolution=0, plot=FALSE)
canusamex.natl.lonlat <- cbind(canusamex.natl$x, canusamex.natl$y)
rm(canusamex.natl)

## Get map of state borders, without including the outer boundaries.
## These outer boundaries are provided by the the high-resolution
## national maps created in the previous step.
state <- map("state", exact=F, boundary=FALSE, resolution=0, plot=FALSE)
state.lonlat <- cbind(state$x, state$y)
rm(state)

## Put it national and state boundaries together.
canusamex <- rbind(canusamex.natl.lonlat, matrix(NA, ncol=2), state.lonlat)
rm(canusamex.natl.lonlat, state.lonlat)

#> note warning message below - doesn't seem to hurt anything
map.lines <- project(canusamex, 
                     from = lonlat.proj.string, 
                     to = lambert.proj.string)
#> [project] 651 failed transformations
colnames(map.lines) <- c("x", "y")
rm(canusamex)

#> Extract PM25_AVG (daily 24hr LST averaged PM2.5 concentrations for 1/1/2017 - 12/31/2017).
mod.name <- "PM25_AVG"
mod.unit <- all.mod.units[all.mod.names==mod.name]
mod.array <- ncvar_get(cctm.in,var=mod.name) 
dim (mod.array)
#[1] 459 299 365

#> Create annual average.
mod.annual.avg <- apply(mod.array,c(1,2),mean)
dim(mod.annual.avg)
#[1] 459 299

#> A 'nice' color palette to try for mapping. 
my.color.palette <- colorRampPalette(c("#E6E6E6","#999999","#56B4E9","#0072B2","#009E73","#F0E442","#E69F00","#D55E00","#A52A2A"))
my.range <- c(0,25)
my.color.bins <- seq(0,25,by=1)
n.colors <- length(my.color.bins)-1
my.colors <- my.color.palette(n.colors)
#> Spatial map of annual average PM25 for 2017.
image.plot(x.proj.coord,y.proj.coord,mod.annual.avg,col=my.colors,breaks=my.color.bins,legend.args=list(text=paste(mod.unit)),xlab="",ylab="",axes=F,main=paste("Annual Average",mod.name)) 
#> Add US, Canada, Mexico map lines
lines(map.lines)

#> Create raster object using projection information extracted from the I/O API file. 
xyz <- data.frame(x=xy.proj.coord.meters[,1],y=xy.proj.coord.meters[,2],z=matrix(mod.annual.avg))
#> Create the raster object from the meter coordinates because they 
#>   the grid is regularly spaced
mod.raster <- rast(xyz, crs = lambert.proj.string) # using terra::rast
#> Now project it on latitude/longitude
mod.raster.lonlat <- project(mod.raster, lonlat.proj.string) # using terra::project

#> Now make a lovely leaflet map using the raster object.
my.pal <- colorNumeric(
  palette =  my.color.palette(n.colors),
  domain = my.range
)

mapStates <- map("state",fill=T,plot=F)
my.leaf <- leaflet(data=mapStates) %>% 
  addTiles(options=tileOptions(opacity=1)) %>% 
  addProviderTiles("OpenStreetMap.Mapnik")  %>% 
  addRasterImage(mod.raster,colors=my.pal,opacity=.5) %>% 
  addLegend("bottomright", pal = my.pal, values = my.range,
            title = paste("Annual Average",mod.name),opacity = 1)

my.leaf <- leaflet(data=mapStates) %>% 
  addTiles(options=tileOptions(opacity=1)) %>% 
  addProviderTiles("OpenStreetMap.Mapnik")  %>% 
  addRasterImage(mod.raster.lonlat,colors=my.pal,opacity=.5) %>% 
  addLegend("bottomright", pal = my.pal, values = my.range,
            title = paste("Annual Average",mod.name),opacity = 1)

#> Save the plot as an interactive html file.
saveWidget(my.leaf, file="PM25_leaflet.html",selfcontained=T)