Commits

Tom Roche  committed cd4e901

rename, tweak prior to reuse/refactoring on CMAQ side

  • Participants
  • Parent commits 8074764

Comments (0)

Files changed (8)

File regrid.surface.elevation.r

-# R code to regrid pre-cropped surface elevation netCDF data to match other data over the same extents (but diferently gridded).
-# TODO: make generic, pass parameters
-
-# constants-----------------------------------------------------------
-
-cat('Setting global constants ...\n') # debugging
-
-# global constants, mostly from environment---------------------------
-
-cat('Loading package=raster ...\n') # debugging
-library(raster)
-# coordinate reference system: here, just global lon-lat
-global.proj <- CRS('+proj=longlat +ellps=WGS84')
-
-# all the following env vars must be set and exported in driver script (TODO: fail if not)
-this.fn <- Sys.getenv('CALL_R_FN')  # FIXME: gets overidden by `source` below
-cat(sprintf('this.fn=%s\n', this.fn)) # debugging
-work.dir <- Sys.getenv('WORK_DIR')
-cat(sprintf('work.dir=%s\n', work.dir)) # debugging
-
-sigdigs <- as.numeric(Sys.getenv('SIGDIGS')) # significant digits for lattice strips
-
-stat.script.fp <- Sys.getenv('STAT_SCRIPT_FP')
-source(stat.script.fp) # produces errant error=
-#> netCDF.stats.to.stdout.r: no arguments supplied, exiting
-# unnecessary here
-# plot.script.fp <- Sys.getenv('PLOT_SCRIPT_FP')
-# source(plot.script.fp)
-
-# local constants, from environment-----------------------------------
-# taking from Rscript call failed, annoyingly and mysteriously :-(
-
-cat('Setting local constants ...\n') # debugging
-
-in.fp <- Sys.getenv('IN_FP')
-in.datavar.name <- Sys.getenv('IN_DATAVAR_NAME')
-in.datavar.band <- Sys.getenv('IN_DATAVAR_BAND')
-in.crs <- global.proj # TODO: get from file!
-
-template.fp <- Sys.getenv('TEMPLATE_FP')
-template.datavar.name <- Sys.getenv('TEMPLATE_DATAVAR_NAME')
-
-raster.brick <- Sys.getenv('RASTER_BRICK')
-raster.rotate <- Sys.getenv('RASTER_ROTATE')
-raster.plot <- Sys.getenv('RASTER_PLOT')
-
-out.fp <- Sys.getenv('OUT_FP')
-out.crs <- in.crs
-out.datavar.name <- Sys.getenv('OUT_DATAVAR_NAME')
-out.datavar.type <- Sys.getenv('OUT_DATAVAR_TYPE')
-out.datavar.band <- Sys.getenv('OUT_DATAVAR_BAND')
-out.datavar.na <- Sys.getenv('OUT_DATAVAR_NA')
-out.datavar.unit <- Sys.getenv('OUT_DATAVAR_UNIT')
-out.datavar.longname <- Sys.getenv('OUT_DATAVAR_LONGNAME')
-out.datavar.coord.x.name <- Sys.getenv('OUT_DATAVAR_COORD_X_NAME')
-out.datavar.coord.y.name <- Sys.getenv('OUT_DATAVAR_COORD_Y_NAME')
-# out.datavar.coord.z.name <- Sys.getenv('OUT_DATAVAR_COORD_Z_NAME')
-# out.datavar.coord.z.unit <- Sys.getenv('OUT_DATAVAR_COORD_Z_UNIT')
-
-out.pdf.fp <- Sys.getenv('OUT_PDF_FP')
-out.pdf.height <- Sys.getenv('OUT_PDF_HEIGHT')
-out.pdf.width <- Sys.getenv('OUT_PDF_WIDTH')
-
-# # start debugging
-# # TODO: check paths
-# cat(sprintf('in.fp==%s\n', in.fp))
-# cat(sprintf('in.datavar.name==%s\n', in.datavar.name))
-# cat(sprintf('in.datavar.band==%s\n', in.datavar.band))
-
-# cat(sprintf('raster.brick==%s\n', raster.brick))
-# cat(sprintf('raster.rotate==%s\n', raster.rotate))
-# cat(sprintf('raster.plot==%s\n', raster.plot))
-
-# cat(sprintf('template.fp==%s\n', template.fp))
-# cat(sprintf('template.datavar.name==%s\n', template.datavar.name))
-
-# cat(sprintf('out.fp==%s\n', out.fp))
-# cat(sprintf('out.datavar.name==%s\n', out.datavar.name))
-# cat(sprintf('out.datavar.type==%s\n', out.datavar.type))
-# cat(sprintf('out.datavar.band==%s\n', out.datavar.band))
-# cat(sprintf('out.datavar.na==%s\n', out.datavar.na))
-# cat(sprintf('out.datavar.unit==%s\n', out.datavar.unit))
-# cat(sprintf('out.datavar.longname==%s\n', out.datavar.longname))
-# cat(sprintf('out.datavar.coord.x.name==%s\n', out.datavar.coord.x.name))
-# cat(sprintf('out.datavar.coord.y.name==%s\n', out.datavar.coord.y.name))
-# # cat(sprintf('out.datavar.coord.z.name==%s\n', out.datavar.coord.z.name))
-# # cat(sprintf('out.datavar.coord.z.unit==%s\n', out.datavar.coord.z.unit))
-# cat(sprintf('out.pdf.fp==%s\n', out.pdf.fp))
-# cat(sprintf('out.pdf.height==%s\n', out.pdf.height))
-# cat(sprintf('out.pdf.width==%s\n', out.pdf.width))
-# #   end debugging
-
-raster.brick <- as.logical(raster.brick)
-raster.rotate <- as.logical(raster.rotate)
-raster.plot <- as.logical(raster.plot)
-in.datavar.band <- as.numeric(in.datavar.band)
-out.datavar.na <- as.numeric(out.datavar.na)
-out.datavar.band <- as.numeric(out.datavar.band)
-out.pdf.height <- as.numeric(out.pdf.height)
-out.pdf.width <- as.numeric(out.pdf.width)
-
-# ensure global scope, so raster.aqmeii can be accessed inside raster.plot conditional? nope, fail
-# raster.crop <- TRUE # kludge
-
-# functions-----------------------------------------------------------
-
-# code----------------------------------------------------------------
-
-# 1. Load the input datavar
-
-if (raster.brick) { # multiple vertical layers
-  in.raster <- brick(in.fp, varname=in.datavar.name, crs=global.proj)
-  layers.n <- nbands(in.raster)
-} else {
-  in.raster <- raster(in.fp, varname=in.datavar.name, crs=global.proj)
-  layers.n <- 1
-}
-# start debugging
-in.raster 
-# e.g.
-# class       : RasterLayer 
-# dimensions  : 2221, 4351, 9663571  (nrow, ncol, ncell)
-# resolution  : 0.01666667, 0.01666667  (x, y)
-# extent      : -131.0083, -58.49167, 22.49167, 59.50833  (xmin, xmax, ymin, ymax)
-# coord. ref. : +proj=longlat +datum=WGS84 
-# data source : ./ETOPO1_Ice_g_gmt4.grd_region_zS.nc 
-# names       : surface.elevation 
-# zvar        : zS 
-
-netCDF.stats.to.stdout(netcdf.fp=in.fp, var.name=in.datavar.name)
-# e.g.
-# For ./ETOPO1_Ice_g_gmt4.grd_region_zS.nc var=zS
-#       cells=9663571
-#       obs=9663571
-#       min=-6.63e+03
-#       q1=-1.36e+03
-#       med=203
-#       mean=-659
-#       q3=569
-#       max=4.16e+03
-#   end debugging
-
-# 2. Rotate the input global datavar if requested:
-#    zero-based longitudes cause problems for wrld_simpl (and probably other data)
-if (raster.rotate) {
-  in.raster <- rotate(in.raster,
-    overwrite=TRUE) # else levelplot does one layer per page?
-  in.raster # debugging
-}
-
-# 3. Create template raster.
-template.raster <-projectExtent(
-  raster(template.fp, varname=template.datavar.name), crs=out.crs)
-# template.raster <- projectExtent(template.in.raster, crs=out.crs)
-# start debugging
-template.raster 
-# e.g.
-# class       : RasterLayer 
-# dimensions  : 21, 30, 630  (nrow, ncol, ncell)
-# resolution  : 2.5, 1.804511  (x, y)
-# extent      : -131.25, -56.25, 22.73684, 60.63158  (xmin, xmax, ymin, ymax)
-# coord. ref. : +proj=longlat +ellps=WGS84 
-
-netCDF.stats.to.stdout(netcdf.fp=template.fp, var.name=template.datavar.name)
-# e.g.
-# For ./2008N2O_restart_region_PS.nc var=PS
-#       cells=630
-#       obs=630
-#       min=7.43e+04
-#       q1=9.59e+04
-#       med=9.93e+04
-#       mean=9.74e+04
-#       q3=1.02e+05
-#       max=1.03e+05
-#   end debugging
-
-# 4. Regrid input using template.
-out.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.crs,
-    method='bilinear', overwrite=TRUE, format='CDF',
-    # args from writeRaster
-    # datavar type will truncate unless set?
-    datatype=out.datavar.type,
-    NAflag=out.datavar.na,
-    varname=out.datavar.name, 
-    # netCDF-specific arguments
-    varunit=out.datavar.unit,
-    longname=out.datavar.longname,
-    xname=out.datavar.coord.x.name,
-    yname=out.datavar.coord.y.name,
-#    zname=out.datavar.coord.z.name,
-#    zunit=out.datavar.coord.z.unit
-    filename=out.fp)
-# if above fails to set CRS
-# out.raster@crs <- CRS(out.crs)
-
-# start debugging
-# regrid.end.time <- system('date', intern=TRUE)
-# cat(sprintf('  end regrid @ %s\n', regrid.end.time))
-# cat(sprintf('%s\n', regrid.start.time))
-
-out.raster 
-# e.g.
-# class       : RasterLayer 
-# dimensions  : 21, 30, 630  (nrow, ncol, ncell)
-# resolution  : 2.5, 1.804511  (x, y)
-# extent      : -131.25, -56.25, 22.73684, 60.63158  (xmin, xmax, ymin, ymax)
-# coord. ref. : +proj=longlat +ellps=WGS84 
-# data source : ./ETOPO1_Ice_g_gmt4.grd_region_zS_regrid.nc 
-# names       : surface.elevation 
-# zvar        : zS 
-
-netCDF.stats.to.stdout(netcdf.fp=out.fp, var.name=out.datavar.name)
-# e.g.
-# For ./ETOPO1_Ice_g_gmt4.grd_region_zS_regrid.nc var=zS
-#       cells=630
-#       obs=580
-#       min=-6.16e+03
-#       q1=-1.54e+03
-#       med=200
-#       mean=-658
-#       q3=565
-#       max=3.82e+03
-# problem! significant decrease from orig max=4.16e+03, but this just regrids over same region ...
-
-system(sprintf('ls -alht %s | head', work.dir))
-system(sprintf('ncdump -h %s', out.fp))
-#   end debugging
-
-# 5. Plot data if requested
-
-# 5.0 Massive annoyance: out.raster is out of scope inside this conditional!
-# Despite numerous attempted kludges, I cannot make plotting conditional!
-# if (raster.plot) {
-
-# 5.1 Fix layer names (which are plot labels in `layerplot`)
-
-if (raster.brick) {
-  # TODO: since using rasterVis::levelplot, get/set names with `layerNames`
-  aqmeii.df <- as.data.frame(out.raster)
-  names(aqmeii.df) <- formatC(
-    as.numeric(sub(x=names(aqmeii.df), pattern="X", replacement="", fixed=TRUE)),
-    format="g", digits=sigdigs)
-  # start debugging
-  cat(sprintf('layers.n==%s\n', layers.n))
-  ## out.raster
-  ## aqmeii.df
-  ## cat('head and tail of aqmeii.df:\n')
-  ## head(aqmeii.df)
-  ## tail(aqmeii.df)
-  cat('head and tail of names(aqmeii.df):\n')
-  head(names(aqmeii.df))
-  tail(names(aqmeii.df))
-  #  end debugging
-} # end if (raster.brick)
-
-# 5.2 Create world map (which rasterVis will autocrop)
-
-  cat('Loading map packages ...\n') # debugging
-  library(maps)
-  library(maptools) # needed for map2SpatialLines
-  data(wrld_simpl)      # from maptools
-
-  # TODO: save/restore after first creation!
-  global.map <- map('world', plot=FALSE)
-  global.map.shp <- map2SpatialLines(global.map, proj4string=global.proj)
-
-# 5.3 Plot and display.
-
-  cat('Loading plot packages ...\n') # debugging
-  library(rasterVis)
-
-  cat(sprintf(
-    '%s: plotting %s may take awhile! (TODO: add progress control)\n',
-    this.fn, out.pdf.fp))
-#  if (raster.plot) {
-  pdf(file=out.pdf.fp, width=out.pdf.width, height=out.pdf.height)
-
-if (raster.brick) {
-    rasterVis::levelplot(out.raster,
-      margin=FALSE,
-      names.attr=names(aqmeii.df),
-      layout=c(1,layers.n), # all layers in one "atmospheric" column
-    ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
-  } else {
-    rasterVis::levelplot(out.raster,
-      margin=FALSE
-    ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
-    # the plot crops the map!
-  } # end if (raster.brick)
-  # teardown
-  dev.off()
-
-# } # end if (raster.plot)
-
-system(sprintf('ls -alh %s', out.pdf.fp)) # debugging
-

File regrid_surface_elevation.sh

-#!/usr/bin/env bash
-# Requires R.
-# Driver for regrid.surface.elevation.r--configure as needed for your platform.
-# Crops the following netCDF {file, datavar} pairs, saving output to separate file:
-# * 2008N2O_restart.nc,N2O
-# * 2008N2O_restart.nc,PS
-# * ETOPO1_Ice_g_gmt4.grd.nc,zS: after renaming
-# ** x->lon (to match 2008N2O_restart.nc)
-# ** y->lat (to match 2008N2O_restart.nc)
-# ** z->zS  (à la 2008N2O_restart.nc:PS)
-
-# constants with some simple manipulations------------------------------
-
-THIS_FN="$(basename $0)"
-# will do the real work
-# TODO: create from THIS_FN, e.g., use `sed -e s/_/./g`
-export CALL_R_FN='regrid.surface.elevation.r'
-
-# modules on EMVL. TODO: fixme! (these fail from here, are currently just reminders.)
-# You may not need them at all!
-module add netcdf-4.1.2
-#module add ncl_ncarg-6.0.0
-#module add nco-4.0.5
-
-# for R on EMVL: don't ask :-(
-R_DIR='/usr/local/R-2.15.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 (restricting and plotting)
-export CALL_R_FP="${WORK_DIR}/${CALL_R_FN}"
-
-# for data display
-export SIGDIGS='4' # significant digits
-
-# for plotting
-PDF_VIEWER='xpdf'  # whatever works on your platform
-# temporally disaggregate multiple plots
-DATE_FORMAT='%Y%m%d_%H%M'
-export PDF_DIR="${WORK_DIR}"
-
-# global lon-lat netCDFs: get from my repository, or download original if available
-
-# template file (with the grid which we want the output to match)
-TEMPLATE_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart_region_PS.nc'
-TEMPLATE_FN="$(basename ${TEMPLATE_URI})"
-# get both parts of the filename (before and after '.')
-TEMPLATE_FN_ROOT="${TEMPLATE_FN%.*}"
-TEMPLATE_FN_EXT="${TEMPLATE_FN##*.}"
-export TEMPLATE_FP="${WORK_DIR}/${TEMPLATE_FN}"
-export TEMPLATE_DATAVAR_NAME='PS'
-
-# input: surface elevation from ETOPO1 cropped with restrict_lon_lat_to_region.sh
-IN_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/ETOPO1_Ice_g_gmt4.grd_region_zS.nc'
-IN_FN="$(basename ${IN_URI})"
-IN_FN_ROOT="${IN_FN%.*}"
-IN_FN_EXT="${IN_FN##*.}"
-export IN_FP="${WORK_DIR}/${IN_FN}"
-export IN_DATAVAR_NAME='zS'
-# TODO: get from file, e.g., with NCL
-export IN_DATAVAR_NA='-2147483648'
-export IN_DATAVAR_BAND='NA'
-IN_DATAVAR_UNIT='m'
-IN_DATAVAR_LONGNAME='surface elevation'
-IN_DATAVAR_COORD_X_NAME='lon'
-IN_DATAVAR_COORD_Y_NAME='lat'
-
-export RASTER_BRICK='FALSE'
-export RASTER_ROTATE='FALSE'
-export RASTER_PLOT='TRUE'  # can't make R use this???
-
-# output: cropped surface elevation from ETOPO1 regridded to match template
-OUT_FN_ROOT="${IN_FN_ROOT}_regrid"
-OUT_FN="${OUT_FN_ROOT}.${IN_FN_EXT}"
-export OUT_FP="${WORK_DIR}/${OUT_FN}"
-export OUT_DATAVAR_NAME="${IN_DATAVAR_NAME}"
-# note type=int in original, but probably should be float for regridding
-# raster::writeRaster seems to autocast down unless instructed otherwise?
-export OUT_DATAVAR_TYPE='FLT4S' # == float: see raster.pdf::dataType
-export OUT_DATAVAR_NA="${IN_DATAVAR_NA}"
-export OUT_DATAVAR_BAND="${IN_DATAVAR_BAND}"
-export OUT_DATAVAR_UNIT="${IN_DATAVAR_UNIT}"
-export OUT_DATAVAR_LONGNAME="${IN_DATAVAR_LONGNAME}"
-export OUT_DATAVAR_COORD_X_NAME="${IN_DATAVAR_COORD_X_NAME}"
-export OUT_DATAVAR_COORD_Y_NAME="${IN_DATAVAR_COORD_Y_NAME}"
-# export OUT_DATAVAR_COORD_Z_NAME="${IN_DATAVAR_COORD_Z_NAME}"
-# export OUT_DATAVAR_COORD_Z_UNIT="${IN_DATAVAR_COORD_Z_UNIT}"
-
-OUT_PDF_FN="${OUT_FN_ROOT}_$(date +${DATE_FORMAT}).pdf"
-export OUT_PDF_FP="${PDF_DIR}/${OUT_PDF_FN}" # path to PDF output
-export OUT_PDF_HEIGHT='15'
-export OUT_PDF_WIDTH='20'
-
-# called scripts. TODO: R-package my code!
-
-# this script does simple statistics, and is copied from ioapi-hack-R (also on github)
-# TODO: R-package my code
-STAT_SCRIPT_URI='https://github.com/TomRoche/GEIA_to_netCDF/raw/master/netCDF.stats.to.stdout.r'
-STAT_SCRIPT_FN="$(basename ${STAT_SCRIPT_URI})"
-export STAT_SCRIPT_FP="${WORK_DIR}/${STAT_SCRIPT_FN}"
-
-# payload---------------------------------------------------------------
-
-# setup-----------------------------------------------------------------
-
-mkdir -p ${WORK_DIR}
-
-# get scripts-----------------------------------------------------------
-
-if [[ ! -r "${STAT_SCRIPT_FP}" ]] ; then
-  for CMD in \
-    "wget --no-check-certificate -c -O ${STAT_SCRIPT_FP} ${STAT_SCRIPT_URI}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# get data------------------------------------------------------------
-
-# template file
-if [[ ! -r "${TEMPLATE_FP}" ]] ; then
-  for CMD in \
-    "wget --no-check-certificate -c -O ${TEMPLATE_FP} ${TEMPLATE_URI}" \
-    "ncdump -h ${TEMPLATE_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# elevation source
-if [[ ! -r "${IN_FP}" ]] ; then
-  for CMD in \
-    "wget --no-check-certificate -c -O ${IN_FP} ${IN_URI}" \
-    "ncdump -h ${IN_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# call R script (the real work)---------------------------------------
-
-# TODO: pass commandline args to Rscript
-# punt: just use envvars :-(
-# for CMD in \
-#   "${RSCRIPT} ${CALL_R_FP}" \
-${RSCRIPT} ${CALL_R_FP}
-#   ; do
-#   echo -e "$ ${CMD}"
-#   eval "${CMD}"
-# done
-
-# After exiting R, show output files and display output PDF.
-for CMD in \
-  "ls -alht ${WORK_DIR} | head" \
-; do
-  echo -e "$ ${CMD}"
-  eval "${CMD}"
-done
-
-if [[ -r "${OUT_PDF_FP}" ]] ; then
-  for CMD in \
-    "${PDF_VIEWER} ${OUT_PDF_FP} &" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-else
-  echo -e "ERROR: ${THIS_FN}: PDF output file='${OUT_PDF_FP}' not readable"
-fi

File regrid_surface_elevation__MOZART.r

+# R code to regrid pre-cropped surface elevation netCDF data to match other data over the same extents (but diferently gridded).
+# TODO: make generic, pass parameters
+
+# constants-----------------------------------------------------------
+
+cat('Setting global constants ...\n') # debugging
+
+# global constants, mostly from environment---------------------------
+
+cat('Loading package=raster ...\n') # debugging
+library(raster)
+# coordinate reference system: here, just global lon-lat
+global.proj <- CRS('+proj=longlat +ellps=WGS84')
+
+# all the following env vars must be set and exported in driver script (TODO: fail if not)
+this.fn <- Sys.getenv('CALL_R_FN')  # FIXME: gets overidden by `source` below
+cat(sprintf('this.fn=%s\n', this.fn)) # debugging
+work.dir <- Sys.getenv('WORK_DIR')
+cat(sprintf('work.dir=%s\n', work.dir)) # debugging
+
+sigdigs <- as.numeric(Sys.getenv('SIGDIGS')) # significant digits for lattice strips
+
+stat.script.fp <- Sys.getenv('STAT_SCRIPT_FP')
+source(stat.script.fp) # produces errant error=
+#> netCDF.stats.to.stdout.r: no arguments supplied, exiting
+# unnecessary here
+# plot.script.fp <- Sys.getenv('PLOT_SCRIPT_FP')
+# source(plot.script.fp)
+
+# local constants, from environment-----------------------------------
+# taking from Rscript call failed, annoyingly and mysteriously :-(
+
+cat('Setting local constants ...\n') # debugging
+
+in.fp <- Sys.getenv('IN_FP')
+in.datavar.name <- Sys.getenv('IN_DATAVAR_NAME')
+in.datavar.band <- Sys.getenv('IN_DATAVAR_BAND')
+in.crs <- global.proj # TODO: get from file!
+
+template.fp <- Sys.getenv('TEMPLATE_FP')
+template.datavar.name <- Sys.getenv('TEMPLATE_DATAVAR_NAME')
+
+raster.brick <- Sys.getenv('RASTER_BRICK')
+raster.rotate <- Sys.getenv('RASTER_ROTATE')
+raster.plot <- Sys.getenv('RASTER_PLOT')
+
+out.fp <- Sys.getenv('OUT_FP')
+out.crs <- in.crs
+out.datavar.name <- Sys.getenv('OUT_DATAVAR_NAME')
+out.datavar.type <- Sys.getenv('OUT_DATAVAR_TYPE')
+out.datavar.band <- Sys.getenv('OUT_DATAVAR_BAND')
+out.datavar.na <- Sys.getenv('OUT_DATAVAR_NA')
+out.datavar.unit <- Sys.getenv('OUT_DATAVAR_UNIT')
+out.datavar.longname <- Sys.getenv('OUT_DATAVAR_LONGNAME')
+out.datavar.coord.x.name <- Sys.getenv('OUT_DATAVAR_COORD_X_NAME')
+out.datavar.coord.y.name <- Sys.getenv('OUT_DATAVAR_COORD_Y_NAME')
+# out.datavar.coord.z.name <- Sys.getenv('OUT_DATAVAR_COORD_Z_NAME')
+# out.datavar.coord.z.unit <- Sys.getenv('OUT_DATAVAR_COORD_Z_UNIT')
+
+out.pdf.fp <- Sys.getenv('OUT_PDF_FP')
+out.pdf.height <- Sys.getenv('OUT_PDF_HEIGHT')
+out.pdf.width <- Sys.getenv('OUT_PDF_WIDTH')
+
+# # start debugging
+# # TODO: check paths
+# cat(sprintf('in.fp==%s\n', in.fp))
+# cat(sprintf('in.datavar.name==%s\n', in.datavar.name))
+# cat(sprintf('in.datavar.band==%s\n', in.datavar.band))
+
+# cat(sprintf('raster.brick==%s\n', raster.brick))
+# cat(sprintf('raster.rotate==%s\n', raster.rotate))
+# cat(sprintf('raster.plot==%s\n', raster.plot))
+
+# cat(sprintf('template.fp==%s\n', template.fp))
+# cat(sprintf('template.datavar.name==%s\n', template.datavar.name))
+
+# cat(sprintf('out.fp==%s\n', out.fp))
+# cat(sprintf('out.datavar.name==%s\n', out.datavar.name))
+# cat(sprintf('out.datavar.type==%s\n', out.datavar.type))
+# cat(sprintf('out.datavar.band==%s\n', out.datavar.band))
+# cat(sprintf('out.datavar.na==%s\n', out.datavar.na))
+# cat(sprintf('out.datavar.unit==%s\n', out.datavar.unit))
+# cat(sprintf('out.datavar.longname==%s\n', out.datavar.longname))
+# cat(sprintf('out.datavar.coord.x.name==%s\n', out.datavar.coord.x.name))
+# cat(sprintf('out.datavar.coord.y.name==%s\n', out.datavar.coord.y.name))
+# # cat(sprintf('out.datavar.coord.z.name==%s\n', out.datavar.coord.z.name))
+# # cat(sprintf('out.datavar.coord.z.unit==%s\n', out.datavar.coord.z.unit))
+# cat(sprintf('out.pdf.fp==%s\n', out.pdf.fp))
+# cat(sprintf('out.pdf.height==%s\n', out.pdf.height))
+# cat(sprintf('out.pdf.width==%s\n', out.pdf.width))
+# #   end debugging
+
+raster.brick <- as.logical(raster.brick)
+raster.rotate <- as.logical(raster.rotate)
+raster.plot <- as.logical(raster.plot)
+in.datavar.band <- as.numeric(in.datavar.band)
+out.datavar.na <- as.numeric(out.datavar.na)
+out.datavar.band <- as.numeric(out.datavar.band)
+out.pdf.height <- as.numeric(out.pdf.height)
+out.pdf.width <- as.numeric(out.pdf.width)
+
+# ensure global scope, so raster.aqmeii can be accessed inside raster.plot conditional? nope, fail
+# raster.crop <- TRUE # kludge
+
+# functions-----------------------------------------------------------
+
+# code----------------------------------------------------------------
+
+# 1. Load the input datavar
+
+if (raster.brick) { # multiple vertical layers
+  in.raster <- brick(in.fp, varname=in.datavar.name, crs=global.proj)
+  layers.n <- nbands(in.raster)
+} else {
+  in.raster <- raster(in.fp, varname=in.datavar.name, crs=global.proj)
+  layers.n <- 1
+}
+# start debugging
+in.raster 
+# e.g.
+# class       : RasterLayer 
+# dimensions  : 2221, 4351, 9663571  (nrow, ncol, ncell)
+# resolution  : 0.01666667, 0.01666667  (x, y)
+# extent      : -131.0083, -58.49167, 22.49167, 59.50833  (xmin, xmax, ymin, ymax)
+# coord. ref. : +proj=longlat +datum=WGS84 
+# data source : ./ETOPO1_Ice_g_gmt4.grd_region_zS.nc 
+# names       : surface.elevation 
+# zvar        : zS 
+
+netCDF.stats.to.stdout(netcdf.fp=in.fp, var.name=in.datavar.name)
+# e.g.
+# For ./ETOPO1_Ice_g_gmt4.grd_region_zS.nc var=zS
+#       cells=9663571
+#       obs=9663571
+#       min=-6.63e+03
+#       q1=-1.36e+03
+#       med=203
+#       mean=-659
+#       q3=569
+#       max=4.16e+03
+#   end debugging
+
+# 2. Rotate the input global datavar if requested:
+#    zero-based longitudes cause problems for wrld_simpl (and probably other data)
+if (raster.rotate) {
+  in.raster <- rotate(in.raster,
+    overwrite=TRUE) # else levelplot does one layer per page?
+  in.raster # debugging
+}
+
+# 3. Create template raster.
+template.raster <-projectExtent(
+  raster(template.fp, varname=template.datavar.name), crs=out.crs)
+# template.raster <- projectExtent(template.in.raster, crs=out.crs)
+# start debugging
+template.raster 
+# e.g.
+# class       : RasterLayer 
+# dimensions  : 21, 30, 630  (nrow, ncol, ncell)
+# resolution  : 2.5, 1.804511  (x, y)
+# extent      : -131.25, -56.25, 22.73684, 60.63158  (xmin, xmax, ymin, ymax)
+# coord. ref. : +proj=longlat +ellps=WGS84 
+
+netCDF.stats.to.stdout(netcdf.fp=template.fp, var.name=template.datavar.name)
+# e.g.
+# For ./2008N2O_restart_region_PS.nc var=PS
+#       cells=630
+#       obs=630
+#       min=7.43e+04
+#       q1=9.59e+04
+#       med=9.93e+04
+#       mean=9.74e+04
+#       q3=1.02e+05
+#       max=1.03e+05
+#   end debugging
+
+# 4. Regrid input using template.
+out.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.crs,
+    method='bilinear', overwrite=TRUE, format='CDF',
+    # args from writeRaster
+    # datavar type will truncate unless set?
+    datatype=out.datavar.type,
+    NAflag=out.datavar.na,
+    varname=out.datavar.name, 
+    # netCDF-specific arguments
+    varunit=out.datavar.unit,
+    longname=out.datavar.longname,
+    xname=out.datavar.coord.x.name,
+    yname=out.datavar.coord.y.name,
+#    zname=out.datavar.coord.z.name,
+#    zunit=out.datavar.coord.z.unit
+    filename=out.fp)
+# if above fails to set CRS
+# out.raster@crs <- CRS(out.crs)
+
+# start debugging
+# regrid.end.time <- system('date', intern=TRUE)
+# cat(sprintf('  end regrid @ %s\n', regrid.end.time))
+# cat(sprintf('%s\n', regrid.start.time))
+
+out.raster 
+# e.g.
+# class       : RasterLayer 
+# dimensions  : 21, 30, 630  (nrow, ncol, ncell)
+# resolution  : 2.5, 1.804511  (x, y)
+# extent      : -131.25, -56.25, 22.73684, 60.63158  (xmin, xmax, ymin, ymax)
+# coord. ref. : +proj=longlat +ellps=WGS84 
+# data source : ./ETOPO1_Ice_g_gmt4.grd_region_zS_regrid.nc 
+# names       : surface.elevation 
+# zvar        : zS 
+
+netCDF.stats.to.stdout(netcdf.fp=out.fp, var.name=out.datavar.name)
+# e.g.
+# For ./ETOPO1_Ice_g_gmt4.grd_region_zS_regrid.nc var=zS
+#       cells=630
+#       obs=580
+#       min=-6.16e+03
+#       q1=-1.54e+03
+#       med=200
+#       mean=-658
+#       q3=565
+#       max=3.82e+03
+# problem! significant decrease from orig max=4.16e+03, but this just regrids over same region ...
+
+system(sprintf('ls -alht %s | head', work.dir))
+system(sprintf('ncdump -h %s', out.fp))
+#   end debugging
+
+# 5. Plot data if requested
+
+# 5.0 Massive annoyance: out.raster is out of scope inside this conditional!
+# Despite numerous attempted kludges, I cannot make plotting conditional!
+# if (raster.plot) {
+
+# 5.1 Fix layer names (which are plot labels in `layerplot`)
+
+if (raster.brick) {
+  # TODO: since using rasterVis::levelplot, get/set names with `layerNames`
+  aqmeii.df <- as.data.frame(out.raster)
+  names(aqmeii.df) <- formatC(
+    as.numeric(sub(x=names(aqmeii.df), pattern="X", replacement="", fixed=TRUE)),
+    format="g", digits=sigdigs)
+  # start debugging
+  cat(sprintf('layers.n==%s\n', layers.n))
+  ## out.raster
+  ## aqmeii.df
+  ## cat('head and tail of aqmeii.df:\n')
+  ## head(aqmeii.df)
+  ## tail(aqmeii.df)
+  cat('head and tail of names(aqmeii.df):\n')
+  head(names(aqmeii.df))
+  tail(names(aqmeii.df))
+  #  end debugging
+} # end if (raster.brick)
+
+# 5.2 Create world map (which rasterVis will autocrop)
+
+  cat('Loading map packages ...\n') # debugging
+  library(maps)
+  library(maptools) # needed for map2SpatialLines
+  data(wrld_simpl)      # from maptools
+
+  # TODO: save/restore after first creation!
+  global.map <- map('world', plot=FALSE)
+  global.map.shp <- map2SpatialLines(global.map, proj4string=global.proj)
+
+# 5.3 Plot and display.
+
+  cat('Loading plot packages ...\n') # debugging
+  library(rasterVis)
+
+  cat(sprintf(
+    '%s: plotting %s may take awhile! (TODO: add progress control)\n',
+    this.fn, out.pdf.fp))
+#  if (raster.plot) {
+  pdf(file=out.pdf.fp, width=out.pdf.width, height=out.pdf.height)
+
+if (raster.brick) {
+    rasterVis::levelplot(out.raster,
+      margin=FALSE,
+      names.attr=names(aqmeii.df),
+      layout=c(1,layers.n), # all layers in one "atmospheric" column
+    ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
+  } else {
+    rasterVis::levelplot(out.raster,
+      margin=FALSE
+    ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
+    # the plot crops the map!
+  } # end if (raster.brick)
+  # teardown
+  dev.off()
+
+# } # end if (raster.plot)
+
+system(sprintf('ls -alh %s', out.pdf.fp)) # debugging
+

File regrid_surface_elevation__MOZART.sh

+#!/usr/bin/env bash
+# Requires R.
+# Driver for regrid.surface.elevation.r--configure as needed for your platform.
+# Crops the following netCDF {file, datavar} pairs, saving output to separate file:
+# * 2008N2O_restart.nc,N2O
+# * 2008N2O_restart.nc,PS
+# * ETOPO1_Ice_g_gmt4.grd.nc,zS: after renaming
+# ** x->lon (to match 2008N2O_restart.nc)
+# ** y->lat (to match 2008N2O_restart.nc)
+# ** z->zS  (à la 2008N2O_restart.nc:PS)
+
+# constants with some simple manipulations------------------------------
+
+THIS_FN="$(basename $0)"
+# will do the real work
+export CALL_R_FN="$( echo -e ${THIS_FN} | sed -e 's/sh$/r/' )"
+
+# modules on EMVL. TODO: fixme! (these fail from here, are currently just reminders.)
+# You may not need them at all!
+module add netcdf-4.1.2
+#module add ncl_ncarg-6.0.0
+#module add nco-4.0.5
+
+# for R on EMVL: don't ask :-(
+R_DIR='/usr/local/R-2.15.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 (restricting and plotting)
+export CALL_R_FP="${WORK_DIR}/${CALL_R_FN}"
+
+# for data display
+export SIGDIGS='4' # significant digits
+
+# for plotting
+PDF_VIEWER='xpdf'  # whatever works on your platform
+# temporally disaggregate multiple plots
+DATE_FORMAT='%Y%m%d_%H%M'
+export PDF_DIR="${WORK_DIR}"
+
+# global lon-lat netCDFs: get from my repository, or download original if available
+
+# template file (with the grid which we want the output to match)
+TEMPLATE_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart_region_PS.nc'
+TEMPLATE_FN="$(basename ${TEMPLATE_URI})"
+# get both parts of the filename (before and after '.')
+TEMPLATE_FN_ROOT="${TEMPLATE_FN%.*}"
+TEMPLATE_FN_EXT="${TEMPLATE_FN##*.}"
+export TEMPLATE_FP="${WORK_DIR}/${TEMPLATE_FN}"
+export TEMPLATE_DATAVAR_NAME='PS'
+
+# input: surface elevation from ETOPO1 cropped with restrict_lon_lat_to_region.sh
+IN_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/ETOPO1_Ice_g_gmt4.grd_region_zS.nc'
+IN_FN="$(basename ${IN_URI})"
+IN_FN_ROOT="${IN_FN%.*}"
+IN_FN_EXT="${IN_FN##*.}"
+export IN_FP="${WORK_DIR}/${IN_FN}"
+export IN_DATAVAR_NAME='zS'
+# TODO: get from file, e.g., with NCL
+export IN_DATAVAR_NA='-2147483648'
+export IN_DATAVAR_BAND='NA'
+IN_DATAVAR_UNIT='m'
+IN_DATAVAR_LONGNAME='surface elevation'
+IN_DATAVAR_COORD_X_NAME='lon'
+IN_DATAVAR_COORD_Y_NAME='lat'
+
+export RASTER_BRICK='FALSE'
+export RASTER_ROTATE='FALSE'
+export RASTER_PLOT='TRUE'  # can't make R use this???
+
+# output: cropped surface elevation from ETOPO1 regridded to match template
+OUT_FN_ROOT="${IN_FN_ROOT}_regrid"
+OUT_FN="${OUT_FN_ROOT}.${IN_FN_EXT}"
+export OUT_FP="${WORK_DIR}/${OUT_FN}"
+export OUT_DATAVAR_NAME="${IN_DATAVAR_NAME}"
+# note type=int in original, but probably should be float for regridding
+# raster::writeRaster seems to autocast down unless instructed otherwise?
+export OUT_DATAVAR_TYPE='FLT4S' # == float: see raster.pdf::dataType
+export OUT_DATAVAR_NA="${IN_DATAVAR_NA}"
+export OUT_DATAVAR_BAND="${IN_DATAVAR_BAND}"
+export OUT_DATAVAR_UNIT="${IN_DATAVAR_UNIT}"
+export OUT_DATAVAR_LONGNAME="${IN_DATAVAR_LONGNAME}"
+export OUT_DATAVAR_COORD_X_NAME="${IN_DATAVAR_COORD_X_NAME}"
+export OUT_DATAVAR_COORD_Y_NAME="${IN_DATAVAR_COORD_Y_NAME}"
+# export OUT_DATAVAR_COORD_Z_NAME="${IN_DATAVAR_COORD_Z_NAME}"
+# export OUT_DATAVAR_COORD_Z_UNIT="${IN_DATAVAR_COORD_Z_UNIT}"
+
+OUT_PDF_FN="${OUT_FN_ROOT}_$(date +${DATE_FORMAT}).pdf"
+export OUT_PDF_FP="${PDF_DIR}/${OUT_PDF_FN}" # path to PDF output
+export OUT_PDF_HEIGHT='15'
+export OUT_PDF_WIDTH='20'
+
+# called scripts. TODO: R-package my code!
+
+# this script does simple statistics, and is copied from ioapi-hack-R (also on github)
+# TODO: R-package my code
+STAT_SCRIPT_URI='https://github.com/TomRoche/GEIA_to_netCDF/raw/master/netCDF.stats.to.stdout.r'
+STAT_SCRIPT_FN="$(basename ${STAT_SCRIPT_URI})"
+export STAT_SCRIPT_FP="${WORK_DIR}/${STAT_SCRIPT_FN}"
+
+# payload---------------------------------------------------------------
+
+# setup-----------------------------------------------------------------
+
+mkdir -p ${WORK_DIR}
+
+# get scripts-----------------------------------------------------------
+
+if [[ ! -r "${STAT_SCRIPT_FP}" ]] ; then
+  for CMD in \
+    "wget --no-check-certificate -c -O ${STAT_SCRIPT_FP} ${STAT_SCRIPT_URI}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+# get data------------------------------------------------------------
+
+# template file
+if [[ ! -r "${TEMPLATE_FP}" ]] ; then
+  for CMD in \
+    "wget --no-check-certificate -c -O ${TEMPLATE_FP} ${TEMPLATE_URI}" \
+    "ncdump -h ${TEMPLATE_FP}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+# elevation source
+if [[ ! -r "${IN_FP}" ]] ; then
+  for CMD in \
+    "wget --no-check-certificate -c -O ${IN_FP} ${IN_URI}" \
+    "ncdump -h ${IN_FP}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+# call R script (the real work)---------------------------------------
+
+# TODO: pass commandline args to Rscript
+# punt: just use envvars :-(
+# for CMD in \
+#   "${RSCRIPT} ${CALL_R_FP}" \
+${RSCRIPT} ${CALL_R_FP}
+#   ; do
+#   echo -e "$ ${CMD}"
+#   eval "${CMD}"
+# done
+
+# After exiting R, show output files and display output PDF.
+for CMD in \
+  "ls -alht ${WORK_DIR} | head" \
+; do
+  echo -e "$ ${CMD}"
+  eval "${CMD}"
+done
+
+if [[ -r "${OUT_PDF_FP}" ]] ; then
+  for CMD in \
+    "${PDF_VIEWER} ${OUT_PDF_FP} &" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+else
+  echo -e "ERROR: ${THIS_FN}: PDF output file='${OUT_PDF_FP}' not readable"
+fi

File restrict.lon.lat.to.region.r

-# Code to restrict global data in a netCDF file to a region/window, save, and plot.
-# TODO: make generic, pass parameters
-
-# constants-----------------------------------------------------------
-
-cat('Setting global constants ...\n') # debugging
-
-# global constants, mostly from environment---------------------------
-
-cat('Loading package=raster ...\n') # debugging
-library(raster)
-# coordinate reference system: here, just global lon-lat
-global.proj <- CRS('+proj=longlat +ellps=WGS84')
-
-# all the following env vars must be set and exported in driver script (TODO: fail if not)
-this.fn <- Sys.getenv('CALL_R_FN')  # FIXME: gets overidden by `source` below
-cat(sprintf('this.fn=%s\n', this.fn)) # debugging
-
-sigdigs <- as.numeric(Sys.getenv('SIGDIGS')) # significant digits for lattice strips
-
-# bounds to which to restrict: 1° outside AQMEII-NA bounds
-# longitudes: -130° to -59.5° -> -131° to -58.5°
-# latitudes: 23.5° to 58.5°   ->  22.5° to 59.5°
-restrict.lat.west.deg <- -131.0  # just west of AQMEII-NA
-restrict.lat.east.deg <- -58.5   # just east of AQMEII-NA
-restrict.lon.south.deg <- 22.5   # just south of AQMEII-NA
-restrict.lon.north.deg <- 59.5   # just north of AQMEII-NA
-
-aqmeii.extent <-
-  extent(restrict.lat.west.deg, restrict.lat.east.deg,
-         restrict.lon.south.deg, restrict.lon.north.deg)
-# aqmeii.extent  # debugging
-# class       : Extent 
-# xmin        : -131 
-# xmax        : -58.5 
-# ymin        : 22.5 
-# ymax        : 59.5 
-
-stat.script.fp <- Sys.getenv('STAT_SCRIPT_FP')
-source(stat.script.fp) # produces errant error=
-#> netCDF.stats.to.stdout.r: no arguments supplied, exiting
-# unnecessary here
-# plot.script.fp <- Sys.getenv('PLOT_SCRIPT_FP')
-# source(plot.script.fp)
-
-# local constants, from environment-----------------------------------
-# taking from Rscript call failed, annoyingly and mysteriously :-(
-
-cat('Setting local constants ...\n') # debugging
-
-raster.brick <- Sys.getenv('RASTER_BRICK')
-raster.rotate <- Sys.getenv('RASTER_ROTATE')
-raster.plot <- Sys.getenv('RASTER_PLOT')
-in.fp <- Sys.getenv('IN_FP')
-in.datavar.name <- Sys.getenv('IN_DATAVAR_NAME')
-in.datavar.type <- Sys.getenv('IN_DATAVAR_TYPE')
-in.datavar.band <- Sys.getenv('IN_DATAVAR_BAND')
-out.fp <- Sys.getenv('OUT_FP')
-out.datavar.name <- Sys.getenv('OUT_DATAVAR_NAME')
-out.datavar.na <- Sys.getenv('OUT_DATAVAR_NA')
-out.datavar.unit <- Sys.getenv('OUT_DATAVAR_UNIT')
-out.datavar.longname <- Sys.getenv('OUT_DATAVAR_LONGNAME')
-out.datavar.coord.x.name <- Sys.getenv('OUT_DATAVAR_COORD_X_NAME')
-out.datavar.coord.y.name <- Sys.getenv('OUT_DATAVAR_COORD_Y_NAME')
-out.datavar.coord.z.name <- Sys.getenv('OUT_DATAVAR_COORD_Z_NAME')
-out.datavar.coord.z.unit <- Sys.getenv('OUT_DATAVAR_COORD_Z_UNIT')
-pdf.fp <- Sys.getenv('PDF_FP')
-pdf.height <- Sys.getenv('PDF_HEIGHT')
-pdf.width <- Sys.getenv('PDF_WIDTH')
-
-# start debugging
-cat(sprintf('raster.brick==%s\n', raster.brick))
-cat(sprintf('raster.rotate==%s\n', raster.rotate))
-cat(sprintf('raster.plot==%s\n', raster.plot))
-cat(sprintf('in.fp==%s\n', in.fp))
-cat(sprintf('in.datavar.name==%s\n', in.datavar.name))
-cat(sprintf('in.datavar.type==%s\n', in.datavar.type))
-cat(sprintf('in.datavar.band==%s\n', in.datavar.band))
-cat(sprintf('out.fp==%s\n', out.fp))
-cat(sprintf('out.datavar.name==%s\n', out.datavar.name))
-cat(sprintf('out.datavar.na==%s\n', out.datavar.na))
-cat(sprintf('out.datavar.unit==%s\n', out.datavar.unit))
-cat(sprintf('out.datavar.longname==%s\n', out.datavar.longname))
-cat(sprintf('out.datavar.coord.x.name==%s\n', out.datavar.coord.x.name))
-cat(sprintf('out.datavar.coord.y.name==%s\n', out.datavar.coord.y.name))
-cat(sprintf('out.datavar.coord.z.name==%s\n', out.datavar.coord.z.name))
-cat(sprintf('out.datavar.coord.z.unit==%s\n', out.datavar.coord.z.unit))
-cat(sprintf('pdf.fp==%s\n', pdf.fp))
-cat(sprintf('pdf.height==%s\n', pdf.height))
-cat(sprintf('pdf.width==%s\n', pdf.width))
-#   end debugging
-
-raster.brick <- as.logical(raster.brick)
-raster.rotate <- as.logical(raster.rotate)
-raster.plot <- as.logical(raster.plot)
-out.datavar.na <- as.numeric(out.datavar.na)
-pdf.height <- as.numeric(pdf.height)
-pdf.width <- as.numeric(pdf.width)
-
-# ensure global scope, so raster.aqmeii can be accessed inside raster.plot conditional? nope, fail
-# raster.crop <- TRUE # kludge
-
-# functions-----------------------------------------------------------
-
-# fails to fix the conditional scoping problem below
-## plot <- function(to.plot, sigdigs, raster.to.plot) { # plot(to.plot, raster.to.plot)
-##   if (!to.plot) return
-##   # 5.1 Fix layer names (which are plot labels in `layerplot`)
-
-##   # TODO: since using rasterVis::levelplot, get/set names with `layerNames`
-##   raster.to.plot.df <- as.data.frame(raster.to.plot)
-##   names(raster.to.plot.df) <- formatC(
-##     as.numeric(sub(x=names(raster.to.plot.df), pattern="X", replacement="", fixed=TRUE)),
-##     format="g", digits=sigdigs)
-##   # start debugging
-##   cat(sprintf('layers.n==%s\n', layers.n))
-##   ## raster.to.plot
-##   ## raster.to.plot.df
-##   ## cat('head and tail of raster.to.plot.df:\n')
-##   ## head(raster.to.plot.df)
-##   ## tail(raster.to.plot.df)
-##   cat('head and tail of names(raster.to.plot.df):\n')
-##   head(names(raster.to.plot.df))
-##   tail(names(raster.to.plot.df))
-##   #  end debugging
-
-##   # 5.2 Create world map (which rasterVis will autocrop)
-
-##   cat('Loading map packages ...\n') # debugging
-##   library(maps)
-##   library(maptools) # needed for map2SpatialLines
-##   data(wrld_simpl)      # from maptools
-
-##   # TODO: save/restore after first creation!
-##   global.map <- map('world', plot=FALSE)
-##   global.map.shp <- map2SpatialLines(global.map, proj4string=global.proj)
-
-##   # 5.3 Plot and display.
-
-##   cat('Loading plot packages ...\n') # debugging
-##   library(rasterVis)
-
-##   cat(sprintf(
-##     '%s: plotting %s may take awhile! (TODO: add progress control)\n',
-##     this.fn, pdf.fp))
-##   pdf(file=pdf.fp, width=pdf.width, height=pdf.height)
-##   rasterVis::levelplot(raster.to.plot,
-##     margin=FALSE,
-##     names.attr=names(raster.to.plot.df),
-##     layout=c(1,layers.n), # all layers in one "atmospheric" column
-##   ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
-##   # teardown
-##   dev.off()
-##   # system(sprintf('xpdf %s &', pdf.fp)) # do in driver
-##   # the plot crops the map!
-## } # end function plot(to.plot, raster.to.plot)
-
-# code----------------------------------------------------------------
-
-# 1. Load the input global datavar
-
-if (raster.brick) { # multiple vertical layers
-  raster.global <- brick(in.fp, varname=in.datavar.name, crs=global.proj)
-  layers.n <- nbands(raster.global)
-} else {
-  raster.global <- raster(in.fp, varname=in.datavar.name, crs=global.proj)
-  layers.n <- 1
-}
-# raster.global # debugging
-# e.g.
-# class       : RasterBrick 
-# dimensions  : 96, 144, 13824, 56  (nrow, ncol, ncell, nlayers)
-# resolution  : 2.5, 1.894737  (x, y)
-# extent      : -1.25, 358.75, -90.94737, 90.94737  (xmin, xmax, ymin, ymax)
-# coord. ref. : +proj=longlat +datum=WGS84 
-# data source : ~/code/regridding/MOZART_global_to_AQMEII-NA/2008N2O_restart.nc 
-# names       : X1.86787999700755, X2.35259043984115, X2.94832093641162, X3.67650110274553, X4.56168595701456, X5.63180120661855, X6.91832136362791, X8.4563922137022, X10.2849211543798, X12.4601498246193, X15.0502501055598, X18.1243494153023, X21.761005744338, X26.0491091758013, X31.0889091342688, ... 
-# hybrid_sigma_pressure: 1.86787999700755, 992.500010610456 (min, max)
-# varname     : N2O 
-
-# 2. Rotate the input global datavar if requested:
-#    zero-based longitudes cause problems for wrld_simpl (and probably other data)
-if (raster.rotate) {
-  raster.global <- rotate(raster.global,
-    overwrite=TRUE) # else levelplot does one layer per page?
-#   raster.global # debugging
-# e.g.
-# class       : RasterBrick 
-# dimensions  : 96, 144, 13824, 56  (nrow, ncol, ncell, nlayers)
-# resolution  : 2.5, 1.894737  (x, y)
-# extent      : -181.25, 178.75, -90.94737, 90.94737  (xmin, xmax, ymin, ymax)
-# coord. ref. : +proj=longlat +datum=WGS84 
-# data source : in memory
-# ...
-}
-
-# 3. Crop datavar to our region.
-# snap='out' -> output includes all cells intersected by the extent (whether or not the extent covers the cell center)
-
-# ensure global scope, so raster.aqmeii can be accessed inside raster.plot conditional? nope, fail
-# if (raster.crop) {
-  raster.aqmeii <- crop(x=raster.global, y=aqmeii.extent, snap='out')
-# ensure global scope, so raster.aqmeii can be accessed inside raster.plot conditional? nope, fail
-# raster.aqmeii <<- crop(x=raster.global, y=aqmeii.extent, snap='out') 
-
-# raster.aqmeii # debugging
-# e.g.
-# class       : RasterBrick 
-# dimensions  : 21, 30, 630, 56  (nrow, ncol, ncell, nlayers)
-# resolution  : 2.5, 1.894737  (x, y)
-# extent      : -131.25, -56.25, 20.84211, 60.63158  (xmin, xmax, ymin, ymax)
-# coord. ref. : +proj=longlat +datum=WGS84 
-# data source : in memory
-# ...
-# } # if (raster.crop)
-
-# 4. Write cropped datavar to output netCDF.
-# save to file, overwriting pre-existing
-cat(sprintf('Writing cropped raster to path=%s ...\n', out.fp)) # debugging
-if (raster.brick) {
-  # ... we have layers, i.e., z
-  writeRaster(
-    x=raster.aqmeii, filename=out.fp, varname=out.datavar.name,
-    format="CDF", NAflag=out.datavar.na, overwrite=TRUE,
-    # datavar type will truncate unless set?
-    datatype=in.datavar.type,
-    # netCDF-specific arguments
-    varunit=out.datavar.unit, longname=out.datavar.longname,
-    xname=out.datavar.coord.x.name, yname=out.datavar.coord.y.name,
-    zname=out.datavar.coord.z.name, zunit=out.datavar.coord.z.unit
-  )
-} else { # it's a RasterLayer
-  writeRaster(
-    x=raster.aqmeii, filename=out.fp, varname=out.datavar.name,
-    format="CDF", NAflag=out.datavar.na, overwrite=TRUE,
-    # datavar type will truncate unless set?
-    datatype=in.datavar.type,
-    # netCDF-specific arguments
-    varunit=out.datavar.unit, longname=out.datavar.longname,
-    xname=out.datavar.coord.x.name, yname=out.datavar.coord.y.name
-  )
-}
-
-# start debugging
-system(sprintf('ls -alh %s', out.fp))
-netCDF.stats.to.stdout(netcdf.fp=in.fp, var.name=in.datavar.name)
-netCDF.stats.to.stdout(netcdf.fp=out.fp, var.name=out.datavar.name)
-system(sprintf('ncdump -h %s', out.fp))
-# e.g.
-# For /home/rtd/code/regridding/MOZART_global_to_AQMEII-NA/2008N2O_restart_region_N2O.nc var=N2O
-#         cells=35280
-#         obs=35280
-#         min=0.0151
-#         q1=0.242
-#         med=0.485
-#         mean=0.367
-#         q3=0.486
-#         max=0.492
-
-# netcdf \2008N2O_restart_region_N2O {
-# dimensions:
-#         lon = 30 ;
-#         lat = 21 ;
-#         lev = UNLIMITED ; // (56 currently)
-# variables:
-#         double lon(lon) ;
-#                 lon:units = "degrees_east" ;
-#                 lon:long_name = "lon" ;
-#         double lat(lat) ;
-#                 lat:units = "degrees_north" ;
-#                 lat:long_name = "lat" ;
-#         int lev(lev) ;
-#                 lev:units = "hybrid_sigma_pressure" ;
-#                 lev:long_name = "lev" ;
-#         double N2O(lev, lat, lon) ;
-#                 N2O:units = "ppmV" ;
-#                 N2O:_FillValue = -1. ;
-#                 N2O:missing_value = -1. ;
-#                 N2O:long_name = "N2O" ;
-#                 N2O:projection = "+proj=longlat +datum=WGS84" ;
-#                 N2O:projection_format = "PROJ.4" ;
-#   end debugging
-
-# 5. Plot data if requested
-
-# 5.0 Massive annoyance: raster.aqmeii is out of scope inside this conditional!
-# Despite numerous attempted kludges, I cannot make plotting conditional!
-# if (raster.plot) {
-
-# 5.1 Fix layer names (which are plot labels in `layerplot`)
-
-if (raster.brick) {
-  # TODO: since using rasterVis::levelplot, get/set names with `layerNames`
-  aqmeii.df <- as.data.frame(raster.aqmeii)
-  names(aqmeii.df) <- formatC(
-    as.numeric(sub(x=names(aqmeii.df), pattern="X", replacement="", fixed=TRUE)),
-    format="g", digits=sigdigs)
-  # start debugging
-  cat(sprintf('layers.n==%s\n', layers.n))
-  ## raster.aqmeii
-  ## aqmeii.df
-  ## cat('head and tail of aqmeii.df:\n')
-  ## head(aqmeii.df)
-  ## tail(aqmeii.df)
-  cat('head and tail of names(aqmeii.df):\n')
-  head(names(aqmeii.df))
-  tail(names(aqmeii.df))
-  #  end debugging
-} # end if (raster.brick)
-
-# 5.2 Create world map (which rasterVis will autocrop)
-
-  cat('Loading map packages ...\n') # debugging
-  library(maps)
-  library(maptools) # needed for map2SpatialLines
-  data(wrld_simpl)      # from maptools
-
-  # TODO: save/restore after first creation!
-  global.map <- map('world', plot=FALSE)
-  global.map.shp <- map2SpatialLines(global.map, proj4string=global.proj)
-
-# 5.3 Plot and display.
-
-  cat('Loading plot packages ...\n') # debugging
-  library(rasterVis)
-
-  cat(sprintf(
-    '%s: plotting %s may take awhile! (TODO: add progress control)\n',
-    this.fn, pdf.fp))
-#  if (raster.plot) {
-  pdf(file=pdf.fp, width=pdf.width, height=pdf.height)
-  if (raster.brick) {
-    rasterVis::levelplot(raster.aqmeii,
-      margin=FALSE,
-      names.attr=names(aqmeii.df),
-      layout=c(1,layers.n), # all layers in one "atmospheric" column
-    ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
-  } else {
-    rasterVis::levelplot(raster.aqmeii,
-      margin=FALSE
-    ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
-    # the plot crops the map!
-  } # end if (raster.brick)
-  # teardown
-  dev.off()
-
-# } # end if (raster.plot)

File restrict_lon_lat_to_region.sh

-#!/usr/bin/env bash
-# Requires R.
-# Driver for restrict.lon.lat.to.region.r--configure as needed for your platform.
-# Crops the following netCDF {file, datavar} pairs, saving output to separate file:
-# * 2008N2O_restart.nc,N2O
-# * 2008N2O_restart.nc,PS
-# * ETOPO1_Ice_g_gmt4.grd.nc,zS: after renaming
-# ** x->lon (to match 2008N2O_restart.nc)
-# ** y->lat (to match 2008N2O_restart.nc)
-# ** z->zS  (à la 2008N2O_restart.nc:PS)
-
-# constants with some simple manipulations------------------------------
-
-THIS_FN="$(basename $0)"
-# will do the real work
-# TODO: create from THIS_FN, e.g., use `sed -e s/_/./g`
-export CALL_R_FN='restrict.lon.lat.to.region.r'
-
-# modules on EMVL. TODO: fixme! (these fail from here, are currently just reminders.)
-# You may not need them at all!
-module add netcdf-4.1.2
-#module add ncl_ncarg-6.0.0
-#module add nco-4.0.5
-
-# for R on EMVL: don't ask :-(
-R_DIR='/usr/local/R-2.15.0/bin'
-R="${R_DIR}/R"
-RSCRIPT="${R_DIR}/Rscript"
-
-# workspace
-export TEST_DIR="$(pwd)" # keep it simple for now: same dir as top of repo
-mkdir -p ${TEST_DIR}
-
-# will do the real work (restricting and plotting)
-export CALL_R_FP="${TEST_DIR}/${CALL_R_FN}"
-
-# for data display
-export SIGDIGS='4' # significant digits
-
-# for plotting
-PDF_VIEWER='xpdf'  # whatever works on your platform
-# temporally disaggregate multiple plots
-DATE_FORMAT='%Y%m%d_%H%M'
-export PDF_DIR="${TEST_DIR}"
-
-# global lon-lat netCDFs: get from my repository, or download original if available
-
-# for N2O, PS
-# TODO: download/process original
-SAIKAWA_URI='http://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart.nc'
-SAIKAWA_FN="$(basename ${SAIKAWA_URI})"
-# get both parts of the filename (before and after '.')
-SAIKAWA_FN_ROOT="${SAIKAWA_FN%.*}"
-SAIKAWA_FN_EXT="${SAIKAWA_FN##*.}"
-export SAIKAWA_FP="${TEST_DIR}/${SAIKAWA_FN}"
-
-# for surface elevation
-ETOPO1_GZ_URI='http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/netcdf/ETOPO1_Ice_g_gmt4.grd.gz'
-ETOPO1_FN=$(echo $(basename ${ETOPO1_GZ_URI%%\?*}) | sed -e 's/gz$/nc/') # ensure standard suffix
-ETOPO1_FN_ROOT="${ETOPO1_FN%.*}"
-ETOPO1_FN_EXT="${ETOPO1_FN##*.}"
-export ETOPO1_FP="${TEST_DIR}/${ETOPO1_FN}"
-
-# metadata for datavars and their attributes. TODO: read from netCDF
-
-# the N2O concentrations: N2O(lev, lat, lon)
-N2O_DATAVAR_NAME='N2O'
-# note type==double in original
-N2O_DATAVAR_TYPE='FLT8S' # == double: see raster.pdf::dataType
-N2O_DATAVAR_UNIT='ppmV'  # but only after pre-processing!
-N2O_DATAVAR_NA='-1.0'    # assigned here, not in netCDF
-N2O_DATAVAR_BAND='NA'    # input has no timestep defined: N2O(lev, lat, lon)
-N2O_REGIONAL_FN_ROOT="${SAIKAWA_FN_ROOT}_region_${N2O_DATAVAR_NAME}" # used below
-N2O_REGIONAL_FN="${N2O_REGIONAL_FN_ROOT}.${SAIKAWA_FN_EXT}"
-N2O_REGIONAL_FP="${TEST_DIR}/${N2O_REGIONAL_FN}"
-N2O_PDF_FN="${N2O_REGIONAL_FN_ROOT}_$(date +${DATE_FORMAT}).pdf"
-N2O_PDF_FP="${PDF_DIR}/${N2O_PDF_FN}" # path to PDF output
-N2O_PDF_HEIGHT='5'
-N2O_PDF_WIDTH='100'
-
-# following datavars required for vertical coordinate conversion (later): hybrid sigma-pressure -> elevation
-
-# double PS(lat, lon)
-SURF_PRESSURE_DATAVAR_NAME='PS'
-# note type==double in original
-SURF_PRESSURE_DATAVAR_TYPE='FLT8S' # == double: see raster.pdf::dataType
-SURF_PRESSURE_DATAVAR_UNIT='Pa'
-SURF_PRESSURE_DATAVAR_LONGNAME='Surface Pressure' # original case
-SURF_PRESSURE_DATAVAR_NA='-1.0' # assigned here, not in netCDF
-SURF_PRESSURE_DATAVAR_BAND='NA' # input has no timestep defined
-SURF_PRESSURE_DATAVAR_COORD_X_NAME='lon'
-SURF_PRESSURE_DATAVAR_COORD_Y_NAME='lat'
-SURF_PRESSURE_GLOBAL_FP="${SAIKAWA_FP}"
-SURF_PRESSURE_REGIONAL_FN_ROOT="${SAIKAWA_FN_ROOT}_region_${SURF_PRESSURE_DATAVAR_NAME}" # used below
-SURF_PRESSURE_REGIONAL_FN="${SURF_PRESSURE_REGIONAL_FN_ROOT}.${SAIKAWA_FN_EXT}"
-SURF_PRESSURE_REGIONAL_FP="${TEST_DIR}/${SURF_PRESSURE_REGIONAL_FN}"
-SURF_PRESSURE_PDF_FN="${SURF_PRESSURE_REGIONAL_FN_ROOT}_$(date +${DATE_FORMAT}).pdf"
-SURF_PRESSURE_PDF_FP="${PDF_DIR}/${SURF_PRESSURE_PDF_FN}" # path to PDF output
-
-# int z(y, x) -> double zS(lat, lon)
-SURF_ELEVATION_DATAVAR_NAME_IN='z'
-SURF_ELEVATION_DATAVAR_NAME_OUT='zS'
-# note type=int in original, but probably should be float for regridding
-SURF_ELEVATION_DATAVAR_TYPE='FLT4S' # == float: see raster.pdf::dataType
-SURF_ELEVATION_DATAVAR_UNIT='m'
-SURF_ELEVATION_DATAVAR_LONGNAME='surface elevation' # assigned here
-SURF_ELEVATION_DATAVAR_NA='-2147483648' # use value from original
-SURF_ELEVATION_DATAVAR_BAND='NA' # input has no timestep defined
-SURF_ELEVATION_DATAVAR_COORD_X_NAME_IN='x'
-SURF_ELEVATION_DATAVAR_COORD_X_NAME_OUT='lon'
-SURF_ELEVATION_DATAVAR_COORD_Y_NAME_IN='y'
-SURF_ELEVATION_DATAVAR_COORD_Y_NAME_OUT='lat'
-SURF_ELEVATION_GLOBAL_FP="${ETOPO1_FP}"
-SURF_ELEVATION_REGIONAL_FN_ROOT="${ETOPO1_FN_ROOT}_region_${SURF_ELEVATION_DATAVAR_NAME_OUT}" # used below
-SURF_ELEVATION_REGIONAL_FN="${SURF_ELEVATION_REGIONAL_FN_ROOT}.${ETOPO1_FN_EXT}"
-SURF_ELEVATION_REGIONAL_FP="${TEST_DIR}/${SURF_ELEVATION_REGIONAL_FN}"
-SURF_ELEVATION_PDF_FN="${SURF_ELEVATION_REGIONAL_FN_ROOT}_$(date +${DATE_FORMAT}).pdf"
-SURF_ELEVATION_PDF_FP="${PDF_DIR}/${SURF_ELEVATION_PDF_FN}" # path to PDF output
-
-# called scripts. TODO: R-package my code!
-
-# # this script drives image.plot, and is copied from ioapi-hack-R (also on github)
-# PLOT_SCRIPT_URI='https://github.com/TomRoche/GEIA_to_netCDF/raw/master/plotLayersForTimestep.r'
-# PLOT_SCRIPT_FN="$(basename ${PLOT_SCRIPT_URI})"
-# export PLOT_SCRIPT_FP="${TEST_DIR}/${PLOT_SCRIPT_FN}"
-
-# this script does simple statistics, and is copied from ioapi-hack-R (also on github)
-# TODO: R-package my code
-STAT_SCRIPT_URI='https://github.com/TomRoche/GEIA_to_netCDF/raw/master/netCDF.stats.to.stdout.r'
-STAT_SCRIPT_FN="$(basename ${STAT_SCRIPT_URI})"
-export STAT_SCRIPT_FP="${TEST_DIR}/${STAT_SCRIPT_FN}"
-
-# payload---------------------------------------------------------------
-
-# get scripts-----------------------------------------------------------
-
-# if [[ ! -r "${PLOT_SCRIPT_FP}" ]] ; then
-#   for CMD in \
-#     "wget --no-check-certificate -c -O ${PLOT_SCRIPT_FP} ${PLOT_SCRIPT_URI}" \
-#   ; do
-#     echo -e "$ ${CMD}"
-#     eval "${CMD}"
-#   done
-# fi
-
-if [[ ! -r "${STAT_SCRIPT_FP}" ]] ; then
-  for CMD in \
-    "wget --no-check-certificate -c -O ${STAT_SCRIPT_FP} ${STAT_SCRIPT_URI}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# get data------------------------------------------------------------
-
-# N2O, PS source
-if [[ ! -r "${SAIKAWA_FP}" ]] ; then
-  for CMD in \
-    "wget --no-check-certificate -c -O ${SAIKAWA_FP} ${SAIKAWA_URI}" \
-    "ncdump -h ${SAIKAWA_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# elevation source: is gzip'ed, must be unpacked first
-if [[ ! -r "${ETOPO1_FP}" ]] ; then
-  for CMD in \
-    "curl -C - ${ETOPO1_GZ_URI} | gunzip -c > ${ETOPO1_FP}" \
-    "ncdump -h ${ETOPO1_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# call R script (the real work)---------------------------------------
-
-echo # newline
-
-# # for each of our 3 tuples={filepath_in, datavar_in, filepath_out, datavar_out} 
-# for PAIR in \
-#   "${N2O_GLOBAL_FP},${N2O_DATAVAR_NAME_IN},${N2O_REGIONAL_FP},${N2O_DATAVAR_NAME_OUT}" \
-#   "${SURF_PRESSURE_GLOBAL_FP},${SURF_PRESSURE_DATAVAR_NAME_IN},${SURF_PRESSURE_REGIONAL_FP},${SURF_PRESSURE_DATAVAR_NAME_OUT}" \
-#   "${SURF_ELEVATION_GLOBAL_FP},${SURF_ELEVATION_DATAVAR_NAME_IN},${SURF_ELEVATION_REGIONAL_FP},${SURF_ELEVATION_DATAVAR_NAME_OUT}" \
-# ; do
-#    "${RSCRIPT} ${CALL_R_FP} in.fp='' in.datavar.name='' out.fp='' out.datavar.name=''" \
-# done
-
-# TODO: make *1* R invocation with all that package loading!
-
-# # Note R-script argument quoting is *very annoying*!
-# # in.fp='${SAIKAWA_FP}' \
-# # in.fp=\"${SAIKAWA_FP}\" \
-# # in.fp=${SAIKAWA_FP} \
-# for CMD in \
-#   "${RSCRIPT} ${CALL_R_FP} \
-# 'in.fp=\"${SAIKAWA_FP}\"' \
-# # fails with
-# # > Error in eval.with.vis(expr, envir, enclos) : unknown argument='in.fp'
-# in.datavar.name=\"N2O\" \
-# in.datavar.band=\"NA\" \
-# out.fp=\"${N2O_REGIONAL_FP}\" \
-# out.datavar.name=\"N2O\" \
-# raster.brick=\"TRUE\" \
-# raster.rotate=\"TRUE\" \
-# raster.plot=\"TRUE\" \
-# out.datavar.na=\"${N2O_DATAVAR_NA}\" \
-# out.datavar.unit=\"${N2O_DATAVAR_UNIT}\" \
-# out.datavar.longname=\"N2O concentration\" \
-# out.datavar.coord.x.name=\"lon\" \
-# out.datavar.coord.y.name=\"lat\" \
-# out.datavar.coord.z.name=\"lev\" \
-# out.datavar.coord.z.unit=\"hybrid_sigma_pressure\"" \
-#   ; do
-#   echo -e "$ ${CMD}"
-#   eval "${CMD}"
-# done
-
-# punt: just use envvars :-(
-
-# take 1: N2O
-export IN_FP="${SAIKAWA_FP}"
-export IN_DATAVAR_NAME="${N2O_DATAVAR_NAME}"
-export IN_DATAVAR_BAND='NA'
-export IN_DATAVAR_TYPE="${N2O_DATAVAR_TYPE}"
-export OUT_FP="${N2O_REGIONAL_FP}"
-export OUT_DATAVAR_NAME='N2O'
-export RASTER_BRICK='TRUE'
-export RASTER_ROTATE='TRUE'
-export RASTER_PLOT='TRUE'  # can't make R use this???
-export OUT_DATAVAR_NA="${N2O_DATAVAR_NA}"
-export OUT_DATAVAR_UNIT="${N2O_DATAVAR_UNIT}"
-export OUT_DATAVAR_LONGNAME='N2O concentration'
-export OUT_DATAVAR_COORD_X_NAME='lon'
-export OUT_DATAVAR_COORD_Y_NAME='lat'
-export OUT_DATAVAR_COORD_Z_NAME='lev'
-export OUT_DATAVAR_COORD_Z_UNIT='hybrid_sigma_pressure'
-export PDF_FP="${N2O_PDF_FP}" # path to PDF output
-export PDF_HEIGHT='100'
-export PDF_WIDTH='5'
-# for CMD in \
-#   "${RSCRIPT} ${CALL_R_FP}" \
-${RSCRIPT} ${CALL_R_FP}
-#   ; do
-#   echo -e "$ ${CMD}"
-#   eval "${CMD}"
-# done
-
-# After exiting R, show output files and display output PDF.
-for CMD in \
-  "ls -alht ${TEST_DIR}" \
-; do
-  echo -e "$ ${CMD}"
-  eval "${CMD}"
-done
-
-if [[ -r "${PDF_FP}" ]] ; then
-  for CMD in \
-    "${PDF_VIEWER} ${PDF_FP} &" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-else
-  echo -e "ERROR: ${THIS_FN}: PDF output file='${PDF_FP}' not readable"
-fi
-
-# take 2: surface pressure
-export IN_FP="${SAIKAWA_FP}"
-export IN_DATAVAR_NAME="${SURF_PRESSURE_DATAVAR_NAME}"
-export IN_DATAVAR_BAND='NA'
-export IN_DATAVAR_TYPE="${SURF_PRESSURE_DATAVAR_TYPE}"
-export OUT_FP="${SURF_PRESSURE_REGIONAL_FP}"
-export OUT_DATAVAR_NAME="${SURF_PRESSURE_DATAVAR_NAME}"
-export RASTER_BRICK='FALSE' # only one layer
-export RASTER_ROTATE='TRUE'
-export RASTER_PLOT='TRUE'  # can't make R use this???
-export OUT_DATAVAR_NA="${SURF_PRESSURE_DATAVAR_NA}"
-export OUT_DATAVAR_UNIT="${SURF_PRESSURE_DATAVAR_UNIT}"
-export OUT_DATAVAR_LONGNAME="${SURF_PRESSURE_DATAVAR_LONGNAME}"
-export OUT_DATAVAR_COORD_X_NAME='lon'
-export OUT_DATAVAR_COORD_Y_NAME='lat'
-export OUT_DATAVAR_COORD_Z_NAME=''
-export OUT_DATAVAR_COORD_Z_UNIT=''
-export PDF_FP="${SURF_PRESSURE_PDF_FP}" # path to PDF output
-export PDF_HEIGHT='15'
-export PDF_WIDTH='20'
-# for CMD in \
-#   "${RSCRIPT} ${CALL_R_FP}" \
-${RSCRIPT} ${CALL_R_FP}
-#   ; do
-#   echo -e "$ ${CMD}"
-#   eval "${CMD}"
-# done
-
-# After exiting R, show output files and display output PDF.
-for CMD in \
-  "ls -alht ${TEST_DIR}" \
-; do
-  echo -e "$ ${CMD}"
-  eval "${CMD}"
-done
-
-if [[ -r "${PDF_FP}" ]] ; then
-  for CMD in \
-    "${PDF_VIEWER} ${PDF_FP} &" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-else
-  echo -e "ERROR: ${THIS_FN}: PDF output file='${PDF_FP}' not readable"
-fi
-
-# take 3: surface elevation
-export IN_FP="${ETOPO1_FP}"
-export IN_DATAVAR_NAME="${SURF_ELEVATION_DATAVAR_NAME_IN}"
-export IN_DATAVAR_BAND='NA'
-export IN_DATAVAR_TYPE="${SURF_ELEVATION_DATAVAR_TYPE}"
-export OUT_FP="${SURF_ELEVATION_REGIONAL_FP}"
-export OUT_DATAVAR_NAME="${SURF_ELEVATION_DATAVAR_NAME_OUT}"
-export RASTER_BRICK='FALSE' # only one layer
-export RASTER_ROTATE='FALSE'
-export RASTER_PLOT='TRUE'  # can't make R use this???
-export OUT_DATAVAR_NA="${SURF_ELEVATION_DATAVAR_NA}"
-export OUT_DATAVAR_UNIT="${SURF_ELEVATION_DATAVAR_UNIT}"
-export OUT_DATAVAR_LONGNAME="${SURF_ELEVATION_DATAVAR_LONGNAME}"
-export OUT_DATAVAR_COORD_X_NAME='lon'
-export OUT_DATAVAR_COORD_Y_NAME='lat'
-# export OUT_DATAVAR_COORD_Z_NAME=''
-# export OUT_DATAVAR_COORD_Z_UNIT=''
-export PDF_FP="${SURF_ELEVATION_PDF_FP}" # path to PDF output
-export PDF_HEIGHT='15'
-export PDF_WIDTH='20'
-# for CMD in \
-#   "${RSCRIPT} ${CALL_R_FP}" \
-${RSCRIPT} ${CALL_R_FP}
-#   ; do
-#   echo -e "$ ${CMD}"
-#   eval "${CMD}"
-# done
-
-# # int z(y, x) -> double zS(lat, lon)
-# SURF_ELEVATION_DATAVAR_NAME_IN='z'
-# SURF_ELEVATION_DATAVAR_NAME_OUT='zS'
-# SURF_ELEVATION_DATAVAR_LONGNAME='surface elevation' # assigned here
-# SURF_ELEVATION_DATAVAR_UNIT='m'
-# # SURF_ELEVATION_DATAVAR_NA='-1.0' # use value from original
-# SURF_ELEVATION_DATAVAR_BAND='NA' # input has no timestep defined
-# SURF_ELEVATION_DATAVAR_COORD_X_NAME_IN='x'
-# SURF_ELEVATION_DATAVAR_COORD_X_NAME_OUT='lon'
-# SURF_ELEVATION_DATAVAR_COORD_Y_NAME_IN='y'
-# SURF_ELEVATION_DATAVAR_COORD_Y_NAME_OUT='lat'
-# # SURF_ELEVATION_GLOBAL_FP="${ETOPO1_FP}"
-# SURF_ELEVATION_REGIONAL_FN_ROOT="${ETOPO1_FN_ROOT}_region_${SURF_ELEVATION_DATAVAR_NAME_OUT}" # used below
-# SURF_ELEVATION_REGIONAL_FN="${SURF_ELEVATION_REGIONAL_FN_ROOT}.${ETOPO1_FN_EXT}"
-# SURF_ELEVATION_REGIONAL_FP="${TEST_DIR}/${SURF_ELEVATION_REGIONAL_FN}"
-# SURF_ELEVATION_PDF_FN="${SURF_ELEVATION_REGIONAL_FN_ROOT}_$(date +${DATE_FORMAT}).pdf"
-# SURF_ELEVATION_PDF_FP="${PDF_DIR}/${SURF_ELEVATION_PDF_FN}" # path to PDF output
-
-# After exiting R, show output files and display output PDF.
-for CMD in \
-  "ls -alht ${TEST_DIR} | head" \
-; do
-  echo -e "$ ${CMD}"
-  eval "${CMD}"
-done
-
-if [[ -r "${PDF_FP}" ]] ; then
-  for CMD in \
-    "${PDF_VIEWER} ${PDF_FP} &" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-else
-  echo -e "ERROR: ${THIS_FN}: PDF output file='${PDF_FP}' not readable"
-fi

File restrict_lon_lat_to_region__MOZART.r

+# Code to restrict global data in a netCDF file to a region/window, save, and plot.
+# TODO: make generic, pass parameters
+
+# constants-----------------------------------------------------------
+
+cat('Setting global constants ...\n') # debugging
+
+# global constants, mostly from environment---------------------------
+
+cat('Loading package=raster ...\n') # debugging
+library(raster)
+# coordinate reference system: here, just global lon-lat
+global.proj <- CRS('+proj=longlat +ellps=WGS84')
+
+# all the following env vars must be set and exported in driver script (TODO: fail if not)
+this.fn <- Sys.getenv('CALL_R_FN')  # FIXME: gets overidden by `source` below
+cat(sprintf('this.fn=%s\n', this.fn)) # debugging
+
+sigdigs <- as.numeric(Sys.getenv('SIGDIGS')) # significant digits for lattice strips
+
+# bounds to which to restrict: 1° outside AQMEII-NA bounds
+# longitudes: -130° to -59.5° -> -131° to -58.5°
+# latitudes: 23.5° to 58.5°   ->  22.5° to 59.5°
+restrict.lat.west.deg <- -131.0  # just west of AQMEII-NA
+restrict.lat.east.deg <- -58.5   # just east of AQMEII-NA
+restrict.lon.south.deg <- 22.5   # just south of AQMEII-NA
+restrict.lon.north.deg <- 59.5   # just north of AQMEII-NA
+
+aqmeii.extent <-
+  extent(restrict.lat.west.deg, restrict.lat.east.deg,
+         restrict.lon.south.deg, restrict.lon.north.deg)
+# aqmeii.extent  # debugging
+# class       : Extent 
+# xmin        : -131 
+# xmax        : -58.5 
+# ymin        : 22.5 
+# ymax        : 59.5 
+
+stat.script.fp <- Sys.getenv('STAT_SCRIPT_FP')
+source(stat.script.fp) # produces errant error=
+#> netCDF.stats.to.stdout.r: no arguments supplied, exiting
+# unnecessary here
+# plot.script.fp <- Sys.getenv('PLOT_SCRIPT_FP')
+# source(plot.script.fp)
+
+# local constants, from environment-----------------------------------
+# taking from Rscript call failed, annoyingly and mysteriously :-(
+
+cat('Setting local constants ...\n') # debugging
+
+raster.brick <- Sys.getenv('RASTER_BRICK')
+raster.rotate <- Sys.getenv('RASTER_ROTATE')
+raster.plot <- Sys.getenv('RASTER_PLOT')
+in.fp <- Sys.getenv('IN_FP')
+in.datavar.name <- Sys.getenv('IN_DATAVAR_NAME')
+in.datavar.type <- Sys.getenv('IN_DATAVAR_TYPE')
+in.datavar.band <- Sys.getenv('IN_DATAVAR_BAND')
+out.fp <- Sys.getenv('OUT_FP')
+out.datavar.name <- Sys.getenv('OUT_DATAVAR_NAME')
+out.datavar.na <- Sys.getenv('OUT_DATAVAR_NA')
+out.datavar.unit <- Sys.getenv('OUT_DATAVAR_UNIT')
+out.datavar.longname <- Sys.getenv('OUT_DATAVAR_LONGNAME')
+out.datavar.coord.x.name <- Sys.getenv('OUT_DATAVAR_COORD_X_NAME')
+out.datavar.coord.y.name <- Sys.getenv('OUT_DATAVAR_COORD_Y_NAME')
+out.datavar.coord.z.name <- Sys.getenv('OUT_DATAVAR_COORD_Z_NAME')
+out.datavar.coord.z.unit <- Sys.getenv('OUT_DATAVAR_COORD_Z_UNIT')
+pdf.fp <- Sys.getenv('PDF_FP')
+pdf.height <- Sys.getenv('PDF_HEIGHT')
+pdf.width <- Sys.getenv('PDF_WIDTH')
+
+# start debugging
+cat(sprintf('raster.brick==%s\n', raster.brick))
+cat(sprintf('raster.rotate==%s\n', raster.rotate))
+cat(sprintf('raster.plot==%s\n', raster.plot))
+cat(sprintf('in.fp==%s\n', in.fp))
+cat(sprintf('in.datavar.name==%s\n', in.datavar.name))
+cat(sprintf('in.datavar.type==%s\n', in.datavar.type))
+cat(sprintf('in.datavar.band==%s\n', in.datavar.band))
+cat(sprintf('out.fp==%s\n', out.fp))
+cat(sprintf('out.datavar.name==%s\n', out.datavar.name))
+cat(sprintf('out.datavar.na==%s\n', out.datavar.na))
+cat(sprintf('out.datavar.unit==%s\n', out.datavar.unit))
+cat(sprintf('out.datavar.longname==%s\n', out.datavar.longname))
+cat(sprintf('out.datavar.coord.x.name==%s\n', out.datavar.coord.x.name))
+cat(sprintf('out.datavar.coord.y.name==%s\n', out.datavar.coord.y.name))
+cat(sprintf('out.datavar.coord.z.name==%s\n', out.datavar.coord.z.name))
+cat(sprintf('out.datavar.coord.z.unit==%s\n', out.datavar.coord.z.unit))
+cat(sprintf('pdf.fp==%s\n', pdf.fp))
+cat(sprintf('pdf.height==%s\n', pdf.height))
+cat(sprintf('pdf.width==%s\n', pdf.width))
+#   end debugging
+
+raster.brick <- as.logical(raster.brick)
+raster.rotate <- as.logical(raster.rotate)
+raster.plot <- as.logical(raster.plot)
+out.datavar.na <- as.numeric(out.datavar.na)
+pdf.height <- as.numeric(pdf.height)
+pdf.width <- as.numeric(pdf.width)
+
+# ensure global scope, so raster.aqmeii can be accessed inside raster.plot conditional? nope, fail
+# raster.crop <- TRUE # kludge
+
+# functions-----------------------------------------------------------
+
+# fails to fix the conditional scoping problem below
+## plot <- function(to.plot, sigdigs, raster.to.plot) { # plot(to.plot, raster.to.plot)
+##   if (!to.plot) return
+##   # 5.1 Fix layer names (which are plot labels in `layerplot`)
+
+##   # TODO: since using rasterVis::levelplot, get/set names with `layerNames`
+##   raster.to.plot.df <- as.data.frame(raster.to.plot)
+##   names(raster.to.plot.df) <- formatC(
+##     as.numeric(sub(x=names(raster.to.plot.df), pattern="X", replacement="", fixed=TRUE)),
+##     format="g", digits=sigdigs)
+##   # start debugging
+##   cat(sprintf('layers.n==%s\n', layers.n))
+##   ## raster.to.plot
+##   ## raster.to.plot.df
+##   ## cat('head and tail of raster.to.plot.df:\n')
+##   ## head(raster.to.plot.df)
+##   ## tail(raster.to.plot.df)
+##   cat('head and tail of names(raster.to.plot.df):\n')
+##   head(names(raster.to.plot.df))
+##   tail(names(raster.to.plot.df))
+##   #  end debugging
+
+##   # 5.2 Create world map (which rasterVis will autocrop)
+
+##   cat('Loading map packages ...\n') # debugging
+##   library(maps)
+##   library(maptools) # needed for map2SpatialLines
+##   data(wrld_simpl)      # from maptools
+
+##   # TODO: save/restore after first creation!
+##   global.map <- map('world', plot=FALSE)
+##   global.map.shp <- map2SpatialLines(global.map, proj4string=global.proj)
+
+##   # 5.3 Plot and display.
+
+##   cat('Loading plot packages ...\n') # debugging
+##   library(rasterVis)
+
+##   cat(sprintf(
+##     '%s: plotting %s may take awhile! (TODO: add progress control)\n',
+##     this.fn, pdf.fp))
+##   pdf(file=pdf.fp, width=pdf.width, height=pdf.height)
+##   rasterVis::levelplot(raster.to.plot,
+##     margin=FALSE,
+##     names.attr=names(raster.to.plot.df),
+##     layout=c(1,layers.n), # all layers in one "atmospheric" column
+##   ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
+##   # teardown
+##   dev.off()
+##   # system(sprintf('xpdf %s &', pdf.fp)) # do in driver
+##   # the plot crops the map!
+## } # end function plot(to.plot, raster.to.plot)
+
+# code----------------------------------------------------------------
+
+# 1. Load the input global datavar
+
+if (raster.brick) { # multiple vertical layers
+  raster.global <- brick(in.fp, varname=in.datavar.name, crs=global.proj)
+  layers.n <- nbands(raster.global)
+} else {
+  raster.global <- raster(in.fp, varname=in.datavar.name, crs=global.proj)
+  layers.n <- 1
+}
+# raster.global # debugging
+# e.g.
+# class       : RasterBrick 
+# dimensions  : 96, 144, 13824, 56  (nrow, ncol, ncell, nlayers)
+# resolution  : 2.5, 1.894737  (x, y)
+# extent      : -1.25, 358.75, -90.94737, 90.94737  (xmin, xmax, ymin, ymax)
+# coord. ref. : +proj=longlat +datum=WGS84 
+# data source : ~/code/regridding/MOZART_global_to_AQMEII-NA/2008N2O_restart.nc 
+# names       : X1.86787999700755, X2.35259043984115, X2.94832093641162, X3.67650110274553, X4.56168595701456, X5.63180120661855, X6.91832136362791, X8.4563922137022, X10.2849211543798, X12.4601498246193, X15.0502501055598, X18.1243494153023, X21.761005744338, X26.0491091758013, X31.0889091342688, ... 
+# hybrid_sigma_pressure: 1.86787999700755, 992.500010610456 (min, max)
+# varname     : N2O 
+
+# 2. Rotate the input global datavar if requested:
+#    zero-based longitudes cause problems for wrld_simpl (and probably other data)
+if (raster.rotate) {
+  raster.global <- rotate(raster.global,
+    overwrite=TRUE) # else levelplot does one layer per page?
+#   raster.global # debugging
+# e.g.
+# class       : RasterBrick 
+# dimensions  : 96, 144, 13824, 56  (nrow, ncol, ncell, nlayers)
+# resolution  : 2.5, 1.894737  (x, y)
+# extent      : -181.25, 178.75, -90.94737, 90.94737  (xmin, xmax, ymin, ymax)
+# coord. ref. : +proj=longlat +datum=WGS84 
+# data source : in memory
+# ...
+}
+
+# 3. Crop datavar to our region.
+# snap='out' -> output includes all cells intersected by the extent (whether or not the extent covers the cell center)
+
+# ensure global scope, so raster.aqmeii can be accessed inside raster.plot conditional? nope, fail
+# if (raster.crop) {
+  raster.aqmeii <- crop(x=raster.global, y=aqmeii.extent, snap='out')
+# ensure global scope, so raster.aqmeii can be accessed inside raster.plot conditional? nope, fail
+# raster.aqmeii <<- crop(x=raster.global, y=aqmeii.extent, snap='out') 
+
+# raster.aqmeii # debugging
+# e.g.
+# class       : RasterBrick 
+# dimensions  : 21, 30, 630, 56  (nrow, ncol, ncell, nlayers)
+# resolution  : 2.5, 1.894737  (x, y)
+# extent      : -131.25, -56.25, 20.84211, 60.63158  (xmin, xmax, ymin, ymax)
+# coord. ref. : +proj=longlat +datum=WGS84 
+# data source : in memory
+# ...
+# } # if (raster.crop)
+
+# 4. Write cropped datavar to output netCDF.
+# save to file, overwriting pre-existing
+cat(sprintf('Writing cropped raster to path=%s ...\n', out.fp)) # debugging
+if (raster.brick) {
+  # ... we have layers, i.e., z
+  writeRaster(
+    x=raster.aqmeii, filename=out.fp, varname=out.datavar.name,
+    format="CDF", NAflag=out.datavar.na, overwrite=TRUE,
+    # datavar type will truncate unless set?
+    datatype=in.datavar.type,
+    # netCDF-specific arguments
+    varunit=out.datavar.unit, longname=out.datavar.longname,
+    xname=out.datavar.coord.x.name, yname=out.datavar.coord.y.name,
+    zname=out.datavar.coord.z.name, zunit=out.datavar.coord.z.unit
+  )
+} else { # it's a RasterLayer
+  writeRaster(
+    x=raster.aqmeii, filename=out.fp, varname=out.datavar.name,
+    format="CDF", NAflag=out.datavar.na, overwrite=TRUE,
+    # datavar type will truncate unless set?
+    datatype=in.datavar.type,
+    # netCDF-specific arguments
+    varunit=out.datavar.unit, longname=out.datavar.longname,
+    xname=out.datavar.coord.x.name, yname=out.datavar.coord.y.name
+  )
+}
+
+# start debugging
+system(sprintf('ls -alh %s', out.fp))
+netCDF.stats.to.stdout(netcdf.fp=in.fp, var.name=in.datavar.name)
+netCDF.stats.to.stdout(netcdf.fp=out.fp, var.name=out.datavar.name)
+system(sprintf('ncdump -h %s', out.fp))
+# e.g.
+# For /home/rtd/code/regridding/MOZART_global_to_AQMEII-NA/2008N2O_restart_region_N2O.nc var=N2O
+#         cells=35280
+#         obs=35280
+#         min=0.0151
+#         q1=0.242
+#         med=0.485
+#         mean=0.367
+#         q3=0.486
+#         max=0.492
+
+# netcdf \2008N2O_restart_region_N2O {
+# dimensions:
+#         lon = 30 ;
+#         lat = 21 ;
+#         lev = UNLIMITED ; // (56 currently)
+# variables:
+#         double lon(lon) ;
+#                 lon:units = "degrees_east" ;
+#                 lon:long_name = "lon" ;
+#         double lat(lat) ;
+#                 lat:units = "degrees_north" ;
+#                 lat:long_name = "lat" ;
+#         int lev(lev) ;
+#                 lev:units = "hybrid_sigma_pressure" ;
+#                 lev:long_name = "lev" ;
+#         double N2O(lev, lat, lon) ;
+#                 N2O:units = "ppmV" ;
+#                 N2O:_FillValue = -1. ;
+#                 N2O:missing_value = -1. ;
+#                 N2O:long_name = "N2O" ;
+#                 N2O:projection = "+proj=longlat +datum=WGS84" ;
+#                 N2O:projection_format = "PROJ.4" ;
+#   end debugging
+
+# 5. Plot data if requested
+
+# 5.0 Massive annoyance: raster.aqmeii is out of scope inside this conditional!
+# Despite numerous attempted kludges, I cannot make plotting conditional!
+# if (raster.plot) {
+
+# 5.1 Fix layer names (which are plot labels in `layerplot`)
+
+if (raster.brick) {
+  # TODO: since using rasterVis::levelplot, get/set names with `layerNames`
+  aqmeii.df <- as.data.frame(raster.aqmeii)
+  names(aqmeii.df) <- formatC(
+    as.numeric(sub(x=names(aqmeii.df), pattern="X", replacement="", fixed=TRUE)),
+    format="g", digits=sigdigs)
+  # start debugging
+  cat(sprintf('layers.n==%s\n', layers.n))
+  ## raster.aqmeii
+  ## aqmeii.df
+  ## cat('head and tail of aqmeii.df:\n')
+  ## head(aqmeii.df)
+  ## tail(aqmeii.df)
+  cat('head and tail of names(aqmeii.df):\n')
+  head(names(aqmeii.df))
+  tail(names(aqmeii.df))
+  #  end debugging
+} # end if (raster.brick)
+
+# 5.2 Create world map (which rasterVis will autocrop)
+
+  cat('Loading map packages ...\n') # debugging
+  library(maps)
+  library(maptools) # needed for map2SpatialLines
+  data(wrld_simpl)      # from maptools
+
+  # TODO: save/restore after first creation!
+  global.map <- map('world', plot=FALSE)
+  global.map.shp <- map2SpatialLines(global.map, proj4string=global.proj)
+
+# 5.3 Plot and display.
+
+  cat('Loading plot packages ...\n') # debugging
+  library(rasterVis)
+
+  cat(sprintf(
+    '%s: plotting %s may take awhile! (TODO: add progress control)\n',
+    this.fn, pdf.fp))
+#  if (raster.plot) {
+  pdf(file=pdf.fp, width=pdf.width, height=pdf.height)
+  if (raster.brick) {
+    rasterVis::levelplot(raster.aqmeii,
+      margin=FALSE,
+      names.attr=names(aqmeii.df),
+      layout=c(1,layers.n), # all layers in one "atmospheric" column
+    ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
+  } else {
+    rasterVis::levelplot(raster.aqmeii,
+      margin=FALSE
+    ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
+    # the plot crops the map!
+  } # end if (raster.brick)
+  # teardown
+  dev.off()
+
+# } # end if (raster.plot)