Source

EDGAR-4.2_minus_soil_and_biomass_to_AQMEII-NA / sum_reunit.r

# R code to
# * sum several EDGAR global/unprojected inventories
# * convert sum's units to molar rate

# If running manually in R console, remember to run setup actions: `source ./sum_reunit_regrid_retemp.sh`

# ----------------------------------------------------------------------
# constants
# ----------------------------------------------------------------------

## data

# kludge for my clumsy namespacing
my_this_fn <- Sys.getenv('CALL_REUNIT_FN')
this_fn <- my_this_fn

# all the following env vars must be set and exported in driver script
work_dir <- Sys.getenv('WORK_DIR')
# in_pdf_fp <- Sys.getenv('EDGAR_NONAG_RAW_PDF_FP')
out_pdf_fp <- Sys.getenv('EDGAR_NONAG_REUNIT_PDF_FP')
pdf_er <- Sys.getenv('PDF_VIEWER')
sigdigs <- as.numeric(Sys.getenv('OUTPUT_SIGNIFICANT_DIGITS'))

# in_fp <- Sys.getenv('EDGAR_NONAG_RAW_FP')
raster_rotate <- as.logical( Sys.getenv('ROTATE_INPUT'))
# in_band <- Sys.getenv('EDGAR_NONAG_REGRID_BAND')

# raw/unsummed datavars and their files
combust_fp <- Sys.getenv('EDGAR_NONAG_COMBUST_FP')
nonroad_fp <- Sys.getenv('EDGAR_NONAG_NONROAD_FP')
road_fp <- Sys.getenv('EDGAR_NONAG_ROAD_FP')
resid_fp <- Sys.getenv('EDGAR_NONAG_RESID_FP')
petro_fp <- Sys.getenv('EDGAR_NONAG_PETRO_FP')
ind_fp <- Sys.getenv('EDGAR_NONAG_IND_FP')
manure_fp <- Sys.getenv('EDGAR_NONAG_MANURE_FP')
runoff_fp <- Sys.getenv('EDGAR_NONAG_RUNOFF_FP')
waste_fp <- Sys.getenv('EDGAR_NONAG_WASTE_FP')
fff_fp <- Sys.getenv('EDGAR_NONAG_FFF_FP')
tropo_fp <- Sys.getenv('EDGAR_NONAG_TROPO_FP')
# all have the same metadata (I believe, I hope)
raw_datavar_name <- Sys.getenv('EDGAR_NONAG_RAW_DATAVAR_NAME')
raw_datavar_long_name <- Sys.getenv('EDGAR_NONAG_RAW_DATAVAR_LONG_NAME')
raw_datavar_standard_name <- Sys.getenv('EDGAR_NONAG_RAW_DATAVAR_STANDARD_NAME')
raw_datavar_cell_method <- Sys.getenv('EDGAR_NONAG_RAW_DATAVAR_CELL_METHOD')
raw_datavar_comment <- Sys.getenv('EDGAR_NONAG_RAW_DATAVAR_COMMENT')
raw_datavar_units <- Sys.getenv('EDGAR_NONAG_RAW_DATAVAR_UNITS')
# raw_datavar_na <- as.numeric(Sys.getenv('EDGAR_NONAG_RAW_DATAVAR_NA'))

## gridcell area data
areas_fp <- Sys.getenv('EDGAR_NONAG_AREAS_FP')
areas_datavar_name <- Sys.getenv('EDGAR_NONAG_AREAS_DATAVAR_NAME')

## unit conversion
molar_mass_n2o <- as.numeric(Sys.getenv('MOLAR_MASS_N2O')) # units=g/mol
g.per.kg <- 1000 # grams/kilogram, unitless

## output == summed + 'massified'
out_fp <- Sys.getenv('EDGAR_NONAG_REUNIT_FP')
out_datavar_name <- Sys.getenv('EDGAR_NONAG_REUNIT_DATAVAR_NAME')
out_datavar_type <- Sys.getenv('EDGAR_NONAG_REUNIT_DATAVAR_TYPE')
out_datavar_units <- Sys.getenv('EDGAR_NONAG_REUNIT_DATAVAR_UNITS')
out_datavar_long_name <- Sys.getenv('EDGAR_NONAG_REUNIT_DATAVAR_LONG_NAME')
out_datavar_coord_x_name <- Sys.getenv('EDGAR_NONAG_REUNIT_DATAVAR_COORD_X_NAME')
out_datavar_coord_y_name <- Sys.getenv('EDGAR_NONAG_REUNIT_DATAVAR_COORD_Y_NAME')

# helper
stat_script_fp <- Sys.getenv('STAT_SCRIPT_FP')

## plotting

# plot.script.fp <- './plotLayersForTimestep.r' # Sys.getenv('PLOT_SCRIPT_FP')
# pdf.main.title <- 'N2O emissions' # Sys.getenv('EDGAR_NONAG_REGRID_PDF_TITLE')
# pdf.sub.title <- regrid.datavar.unit
# pdf.title <- sprintf('%s\n%s', pdf.title, regrid.datavar.unit)

global_proj4 <- Sys.getenv('GLOBAL_PROJ4')

# ----------------------------------------------------------------------
# functions
# ----------------------------------------------------------------------

visualize.layer <- function(
  nc.fp,         # path to netCDF datasource ...
  layer,         # ... for data (a RasterLayer)
  datavar.name,  # name of the netCDF data variable
  map.list,      # maps to overlay on data: list of SpatialLines or objects castable thereto 
  pdf.fp,        # path to PDF output
  pdf.height,
  pdf.width
) {

  system(sprintf('ncdump -h %s', nc.fp))
  # Now get stats on data for each actual month in the input: no!
  # until netCDF.stats.to.stdout.by.timestep supports range restriction, just get all timesteps
  netCDF.stats.to.stdout(netcdf.fp=nc.fp, data.var.name=datavar.name)
#  # TODO: make these work as if evaluated @ top level
#  layer          # show it
#  summary(layer) # compare to netCDF.stats above

  # start plot driver
  cat(sprintf(
    '%s: plotting %s (may take awhile! TODO: add progress control)\n',
    'visualize.layer', pdf.fp))
  pdf(file=pdf.fp, width=pdf.width, height=pdf.height)

  library(rasterVis)

  # I want to overlay the data with each map in the list, iteratively.
  # But the following does not work: only the last map in the list shows.
#  plots <- rasterVis::levelplot(layer,
#  #  layers,
#    margin=FALSE,
#  #  names.attr=names(global.df),
#    layout=c(1,length(names(layer))))
#  for (i.map in 1:length(map.list)) {
#    plots <- plots + latticeExtra::layer(
#      # why does this fail if map.shp is local? see 'KLUDGE' in callers :-(
#      sp::sp.lines(
#        as(map.list[[i.map]], "SpatialLines"),
#        lwd=0.8, col='darkgray'))
#  } # end for (i.map in 1:length(map.list))
#  plot(plots)

  # For now, kludge :-( handle lists of length 1 and 2 separately:
  if        (length(map.list) == 1) {

    plot(
      rasterVis::levelplot(layer,
        margin=FALSE,
#        layout=c(1,length(names(layer)))
      ) + latticeExtra::layer(
        sp::sp.lines(
          as(map.list[[1]], "SpatialLines"),
          lwd=0.8, col='darkgray')
      )
    )

  } else if (length(map.list) == 2) {

    plot(
      rasterVis::levelplot(layer,
        margin=FALSE,
#        layout=c(1,length(names(layer)))
      ) + latticeExtra::layer(
        sp::sp.lines(
          as(map.list[[1]], "SpatialLines"),
          lwd=0.8, col='darkgray')
      ) + latticeExtra::layer(
        sp::sp.lines(
          as(map.list[[2]], "SpatialLines"),
          lwd=0.8, col='darkgray')
      )
    )

  } else {
    stop(sprintf('visualize.layer: ERROR: cannot handle (length(map.list)==%i', length(map.list)))
  }
  
  dev.off()

} # end visualize.layer <- function

# ----------------------------------------------------------------------
# code
# ----------------------------------------------------------------------

# ----------------------------------------------------------------------
# setup
# ----------------------------------------------------------------------

# accelerate R graphics over SSH, per Adam Wilson
# http://planetflux.adamwilson.us/2012/03/r-graphics-via-ssh.html
X11.options(type="Xlib")

source(stat_script_fp) # in script, produces errant error=
#> netCDF.stats.to.stdout.r: no arguments supplied, exiting

# get a global map
library(maps)
map.world.unproj <- maps::map('world', plot=FALSE)

# get a North American map
library(maptools)
data(wrld_simpl) # from maptools

# coordinate reference system:
global.crs <- sp::CRS(global_proj4)

# "create" world map
map.world.unproj.shp <-
  maptools::map2SpatialLines(map.world.unproj, proj4string=global.crs)
# summary(map.world.unproj.shp) # debugging

# ----------------------------------------------------------------------
# payload
# ----------------------------------------------------------------------

library(raster)

### load inputs as separate layers, then stack them
combust.raster <- raster::raster(combust_fp, varname=raw_datavar_name)
nonroad.raster <- raster::raster(nonroad_fp, varname=raw_datavar_name)
road.raster <- raster::raster(road_fp, varname=raw_datavar_name)
resid.raster <- raster::raster(resid_fp, varname=raw_datavar_name)
petro.raster <- raster::raster(petro_fp, varname=raw_datavar_name)
ind.raster <- raster::raster(ind_fp, varname=raw_datavar_name)
manure.raster <- raster::raster(manure_fp, varname=raw_datavar_name)
runoff.raster <- raster::raster(runoff_fp, varname=raw_datavar_name)
waste.raster <- raster::raster(waste_fp, varname=raw_datavar_name)
fff.raster <- raster::raster(fff_fp, varname=raw_datavar_name)
tropo.raster <- raster::raster(tropo_fp, varname=raw_datavar_name)

cat(sprintf('%s: stacking raw inputs\n', this_fn)) # debugging
in.stack <- raster::stack( # will fail if !all have same grid definitions
  combust.raster, nonroad.raster, road.raster, resid.raster,
  petro.raster, ind.raster, manure.raster, runoff.raster,
  waste.raster, fff.raster, tropo.raster)

### sum raw inputs
cat(sprintf(  # debugging
  '%s: summing fluxes of stacked raw inputs (may take awhile)\n',this_fn))
# sum_raster <- raster::sum(in.stack) # 'sum' is not an exported object from 'namespace:raster'
sum.raster <- sum(in.stack)

# start debugging-----------------------------------------------------

sum.raster
# class       : RasterLayer 
# dimensions  : 1800, 3600, 6480000  (nrow, ncol, ncell)
# resolution  : 0.1, 0.1  (x, y)
# extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 
# data source : in memory
# names       : layer 
# values      : 0, 1.041793e-08  (min, max)

#   end debugging-----------------------------------------------------

### "massify" raw inputs: convert flux-rate values to molar-mass-rate values
cat(sprintf( # debugging
  '%s: converting summed raster from flux rate to molar-mass rate (may take awhile)\n',
  this_fn))

## flux rate (kg/m^2/s) * gridcell area -> mass rate (kg/s)
# fortunately, gridcell areas are provided by EDGAR
areas.raster <- raster::raster(areas_fp, varname=areas_datavar_name)
# re `prod`:
# * does cell-by-cell matrix multiply, not linear-algebra-style
# * also generic
mass.raster <- prod(sum.raster, areas.raster)

## mass rate (kg/s) -> molar rate (mol/s)
# (kg/s) * (g/kg) * (mol/g) -> (mol/s)
out.raster <- (mass.raster * g.per.kg) / molar_mass_n2o

# start debugging-----------------------------------------------------

out.raster


#   end debugging-----------------------------------------------------

### 0-360 longitudes -> -180-+180
if (raster_rotate) {
  out.raster <- rotate(
    out.raster,      # )
    overwrite=TRUE) # else levelplot does one layer per page?
}

### write sum/massified output
cat(sprintf('%s: writing summed-and-massed raster to path=%s ...\n', # debugging
  this_fn, out_fp))
writeRaster(
  x=out.raster, filename=out_fp, varname=out_datavar_name,
  format="CDF", overwrite=TRUE,
# let raster choose default, since "there shouldn't be any"
#  NAflag=out_datavar_na,
  # datavar type will truncate unless set?
  datatype=out_datavar_type,
  # netCDF-specific arguments
  varunit=out_datavar_units, longname=out_datavar_long_name,
  xname=out_datavar_coord_x_name, yname=out_datavar_coord_y_name
)

### visualize output

# Why does '=' fail and '<-' succeed in the arg list?
visualize.layer(
  nc.fp=out_fp,
  layer=out.raster,
  datavar.name=out_datavar_name,
#  map.list=list(map.world.unproj.shp),
  map.list <- list(map.world.unproj.shp),
  pdf.fp=out_pdf_fp,
  pdf.height=5,
  pdf.width=5
)