Commits

Tom Roche committed a87e6d7

initial commit of base function + visualization

* regrids from global/unprojected to AQMEII-NA/LCC
* all plotting via rasterVis::levelplot
* all parameters passed by environment variables

TODO:
* add Canada and Mexico to maps (currently just CONUS states)
* remove dummy months={1,14}
* "retemporalize" from monthly to hourly

Tested on terrae. (Requires ncdf4, not yet available on tlrPanP5.)

Comments (0)

Files changed (5)

+*~
+*.nc      # output netCDF
+*.pdf     # graphical output
+Uses [R][R @ wikipedia] code to
+
+1.  *(done)* 2D-regrid [CLM-CN][] data (in [netCDF][] format) from global/unprojected to a projected subdomain ([AQMEII-NA][]).
+2. *(planned)* "retemporalize" netCDF from monthly (original) to hourly
+
+Currently does not provide a clean or general-purpose solution! but merely shows how to do these tasks using R and packages including
+
+* [ncdf4][]
+* [raster][]
+* [rasterVis][]
+
+[R @ wikipedia]: http://en.wikipedia.org/wiki/R_%28programming_language%29
+[NCL @ wikipedia]: http://en.wikipedia.org/wiki/NCAR_Command_Language
+[netCDF]: http://en.wikipedia.org/wiki/NetCDF#Format_description
+[CLM-CN]: https://github.com/TomRoche/cornbeltN2O/wiki/Simulation-of-N2O-Production-and-Transport-in-the-US-Cornbelt-Compared-to-Tower-Measurements#wiki-CLM-CN-N2O
+[AQMEII-NA]: https://github.com/TomRoche/cornbeltN2O/wiki/AQMEII-North-American-domain
+[ncdf4]: http://cran.r-project.org/web/packages/ncdf4/
+[raster]: http://cran.r-project.org/web/packages/raster/
+[rasterVis]: http://cran.r-project.org/web/packages/rasterVis/
+
+To run this code:
+
+1. `git clone` this repo.
+1. `cd` to its working directory (where you cloned it to).
+1. Open the driver (bash) script `regrid_global_to_AQMEII.sh` in an editor! You will probably need to edit it to make it work on your platform. Notably you will probably want to point it to your R executable and PDF viewer.
+2. Run the driver:
+    `$ ./regrid_global_to_AQMEII.sh`
+3. This will download input, then run
+    * an R script to display the input netCDF data, regrid it, and plot the output. The driver should display [an input PDF][global plots] and [an output PDF][AQMEII plots] if properly configured.
+
+[global plots]: https://bitbucket.org/tlroche/clm_cn_global_to_aqmeii-na/downloads/2008PTONCLMCNN2O.pdf
+[AQMEII plots]: https://bitbucket.org/tlroche/clm_cn_global_to_aqmeii-na/downloads/2008PTONCLMCNN2O_regrid.pdf

netCDF.stats.to.stdout.r

+# R code to write simple stats (min, mean, median, max) for an input file.
+# Run this like
+# $ Rscript ./netCDF.stats.to.stdout.r netcdf.fp=./GEIA_N2O_oceanic.nc var.name=emi_n2o
+# or
+# > source('./netCDF.stats.to.stdout.r')
+# > netCDF.stats.to.stdout(...)
+
+# ----------------------------------------------------------------------
+# constants
+# ----------------------------------------------------------------------
+
+this.fn <- 'netCDF.stats.to.stdout.r'      # TODO: get like $0
+
+# ----------------------------------------------------------------------
+# functions
+# ----------------------------------------------------------------------
+
+# syntactic sugar
+q1 <- function(vec) { quantile(vec, 0.25) } # first quartile
+q3 <- function(vec) { quantile(vec, 0.75) } # third quartile
+
+# the main event
+netCDF.stats.to.stdout <- function(
+  netcdf.fp, # /path/to/netcdf/file, can be relative or FQ
+  data.var.name,   # name of data variable (datavar) of interest
+  stats.precision=3 # sigdigs to use for min, median, max of obs  
+) {
+
+  # TODO: test arguments!
+
+#  # start debug
+#  cat(sprintf(
+#    '%s: netcdf.fp==%s, data.var.name==%s\n',
+#    this.fn, netcdf.fp, data.var.name))
+#  system(sprintf('ls -alth %s', netcdf.fp))
+#  system(sprintf('ncdump -h %s', netcdf.fp))
+#  #   end debug
+
+  # create strings for use in output below (far below)
+  # double-sprintf-ing to set precision by constant: cool or brittle?
+  stat.str.template <- sprintf('%%.%ig', stats.precision)
+  # use these in function=subtitle.stats as sprintf inputs
+  max.str.template <- sprintf('max=%s', stat.str.template)
+  mea.str.template <- sprintf('mean=%s', stat.str.template)
+  med.str.template <- sprintf('med=%s', stat.str.template)  # median==2nd quartile
+  min.str.template <- sprintf('min=%s', stat.str.template)
+  q1.str.template <- sprintf('q1=%s', stat.str.template)    # 1st quartile
+  q3.str.template <- sprintf('q3=%s', stat.str.template)    # 3rd quartile
+
+  # needed to parse netCDF
+  library(ncdf4)
+
+  # open netCDF file, uncautiously
+  # NOTE: you must assign when you nc_open!
+  netcdf.file <- ncdf4::nc_open(
+    filename=netcdf.fp,
+    write=FALSE,    # will only read below
+    readunlim=TRUE) # it's a small file
+
+  # uncautiously get the data out of the datavar
+  data.var.data <- ncvar_get(
+    nc=netcdf.file,
+    varid=data.var.name)
+
+#  # start debug
+#  cat(sprintf('%s: data.var.name==%s has size==\n',
+#    this.fn, data.var.name))
+#  print(dim(data.var.data))
+#  #   end debug
+
+  if (is.numeric(data.var.data) && sum(!is.na(data.var.data))) {
+#    unsparse.data <- data.var.data[!is.na(data.var.data)]
+    # collapse its structure
+    unsparse.data <- c(data.var.data[!is.na(data.var.data)])
+    obs <- length(unsparse.data)
+    if (obs > 0) {
+      cells.str <- sprintf('cells=%i', length(data.var.data))
+      obs.str <- sprintf('obs=%i', obs)
+      min.str <- sprintf(min.str.template, min(unsparse.data))
+      q1.str <- sprintf(q1.str.template, q1(unsparse.data))
+      mea.str <- sprintf(mea.str.template, mean(unsparse.data))
+      med.str <- sprintf(med.str.template, median(unsparse.data))
+      q3.str <- sprintf(q3.str.template, q3(unsparse.data))
+      max.str <- sprintf(max.str.template, max(unsparse.data))
+
+      # at last: output!
+      cat(sprintf('For %s var=%s\n', netcdf.fp, data.var.name))
+      cat(sprintf('\t%s\n', cells.str))
+      cat(sprintf('\t%s\n', obs.str))
+      # 6-number summary
+      cat(sprintf('\t%s\n', min.str))
+      cat(sprintf('\t%s\n', q1.str))
+      cat(sprintf('\t%s\n', med.str))
+      cat(sprintf('\t%s\n', mea.str))
+      cat(sprintf('\t%s\n', q3.str))
+      cat(sprintf('\t%s\n', max.str))
+
+    } else {
+      cat(sprintf('%s: %s var=%s has no non-NA data',
+        this.fn, data.var.name, netcdf.fp))
+    }
+  } else {
+    cat(sprintf('%s: %s var=%s has no numeric non-NA data',
+      this.fn, data.var.name, netcdf.fp))
+  }
+
+  # teardown
+  nc_close(netcdf.file)
+
+} # end function netCDF.stats.to.stdout
+
+netCDF.stats.to.stdout.by.timestep <- function(
+  netcdf.fp, # /path/to/netcdf/file, can be relative or FQ
+  data.var.name,   # name of data variable (datavar) of interest
+  time.dim.name='time', # of time dimension in datavar args
+  stats.precision=3 # sigdigs to use for min, median, max of obs
+) {
+
+  # TODO: test arguments!
+
+#  # start debug
+#  cat(sprintf(
+#    '%s: netcdf.fp==%s, data.var.name==%s, time.dim.name=%s\n',
+#    this.fn, netcdf.fp, data.var.name, time.dim.name))
+#  system(sprintf('ls -alth %s', netcdf.fp))
+#  system(sprintf('ncdump -h %s', netcdf.fp))
+#  #   end debug
+
+  # create strings for use in output below (far below)
+  # double-sprintf-ing to set precision by constant: cool or brittle?
+  stat.str.template <- sprintf('%%.%ig', stats.precision)
+  # use these in function=subtitle.stats as sprintf inputs
+  max.str.template <- sprintf('max=%s', stat.str.template)
+  mea.str.template <- sprintf('mean=%s', stat.str.template)
+  med.str.template <- sprintf('med=%s', stat.str.template)  # median==2nd quartile
+  min.str.template <- sprintf('min=%s', stat.str.template)
+  q1.str.template <- sprintf('q1=%s', stat.str.template)    # 1st quartile
+  q3.str.template <- sprintf('q3=%s', stat.str.template)    # 3rd quartile
+
+  # needed to parse netCDF
+  library(ncdf4)
+
+  # open netCDF file, uncautiously
+  # NOTE: you must assign when you nc_open!
+  netcdf.file <- ncdf4::nc_open(
+    filename=netcdf.fp,
+    write=FALSE,    # will only read below
+    readunlim=TRUE) # it's a small file? TODO: read by timestep!
+
+  # Find the index of the datavar of interest
+  # (in the list of all non-coordinate variables).
+  netcdf.file.vars.n <- netcdf.file$nvars
+  data.var.index <- -1 # an invalid index
+  for (i in 1:netcdf.file.vars.n) {
+    netcdf.file.var <- netcdf.file$var[[i]]
+    if (netcdf.file.var$name == data.var.name) {
+      data.var.index <- i
+    }
+  }
+  if (data.var.index == -1) {
+    stop(sprintf('%s: ERROR: failed to find data.var.index\n', this.fn))
+  }
+
+#  # start debug
+#  cat(sprintf('%s: data.var.name==%s has var posn=%i and size==\n',
+#    this.fn, data.var.name, data.var.index))
+#  print(dim(data.var.data))
+#  #   end debug
+
+  # can only compute n.timesteps once we have the datavar's dims ...
+  data.var <- netcdf.file$var[[data.var.index]]
+  data.var.dims <- data.var$size
+  data.var.dims.n <- data.var$ndims
+  # ... and particularly the position of dimension$name==time.dim.name
+  time.dim.index <- -1 # an invalid index
+  for (i in 1:data.var.dims.n) {
+    data.var.dim <- data.var$dim[[i]]
+    if (data.var.dim$name == time.dim.name) {
+      time.dim.index <- i
+    }
+  }
+  if (time.dim.index == -1) {
+    stop(sprintf('%s: ERROR: failed to find time.dim.index\n', this.fn))
+  } # else
+  n.timesteps <- data.var.dims[time.dim.index]
+
+#  # start debug
+#  cat(sprintf('%s: time.dim.name==%s has dim posn=%i and size=%i\n',
+#    this.fn, time.dim.name, time.dim.index, n.timesteps))
+#  #   end debug
+
+# TODO: ncdf4 bug?
+#  # compute read vectors for use in `ncvar_get`: read all timestep values
+#  vec.start.template <- rep(1, data.var.dims.n)   # ASSERT: constant
+##  vec.count.template <- rep(-1, data.var.dims.n)  # ASSERT: constant
+#  # see bug below, try instead
+#  vec.count.template <- data.var.dims             # ASSERT: constant
+
+  # iterate over the timesteps
+  for (i in 1:n.timesteps) {
+
+#    vec.start <- vec.start.template
+#    vec.start[time.dim.index] <- i # only get the i'th timestep
+#    vec.count <- vec.count.template
+#    vec.count[time.dim.index] <- i # only get the i'th timestep
+#
+#    # start debug
+#    cat(sprintf('%s: vec.start==\n', this.fn))
+#    print(vec.start)
+#    cat(sprintf('%s: vec.count==\n', this.fn))
+#    print(vec.count)
+#    #   end debug
+
+    data.var.data <- ncdf4::ncvar_get(
+      nc=netcdf.file,
+      varid=data.var.name
+# TODO: {debug, bug report} why this fails
+#      varid=data.var.name,
+#      start=vec.start,
+#      count=vec.count
+    )
+# each loop call increases time dimension! e.g.,
+# > netCDF.stats.to.stdout.r: dim(data.var.data)==
+# > [1] 144  96
+# ...
+# > netCDF.stats.to.stdout.r: dim(data.var.data)==
+# > [1] 144  96   2
+# ...
+# > netCDF.stats.to.stdout.r: dim(data.var.data)==
+# > [1] 144  96   3
+    # workaround: read unlimited (above), then just take the slice of interest
+    data.var.timestep <- data.var.data[,,i]
+                                  
+#    # start debug
+#    cat(sprintf('%s: dim(data.var.timestep)==\n', this.fn))
+#    print(dim(data.var.timestep))
+#    cat(sprintf('%s: summary(data.var.timestep)==\n', this.fn))
+#    print(summary(c(data.var.timestep))) # collapse its structure
+#    #   end debug
+
+    if (is.numeric(data.var.timestep) && sum(!is.na(data.var.timestep))) {
+#      unsparse.data <- data.var.timestep[!is.na(data.var.timestep)]
+      # collapse its structure
+      unsparse.data <- c(data.var.timestep[!is.na(data.var.timestep)])
+      obs <- length(unsparse.data)
+      if (obs > 0) {
+        cells.str <- sprintf('cells=%i', length(data.var.timestep))
+        obs.str <- sprintf('obs=%i', obs)
+        min.str <- sprintf(min.str.template, min(unsparse.data))
+        q1.str <- sprintf(q1.str.template, q1(unsparse.data))
+        mea.str <- sprintf(mea.str.template, mean(unsparse.data))
+        med.str <- sprintf(med.str.template, median(unsparse.data))
+        q3.str <- sprintf(q3.str.template, q3(unsparse.data))
+        max.str <- sprintf(max.str.template, max(unsparse.data))
+
+        # at last: output!
+        cat(sprintf('For %s data.var=%s, timestep=%i of %i\n',
+          netcdf.fp, data.var.name, i, n.timesteps))
+        cat(sprintf('\t%s\n', cells.str))
+        cat(sprintf('\t%s\n', obs.str))
+        # 6-number summary
+        cat(sprintf('\t%s\n', min.str))
+        cat(sprintf('\t%s\n', q1.str))
+        cat(sprintf('\t%s\n', med.str))
+        cat(sprintf('\t%s\n', mea.str))
+        cat(sprintf('\t%s\n', q3.str))
+        cat(sprintf('\t%s\n', max.str))
+
+#        rm( # debugging why output is same for each timestep
+#          cells.str,
+#          obs.str,
+#          min.str,
+#          q1.str,
+#          med.str,
+#          mea.str,
+#          q3.str,
+#          max.str)
+
+      } else {
+        cat(sprintf('%s: %s data.var=%s has no non-NA data',
+          this.fn, data.var.name, netcdf.fp))
+      }
+    } else {
+      cat(sprintf('%s: %s data.var=%s has no numeric non-NA data',
+        this.fn, data.var.name, netcdf.fp))
+    }
+
+    rm( # debugging why output is same for each timestep
+#      data.var.data,
+      data.var.timestep,
+      unsparse.data)
+
+  } # end for (i in 1:n.timesteps)
+
+  # teardown
+  ncdf4::nc_close(netcdf.file)
+  rm(this.fn) # hopefully prevents overwriting `this.fn` in prior namespaces
+} # end function netCDF.stats.to.stdout.by.timestep
+
+# ----------------------------------------------------------------------
+# code
+# ----------------------------------------------------------------------
+
+# if this is called as a script, provide a main(): see
+# https://stat.ethz.ch/pipermail/r-help/2012-September/323551.html
+# https://stat.ethz.ch/pipermail/r-help/2012-September/323559.html
+if (!interactive()) {
+
+# start debug
+#  cat(sprintf('%s: interactive()==TRUE\n', this.fn))
+#   end debug
+  
+  # TODO: fix `strsplit` regexp below to make this unnecessary
+  library(stringr)
+
+  # pass named arguments: var above separated by '='
+
+  args <- commandArgs(TRUE)
+  # args is now a list of character vectors
+  # First check to see if any arguments were passed, then evaluate each argument:
+  # assign val (RHS) to key (LHS) for arguments of the (required) form 'key=val'
+  if (length(args)==0) {
+    cat(sprintf('%s: no arguments supplied, exiting\n', this.fn))
+#    q(status=1) # KLUDGE:
+# Currently this is not seeing arguments when called from Rscript,
+# so this exit also exits the caller :-(    
+  } else {
+  # simple positional args work
+  # TODO: also support positional usage
+  #  netcdf.fp <- args[1]
+  #  var.name <- args[2]
+    # TODO: test arg length: 2 is required!
+
+# start debug
+#    cat(sprintf('%s: got length(args)==%i\n', this.fn, length(args)))
+#   end debug
+
+    for (i in 1:length(args)) {
+#       eval(parse(text=args[[i]]))
+      # `eval(parse())` is unsafe and requires awkward quoting:
+      # e.g., of the following (bash) commandlines
+
+      # - Rscript ./netCDF.stats.to.stdout.r netcdf.fp="GEIA_N2O_oceanic.nc" var.name="emi_n2o"
+      #   fails
+
+      # + Rscript ./netCDF.stats.to.stdout.r 'netcdf.fp="GEIA_N2O_oceanic.nc"' 'var.name="emi_n2o"'
+      #   succeeds
+
+      # so instead
+      # TODO: use package `optparse` or `getopt`
+      args.keyval.list <-
+        strsplit(as.character(parse(text=args[[i]])),
+          split='[[:blank:]]*<-|=[[:blank:]]*', fixed=FALSE)
+  #                            split='[ \t]*<-|=[ \t]*', fixed=FALSE)
+      args.keyval.vec <- unlist(args.keyval.list, recursive=FALSE, use.names=FALSE)
+      # TODO: test vector elements!
+      # Neither wants to remove all whitespace from around arguments :-( so
+      args.key <- str_trim(args.keyval.vec[1], side="both")
+      args.val <- str_trim(args.keyval.vec[2], side="both")
+
+# start debug
+#       cat(sprintf('%s: got\n', this.fn))
+#       cat('\targs.keyval.list==\n')
+#       print(args.keyval.list)
+#       cat('\targs.keyval.vec==\n')
+#       print(args.keyval.vec)
+#       cat(sprintf('\targs.key==%s\n', args.key))
+#       cat(sprintf('\targs.val==%s\n', args.val))
+#   end debug
+
+      # A real case statement would be nice to have
+      if        (args.key == 'netcdf.fp') {
+        netcdf.fp <- args.val
+      } else if (args.key == 'var.name') {
+        var.name <- args.val
+      } else {
+        stop(sprintf("unknown argument='%s'", args.key))
+        # TODO: show usage
+        q(status=1) # exit with error
+      }
+    } # end for loop over arguments
+
+    # payload!
+    netCDF.stats.to.stdout(netcdf.fp, var.name)
+
+  } # end if testing number of arguments
+} # end if (!interactive())

regrid_global_to_AQMEII.r

+# R code to 2D-regrid CLM-CN global/unprojected netCDF data to AQMEII-NA, an LCC-projected subdomain
+
+# If running manually in R console, remember to run setup actions: `source ./regrid_global_to_AQMEII.sh`
+
+# ----------------------------------------------------------------------
+# constants
+# ----------------------------------------------------------------------
+
+## data
+
+# kludge for my clumsy namespacing
+my.this.fn <- Sys.getenv('CALL_REGRID_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('CLM_CN_RAW_PDF_FP')
+out.pdf.fp <- Sys.getenv('CLM_CN_REGRID_PDF_FP')
+pdf.er <- Sys.getenv('PDF_VIEWER')
+sigdigs <- as.numeric(Sys.getenv('OUTPUT_SIGNIFICANT_DIGITS'))
+
+in.fp <- Sys.getenv('CLM_CN_RAW_FP')
+raster.rotate <- as.logical( Sys.getenv('ROTATE_INPUT'))
+# in.band <- Sys.getenv('CLM_CN_REGRID_BAND')
+
+raw.datavar.name <- Sys.getenv('CLM_CN_RAW_DATAVAR_NAME')
+raw.datavar.longname <- Sys.getenv('CLM_CN_RAW_DATAVAR_LONGNAME')
+raw.datavar.unit <- Sys.getenv('CLM_CN_RAW_DATAVAR_UNIT')
+# raw.datavar.na <- as.numeric(Sys.getenv('CLM_CN_RAW_DATAVAR_NA'))
+time.dim.name <- Sys.getenv('CLM_CN_RAW_TIME_VAR_NAME')
+
+template.var.name <- Sys.getenv('TEMPLATE_VAR_NAME')
+template.in.fp <- Sys.getenv('TEMPLATE_INPUT_FP')
+
+regrid.datavar.name <- raw.datavar.name
+regrid.datavar.longname <- raw.datavar.longname
+regrid.datavar.unit <- raw.datavar.unit
+# regrid.datavar.na <- raw.datavar.na
+
+out.fp <- Sys.getenv('CLM_CN_REGRID_FP')
+x.var.name <- Sys.getenv('CLM_CN_REGRID_X_VAR_NAME')
+y.var.name <- Sys.getenv('CLM_CN_REGRID_Y_VAR_NAME')
+z.var.name <- Sys.getenv('CLM_CN_REGRID_TIME_VAR_NAME')
+z.var.unit <- Sys.getenv('CLM_CN_REGRID_TIME_VAR_UNIT')
+
+stat.script.fp <- './netCDF.stats.to.stdout.r' # Sys.getenv('STAT_SCRIPT_FP')
+
+## plotting
+
+# plot.script.fp <- './plotLayersForTimestep.r' # Sys.getenv('PLOT_SCRIPT_FP')
+# pdf.main.title <- 'N2O emissions' # Sys.getenv('CLM_CN_REGRID_PDF_TITLE')
+# pdf.sub.title <- regrid.datavar.unit
+# pdf.title <- sprintf('%s\n%s', pdf.title, regrid.datavar.unit)
+
+# get a global map
+library(maps)
+map.world.unproj <- maps::map('world', plot=FALSE)
+global.proj4 <- '+proj=longlat +ellps=WGS84'
+
+# get a North American map
+library(maptools)
+data(wrld_simpl) # from maptools
+
+# # package=M3 map for fields::image.plot
+# map.table.fp <- './map.CMAQm.world.dat' # Sys.getenv('MAP_TABLE_FP')
+# map.cmaq <- read.table(map.table.fp, sep=",")
+# palette.vec <- c(
+# # original from KMF, 3 colors added to get deciles in probabilities.vec
+#   #              R color
+#   #                code
+#   "grey",         # 260
+#   "purple",       # 547
+#   "deepskyblue2", # 123  
+#   "green",        # 254
+#   "greenyellow",  # 259
+#   "yellow",       # 652
+#   "orange",       # 498
+#   "orangered",    # 503
+#   "red",          # 552
+#   "red4",         # 556
+#   "brown"         #  32
+# )
+# colors <- colorRampPalette(palette.vec)
+# # used for quantiling legend
+# probabilities.vec <- seq(0, 1, 1.0/(length(palette.vec) - 1))
+# # probabilities.vec
+# # [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
+
+# ----------------------------------------------------------------------
+# functions
+# ----------------------------------------------------------------------
+
+visualize.layers <- function(
+  nc.fp,           # path to netCDF datasource ...
+  brick,           # ... for data (a RasterBrick)
+  datavar.name,    # name of the netCDF data variable
+  layer.dim.name,  # name of the datavar dimension indexing the layers (e.g., 'time')
+# TODO: take vector of maps
+  map.shp,         # SpatialLines containing map to overlay on data
+  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.by.timestep(
+    netcdf.fp=nc.fp, data.var.name=datavar.name, time.dim.name=layer.dim.name)
+#  # TODO: make these work as if evaluated @ top level
+#  brick          # show it
+#  summary(brick) # compare to netCDF.stats above
+
+  # start plot driver
+  cat(sprintf(
+    '%s: plotting %s (may take awhile! TODO: add progress control)\n',
+    'visualize.layers', pdf.fp))
+  pdf(file=pdf.fp, width=pdf.width, height=pdf.height)
+
+  library(rasterVis)
+#  rasterVis::levelplot(brick,
+  # since we're inside a function
+  plot(rasterVis::levelplot(brick,
+  #  layers,
+    margin=FALSE,
+  #  names.attr=names(global.df),
+    layout=c(1,length(names(brick)))
+  ) +
+    latticeExtra::layer(
+# why does this fail if map.shp is local? see 'KLUDGE' in callers :-(
+      sp::sp.lines(map.shp, lwd=0.8, col='darkgray')))
+# ... and this works?
+#      sp::sp.lines(map.world.unproj.shp, lwd=0.8, col='darkgray')))
+  dev.off()
+
+} # end visualize.layers <- function
+
+# ----------------------------------------------------------------------
+# payload
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# 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
+# source(plot.script.fp)
+
+# coordinate reference system:
+# use package=M3 to get CRS from template file
+library(M3)
+out.proj4 <- M3::get.proj.info.M3(template.in.fp)
+# cat(sprintf('out.proj4=%s\n', out.proj4)) # debugging
+# out.proj4=+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000
+out.crs <- sp::CRS(out.proj4)
+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
+
+# "create" unprojected North American map
+map.noram.shp.unproj <- wrld_simpl[wrld_simpl$ISO3 %in% c('CAN', 'MEX', 'USA'),]
+# > class(map.noram.shp.unproj)
+# [1] "SpatialPolygonsDataFrame"
+# attr(,"package")
+# [1] "sp"
+# > slotNames(map.noram.shp.unproj)
+# [1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"
+# > class(map.noram.shp.unproj@data)
+# [1] "data.frame"
+## > class(map.noram.shp.unproj@polygons)
+## [1] "list"
+# > head(map.noram.shp.unproj@polygons)
+# DON'T DO THIS--it's huge
+
+# ----------------------------------------------------------------------
+# (minimally) process input
+# ----------------------------------------------------------------------
+
+library(raster)
+# in.raster <- raster::raster(in.fp, varname=raw.datavar.name, band=in.band)
+in.raster <- raster::brick(in.fp, varname=raw.datavar.name)
+# replace the layer names
+in.raster.layers.n <- length(names(in.raster))
+names(in.raster) <- as.character(c(1:in.raster.layers.n))
+# correct zero-based longitudes. TODO: test first!
+if (raster.rotate) {
+  in.raster <- rotate(
+    in.raster,      # )
+    overwrite=TRUE) # else levelplot does one layer per page?
+}
+
+# # start debugging
+# in.raster
+# summary(in.raster) # compare to netCDF.stats following
+# #   end debugging
+
+# > in.raster
+# class       : RasterBrick 
+# dimensions  : 96, 144, 13824, 14  (nrow, ncol, ncell, nlayers)
+# resolution  : 2.5, 1.875  (x, y)
+# extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
+# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
+# data source : in memory
+# names       :        X1,        X2,        X3,        X4,        X5,        X6,        X7,        X8,        X9,       X10,       X11,       X12,       X13,       X14 
+# min values  :         0,         0,         0,         0,         0,         0,         0,         0,         0,         0,         0,         0,         0,         0 
+# max values  :  451.1202,  451.1202,  526.0969,  479.2770,  523.1622,  662.8298, 1183.8647, 1011.5519, 1961.3065,  583.5884,  590.4103,  507.7597,  904.2662,  904.2662 
+# z-value     : 732892, 732923, 732954, 732983, 733014, 733044, 733075, 733105, 733136, 733167, 733197, 733228, 733258, 733289 
+
+# > summary(in.raster) # zeros reformatted for clarity
+#                 [,1]         [,2]         [,3]         [,4]         [,5]
+# Min.    0            0            0            0            0           
+# 1st Qu. 0            0            0            0            0           
+# Median  0            0            0            0            0           
+# 3rd Qu. 1.122448e-08 1.122448e-08 3.469217e-09 4.834096e-08 1.429877e-05
+# Max.    4.511202e+02 4.511202e+02 5.260969e+02 4.792770e+02 5.231622e+02
+# NA's    0            0            0            0            0           
+#                 [,6]         [,7]         [,8]         [,9]        [,10]
+# Min.    0            0            0            0            0           
+# 1st Qu. 0            0            0            0            0           
+# Median  0            0            0            0            0           
+# 3rd Qu. 1.534622e-04 1.763448e-03 5.498378e-03 5.348118e-03 2.348898e-03
+# Max.    6.628298e+02 1.183865e+03 1.011552e+03 1.961307e+03 5.835884e+02
+# NA's    0            0            0            0            0           
+#                [,11]        [,12]        [,13]        [,14]
+# Min.    0            0            0            0           
+# 1st Qu. 0            0            0            0           
+# Median  0            0            0            0           
+# 3rd Qu. 7.738026e-04 3.494962e-07 4.448062e-08 4.448062e-08
+# Max.    5.904103e+02 5.077597e+02 9.042662e+02 9.042662e+02
+# NA's    0            0            0            0           
+
+# ----------------------------------------------------------------------
+# visualize input
+# ----------------------------------------------------------------------
+
+# KLUDGE!
+map.shp <- map.world.unproj.shp
+
+visualize.layers(
+  nc.fp=in.fp,
+  brick=in.raster,
+  datavar.name=raw.datavar.name,
+  layer.dim.name=time.dim.name,
+  map.shp=map.world.unproj.shp,
+  pdf.fp=in.pdf.fp,
+  pdf.height=25,
+  pdf.width=5
+)
+
+# globe displays normally, dim=lon-lat
+# interesting Asian spikes 6-9
+
+# ----------------------------------------------------------------------
+# regrid monthly
+# ----------------------------------------------------------------------
+
+# use M3 to get extents from template file (thanks CGN!)
+extents.info <- M3::get.grid.info.M3(template.in.fp)
+extents.xmin <- extents.info$x.orig
+extents.xmax <- max(
+  M3::get.coord.for.dimension(
+    file=template.in.fp, dimension="col", position="upper", units="m")$coords)
+extents.ymin <- extents.info$y.orig
+extents.ymax <- max(
+  M3::get.coord.for.dimension(
+    file=template.in.fp, dimension="row", position="upper", units="m")$coords)
+grid.res <- c(extents.info$x.cell.width, extents.info$y.cell.width) # units=m
+
+template.extents <-
+  raster::extent(extents.xmin, extents.xmax, extents.ymin, extents.ymax)
+# template.extents # debugging
+# class       : Extent 
+# xmin        : -2556000 
+# xmax        : 2952000 
+# ymin        : -1728000 
+# ymax        : 1860000 
+
+template.in.raster <-
+#  raster::raster(template.in.fp, varname=template.var.name, band=template.band)
+  raster::raster(template.in.fp, varname=template.var.name)
+template.raster <- raster::projectExtent(template.in.raster, crs=out.crs)
+#> Warning message:
+#> In projectExtent(template.in.raster, out.proj4) :
+#>   158 projected point(s) not finite
+# is that "projected point(s) not finite" warning important? Probably not, per Hijmans
+
+# without this, extents aren't correct!
+template.raster@extent <- template.extents
+# should resemble the domain specification @
+# https://github.com/TomRoche/cornbeltN2O/wiki/AQMEII-North-American-domain#wiki-EPA
+
+# template.raster # debugging
+# class       : RasterLayer 
+# dimensions  : 299, 459, 137241  (nrow, ncol, ncell)
+# resolution  : 12000, 12000  (x, y)
+# extent      : -2556000, 2952000, -1728000, 1860000  (xmin, xmax, ymin, ymax)
+# coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000 
+
+# at last: do the regridding
+out.raster <-
+  raster::projectRaster(
+    # give a template with extents--fast, but gotta calculate extents
+    from=in.raster, to=template.raster, crs=out.crs,
+    # give a resolution instead of a template? no, that hangs
+#    from=in.raster, res=grid.res, crs=out.proj4,
+    method='bilinear', overwrite=TRUE, format='CDF',
+    # args from writeRaster
+#    NAflag=regrid.datavar.na,
+    varname=regrid.datavar.name, 
+    varunit=regrid.datavar.unit,
+    longname=regrid.datavar.longname,
+    xname=x.var.name,
+    yname=y.var.name,
+    zname=z.var.name,
+    zunit=z.var.unit,
+    filename=out.fp)
+
+# # start debugging
+# out.raster
+# system(sprintf('ls -alht %s | head', work.dir))
+# system(sprintf('ncdump -h %s', out.fp))
+# netCDF.stats.to.stdout.by.timestep(
+#   netcdf.fp=out.fp, data.var.name=regrid.datavar.name, time.dim.name='time')
+# #   end debugging
+
+# > out.raster
+# class       : RasterBrick 
+# dimensions  : 299, 459, 137241, 14  (nrow, ncol, ncell, nlayers)
+# resolution  : 12000, 12000  (x, y)
+# extent      : -2556000, 2952000, -1728000, 1860000  (xmin, xmax, ymin, ymax)
+# coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000 
+# data source : /home/rtd/code/regridding/CLM_CN_global_to_AQMEII-NA/2008PTONCLMCNN2O_regrid.nc 
+# names       : X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14 
+# z-value     : 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 
+# varname     : n2oemissions 
+
+# note out.raster@extent in meters, zero-centered
+
+# > system(sprintf('ncdump -h %s', out.fp))
+# netcdf \2008PTONCLMCNN2O_regrid {
+# dimensions:
+#   COL = 459 ;
+#   ROW = 299 ;
+#   time = UNLIMITED ; // (14 currently)
+# variables:
+#   double COL(COL) ;
+#     COL:units = "meter" ;
+#     COL:long_name = "COL" ;
+#   double ROW(ROW) ;
+#     ROW:units = "meter" ;
+#     ROW:long_name = "ROW" ;
+#   int time(time) ;
+#     time:units = "month (1=2, 2=Jan, ... 13=Dec, 13=14)" ;
+#     time:long_name = "time" ;
+#   float n2oemissions(time, ROW, COL) ;
+#     n2oemissions:units = "mgN/m2/month" ;
+#     n2oemissions:_FillValue = -3.4e+38 ;
+#     n2oemissions:missing_value = -3.4e+38 ;
+#     n2oemissions:long_name = "N2O emissions" ;
+#     n2oemissions:projection = "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000" ;
+#     n2oemissions:projection_format = "PROJ.4" ;
+#     n2oemissions:min = 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. ;
+#     n2oemissions:max = 5.11556934631735, 5.11556934631735, 7.5402387212348, 33.9675487444286, 35.576195359168, 29.08862056211, 149.982927003544, 174.181162633967, 44.017017369279, 53.2281713126664, 25.7052935681188, 14.5454227669125, 15.1988578410051, 15.1988578410051 ;
+
+# ----------------------------------------------------------------------
+# visualize output
+# ----------------------------------------------------------------------
+
+## project North American map
+
+map.noram.shp.proj <- rgdal::spTransform(map.noram.shp.unproj, CRS=out.crs)
+# TODO: gotta have SpatialLines, not SpatialPolygons
+# # start debugging
+# class(map.noram.shp.proj)
+# bbox(map.noram.shp.proj)
+# #   end debugging
+
+# > class(map.noram.shp.proj)
+# [1] "SpatialPolygonsDataFrame"
+# attr(,"package")
+# [1] "sp"
+# > bbox(map.noram.shp.proj)
+#        min     max
+# x -6930127 3233667
+# y -2879573 5675464
+
+# compare to (above)
+# > template.extents
+# class       : Extent 
+# xmin        : -2556000 
+# xmax        : 2952000 
+# ymin        : -1728000 
+# ymax        : 1860000 
+
+## Get US state boundaries in projection units
+
+# see http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot
+state.map <- maps::map(
+  database="state", projection="lambert", par=c(33,45), plot=FALSE)
+#                  parameters to lambert: ^^^^^^^^^^^^
+#                  see mapproj::mapproject
+# note class(state.map) == "map"
+
+# replace its coordinates with the following:
+# metadata.coords.IOAPI.list <- M3::get.grid.info.M3(out.fp)
+# no, out.fp lacks the IOAPI metadata
+metadata.coords.IOAPI.list <- M3::get.grid.info.M3(template.in.fp)
+metadata.coords.IOAPI.x.orig <- metadata.coords.IOAPI.list$x.orig
+metadata.coords.IOAPI.y.orig <- metadata.coords.IOAPI.list$y.orig
+metadata.coords.IOAPI.x.cell.width <- metadata.coords.IOAPI.list$x.cell.width
+metadata.coords.IOAPI.y.cell.width <- metadata.coords.IOAPI.list$y.cell.width
+state.map.lines <- M3::get.map.lines.M3.proj(
+#  file=out.fp, database="state", units="m")
+  file=template.in.fp, database="state", units="m")
+state.map.lines.coords.IOAPI.x <-
+  (state.map.lines$coords[,1] - metadata.coords.IOAPI.x.orig)
+# no! here, dimensions are in meters, not cells
+#    /metadata.coords.IOAPI.x.cell.width
+state.map.lines.coords.IOAPI.y <-
+  (state.map.lines$coords[,2] - metadata.coords.IOAPI.y.orig)
+#    /metadata.coords.IOAPI.y.cell.width
+state.map.lines.coords.IOAPI <- 
+  cbind(state.map.lines.coords.IOAPI.x, state.map.lines.coords.IOAPI.y)
+
+# start debugging
+# class(state.map.lines.coords.IOAPI)
+# summary(state.map.lines.coords.IOAPI)
+# #   end debugging
+
+# > class(state.map.lines.coords.IOAPI)
+# [1] "matrix"
+# > summary(state.map.lines.coords.IOAPI)
+#  state.map.lines.coords.IOAPI.x state.map.lines.coords.IOAPI.y
+#  Min.   : 283762                Min.   : 160844               
+#  1st Qu.:2650244                1st Qu.:1054047               
+#  Median :3469204                Median :1701052               
+#  Mean   :3245997                Mean   :1643356               
+#  3rd Qu.:4300969                3rd Qu.:2252531               
+#  Max.   :4878260                Max.   :2993778               
+#  NA's   :168                    NA's   :168             
+
+# out.raster@extent in meters, zero-centered:
+# extent : -2556000, 2952000, -1728000, 1860000  (xmin, xmax, ymin, ymax)
+# derive out.raster@range(x)==5508000, range(y)==3588000
+# but map@extent is not zero-based. Sooo ... just add (xmin, ymin) ???
+
+state.map.IOAPI <- state.map # copy
+state.map.IOAPI$x <- state.map.lines.coords.IOAPI.x + extents.xmin
+state.map.IOAPI$y <- state.map.lines.coords.IOAPI.y + extents.ymin
+state.map.IOAPI$range <- c(
+  min(state.map.IOAPI$x),
+  max(state.map.IOAPI$x),
+  min(state.map.IOAPI$y),
+  max(state.map.IOAPI$y))
+state.map.IOAPI.shp <-
+  maptools::map2SpatialLines(state.map.IOAPI, proj4string=out.crs)
+# # start debugging
+# class(state.map.IOAPI.shp)
+# bbox(state.map.IOAPI.shp)
+# # summary(do.call("rbind",    # thanks, Felix Andrews!
+# #   unlist(coordinates(state.map.IOAPI.shp), recursive=FALSE)))
+# #   end debugging
+
+# > class(state.map.IOAPI.shp)
+# [1] "SpatialLines"
+# attr(,"package")
+# [1] "sp"
+# > bbox(state.map.IOAPI.shp)
+#        min     max
+# x -2272238 2322260
+# y -1567156 1265778
+
+# compare to (above)
+# > template.extents
+# class       : Extent 
+# xmin        : -2556000 
+# xmax        : 2952000 
+# ymin        : -1728000 
+# ymax        : 1860000 
+
+# KLUDGE!
+map.shp <- state.map.IOAPI.shp
+#map.shp <- map.noram.shp.proj
+
+visualize.layers(
+  nc.fp=out.fp,
+  brick=out.raster,
+  datavar.name=regrid.datavar.name,
+  layer.dim.name=time.dim.name,
+#  map.shp=state.map.IOAPI.shp,
+  map.shp=map.noram.shp.proj,
+  pdf.fp=out.pdf.fp,
+  pdf.height=25,
+  pdf.width=5
+)

regrid_global_to_AQMEII.sh

+#!/usr/bin/env bash
+
+# Requires R.
+# Driver for [regrid_global_to_AQMEII.r,retemp_monthly_to_hourly.ncl], converting
+# * global/unprojected CLM-CN data to AQMEII-NA, an LCC-projected subdomain
+# * monthly to hourly data (over each month)
+# Configure as needed for your platform.
+
+# ----------------------------------------------------------------------
+# constants with some simple manipulations
+# ----------------------------------------------------------------------
+
+THIS_FN="$(basename $0)"
+# will do the real work
+export CALL_REGRID_FN='regrid_global_to_AQMEII.r'
+export CALL_RETEMP_FN='retemp_monthly_to_hourly.ncl'
+
+# for R on EMVL: don't ask :-(
+R_DIR='/usr/local/apps/R-2.15.2/intel-13.0/bin'
+R="${R_DIR}/R"
+RSCRIPT="${R_DIR}/Rscript"
+
+# workspace
+export WORK_DIR="$(pwd)" # keep it simple for now: same dir as top of repo
+
+# will do the real work
+export CALL_REGRID_FP="${WORK_DIR}/${CALL_REGRID_FN}"
+export CALL_RETEMP_FP="${WORK_DIR}/${CALL_RETEMP_FN}"
+
+# for visualization (generally)
+export OUTPUT_SIGNIFICANT_DIGITS='3'
+
+# for plotting (specifically)
+PDF_VIEWER='xpdf'  # whatever works on your platform
+# temporally disaggregate multiple plots
+DATE_FORMAT='%Y%m%d_%H%M'
+export PDF_DIR="${WORK_DIR}"
+
+# Raw input data is global/unprojected CLM-CN netCDF, with coordvars=
+# lat==(-90,+90)
+# lon==(0,+360), hence need to
+export ROTATE_INPUT="true"
+# time: 14 months (not "decimal year" per `time:units`:
+#       month[0]== month[1], month[12]== month[13] (per Saikawa)
+# Get from my repository.
+CLM_CN_RAW_URI='https://bitbucket.org/tlroche/clm_cn_global_to_aqmeii-na/downloads/2008PTONCLMCNN2O.nc'
+CLM_CN_RAW_FN="$(basename ${CLM_CN_RAW_URI})"
+CLM_CN_RAW_FN_ROOT="${CLM_CN_RAW_FN%.*}" # everything left of the dot
+export CLM_CN_RAW_FP="${WORK_DIR}/${CLM_CN_RAW_FN}"
+
+CLM_CN_RAW_PDF_FN="${CLM_CN_RAW_FN_ROOT}_$(date +${DATE_FORMAT}).pdf"
+export CLM_CN_RAW_PDF_FP="${WORK_DIR}/${CLM_CN_RAW_PDF_FN}" # path to PDF output
+
+# TODO: automate getting this metadata from the netCDF
+export CLM_CN_RAW_DATAVAR_NAME='n2oemissions'
+export CLM_CN_RAW_DATAVAR_LONGNAME='N2O emissions'
+export CLM_CN_RAW_DATAVAR_UNIT='mgN/m2/month'
+# export CLM_CN_RAW_DATAVAR_NA='-999.0' # not used--hopefully they have no missing data
+export CLM_CN_RAW_TIME_VAR_NAME='time'
+
+# Template for `raster` regridding.
+
+# The template input is a copy of some meteorological data with 52 "real" datavars.
+# That data is in the IOAPI format, which here is basically a wrapper around netCDF.
+# (IOAPI can be used with other data, and similar wrappers exist for other models.)
+# I removed all but one of the datavars (with NCO 'ncks'). TODO: script that!
+TEMPLATE_INPUT_URI='https://bitbucket.org/tlroche/geia_regrid/downloads/emis_mole_all_20080101_12US1_cmaq_cb05_soa_2008ab_08c.EXTENTS_INPUT.nc' # no need to make another 14 MB copy
+TEMPLATE_INPUT_FN="$(basename ${TEMPLATE_INPUT_URI})"
+export TEMPLATE_INPUT_FP="${WORK_DIR}/${TEMPLATE_INPUT_FN}"
+# export TEMPLATE_INPUT_BAND='1' # `ncks` makes dim=TSTEP first
+export TEMPLATE_VAR_NAME='emi_n2o'
+
+# Intermediate product is netCDF regridded to AQMEII-NA
+
+CLM_CN_REGRID_URI='https://bitbucket.org/tlroche/clm_cn_global_to_aqmeii-na/downloads/2008PTONCLMCNN2O_regrid.nc'
+CLM_CN_REGRID_FN="$(basename ${CLM_CN_REGRID_URI})"
+CLM_CN_REGRID_FN_ROOT="${CLM_CN_REGRID_FN%.*}" # everything left of the dot
+export CLM_CN_REGRID_FP="${WORK_DIR}/${CLM_CN_REGRID_FN}"
+
+export CLM_CN_REGRID_X_VAR_NAME='COL'
+export CLM_CN_REGRID_Y_VAR_NAME='ROW'
+export CLM_CN_REGRID_TIME_VAR_NAME="${CLM_CN_RAW_TIME_VAR_NAME}"
+export CLM_CN_REGRID_TIME_VAR_UNIT='month (1=2, 2=Jan, ... 13=Dec, 13=14)'
+
+CLM_CN_REGRID_PDF_FN="${CLM_CN_REGRID_FN_ROOT}_$(date +${DATE_FORMAT}).pdf"
+export CLM_CN_REGRID_PDF_FP="${WORK_DIR}/${CLM_CN_REGRID_PDF_FN}" # path to PDF output
+
+# These helpers should have been cloned into the cwd, but JIC (and for later use). TODO: R-package my code
+
+# STAT_SCRIPT_URI='https://bitbucket.org/tlroche/geia_regrid/raw/6ecdb8ae7d6fd661f7c4616a14b817c5aa9c7f90/netCDF.stats.to.stdout.r'
+# STAT_SCRIPT_FN="$(basename ${STAT_SCRIPT_URI})"
+# export STAT_SCRIPT_FP="${WORK_DIR}/${STAT_SCRIPT_FN}"
+
+# PLOT_SCRIPT_URI='https://bitbucket.org/tlroche/geia_regrid/raw/6ecdb8ae7d6fd661f7c4616a14b817c5aa9c7f90/plotLayersForTimestep.r'
+# PLOT_SCRIPT_FN="$(basename ${PLOT_SCRIPT_URI})"
+# export PLOT_SCRIPT_FP="${WORK_DIR}/${PLOT_SCRIPT_FN}"
+
+# ----------------------------------------------------------------------
+# setup
+# ----------------------------------------------------------------------
+
+# for netcdf4 (ncdf4.so, libnetcdf.so.7) on EMVL:
+# I may need to run this from CLI, not here, and
+# you may not need this at all!
+module add netcdf-4.1.2
+
+mkdir -p ${WORK_DIR}
+
+if [[ -z "${CLM_CN_RAW_FP}" ]] ; then
+  echo -e "${THIS_FN}: ERROR: CLM_CN_RAW_FP not defined"
+  exit 4
+fi
+if [[ ! -r "${CLM_CN_RAW_FP}" ]] ; then
+  for CMD in \
+    "wget --no-check-certificate -c -O ${CLM_CN_RAW_FP} ${CLM_CN_RAW_URI}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+if [[ -z "${TEMPLATE_INPUT_FP}" ]] ; then
+  echo -e "${THIS_FN}: ERROR: TEMPLATE_INPUT_FP not defined"
+  exit 5
+fi
+if [[ ! -r "${TEMPLATE_INPUT_FP}" ]] ; then
+  for CMD in \
+    "wget --no-check-certificate -c -O ${TEMPLATE_INPUT_FP} ${TEMPLATE_INPUT_URI}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+# only for testing. TODO: flag control
+if [[ -z "${CLM_CN_REGRID_FP}" ]] ; then
+  echo -e "${THIS_FN}: ERROR: CLM_CN_REGRID_FP not defined"
+  exit 7
+fi
+# if [[ ! -r "${CLM_CN_REGRID_FP}" ]] ; then
+#   for CMD in \
+#     "wget --no-check-certificate -c -O ${CLM_CN_REGRID_FP} ${CLM_CN_REGRID_URI}" \
+#   ; do
+#     echo -e "$ ${CMD}"
+#     eval "${CMD}"
+#   done
+# fi
+
+# ----------------------------------------------------------------------
+# payload
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# convert from global extents to AQMEII
+# ----------------------------------------------------------------------
+
+# If this fails ...
+"${RSCRIPT}" "${CALL_REGRID_FP}"
+# ... just start R ...
+# "${R}"
+# ... and source ${CALL_REGRID_FP}, e.g.,
+# source('....r')
+
+# After exiting R, show cwd and display PDFs (if plotted)
+if [[ -r "${CLM_CN_RAW_PDF_FP}" ]] ; then
+  if [[ -r "${CLM_CN_REGRID_PDF_FP}" ]] ; then
+#      "${PDF_VIEWER} ${CLM_CN_RAW_PDF_FP} &" \
+    for CMD in \
+      "ls -alht ${WORK_DIR}" \
+      "${PDF_VIEWER} ${CLM_CN_REGRID_PDF_FP} &" \
+    ; do
+      echo -e "$ ${CMD}"
+      eval "${CMD}"
+    done
+  else
+    echo -e "${THIS_FN}: ERROR: regrid not plotted: ${CLM_CN_REGRID_PDF_FP}"
+  fi
+else
+  echo -e "${THIS_FN}: ERROR: input not plotted: ${CLM_CN_RAW_PDF_FP}"
+fi