Commits

Tom Roche committed c2236d4

refactored, no new function:

* single driver=GEIA_global_to_AQMEII_netCDF.sh consolidating previous *.sh
* *.r s/./_/g
* all parameters passed by environment variables

tested on terrae

  • Participants
  • Parent commits 6ecdb8a

Comments (0)

Files changed (7)

File GEIA.to.netCDF.r

-# R code to write GEIA emissions from the distributed text format
-# (space separated values, 10-line header, no comments)
-# to netCDF format. Unfortunately all hardcoded.
-
-# libraries also read separately below, for input and output processing
-library(maps)   # on tlrPanP5 as well as clusters
-library(ncdf4)  # on clusters only
-library(fields) # on tlrPanP5 as well as clusters
-
-# constants-----------------------------------------------------------
-
-this.fn <- 'GEIA.to.netCDF.r'      # TODO: get like $0
-
-test.dir <- Sys.getenv('TEST_DIR') # set in driver script; MUST exist
-pdf.dir <- Sys.getenv('PDF_DIR')   # set in driver script; MUST exist
-pdf.fp <- Sys.getenv('PDF_FP')     # set in driver script
-
-# input metadata
-# for details, see http://geiacenter.org/presentData/geiadfrm.html
-GEIA.emis.txt.dir <- test.dir      # folder containing ...
-GEIA.emis.txt.fn <- 'N2OOC90Y.1A'  # ... text file from GEIA ...
-# GEIA.emis.txt.fn <- 'N2OOC90Y.1A.short'  # ... text file from GEIA ...
-# GEIA.emis.txt.fn <- 'N2OOC90Y.1A.sorted'  # ... text file from GEIA ...
-GEIA.emis.txt.n.header <- 10       # ... with n.header lines to skip @ top
-
-# GEIA data is 1° x 1° so
-grid.lat.degree.start <- -90.0     # 90 S, scalar
-grid.lat.degree.per.cell <- 1.0
-grid.lon.degree.start <- -180.0    # 180 W, scalar
-grid.lon.degree.per.cell <- 1.0
-GEIA.emis.lat.dim <- 180.0 / grid.lat.degree.per.cell # total lats, scalar
-GEIA.emis.lon.dim <- 360.0 / grid.lon.degree.per.cell # total lons, scalar
-GEIA.emis.grids.dim <-
-  GEIA.emis.lat.dim * GEIA.emis.lon.dim # total gridcells, scalar
-
-# output metadata
-netcdf.dir <- test.dir             # folder containing ...
-netcdf.fn <- 'GEIA_N2O_oceanic.nc' # ... netCDF emissions
-netcdf.title <- 'GEIA annual oceanic N2O emissions'
-netcdf.source_file <-
-  "http://geiacenter.org/data/n2ooc90y.1a.zip, 'GEIA Inventory n2ocg90yr1.1a  18 Dec 95'"
-netcdf.Conventions <-
-  "Contact: A.F. Bouwman, RIVM; Box 1, 3720 BA Bilthoven, Netherlands; tel.31-30-2743635; fax.31-30-2744417; email lex.bouwman@rivm.nl"
-time.var.name <- "time"                # like CLM-CN, EDGAR, GFED
-time.var.long_name <- "time"           # like CLM-CN, EDGAR, GFED
-time.var.units <- "year"               # for now
-time.var.calendar <- "none: emissions attributed to no particular year"
-lat.var.name <- "lat"                  # like CLM-CN, EDGAR, GFED
-lat.var.long_name <- "latitude"        # like CLM-CN, EDGAR, GFED
-lat.var.units <- "degrees_north"       # like CLM-CN, EDGAR, GFED
-lat.var.comment <- "center_of_cell"    # like EDGAR
-lon.var.name <- "lon"                  # like CLM-CN, EDGAR, GFED
-lon.var.long_name <- "longitude"       # like CLM-CN, EDGAR, GFED
-lon.var.units <- "degrees_east"        # like CLM-CN, EDGAR, GFED
-lon.var.comment <- "center_of_cell"    # like EDGAR
-emis.var.name <- "emi_n2o"             # like EDGAR
-emis.var.long_name <- "N2O emissions"  # like CLM-CN
-emis.var.total_emi_n2o <-              # EDGAR-style,
-  3.5959E+06                           # value from header(GEIA.emis.txt.fn)
-emis.var.units <- "ton N2O-N/yr"       # value from header(GEIA.emis.txt.fn)
-emis.var._FillValue <- -999.0          # like GFED
-
-# utilities
-netCDF.stats.fp <- sprintf('%s/netCDF.stats.to.stdout.r', test.dir)
-
-# functions-----------------------------------------------------------
-
-# for metadata, see http://geiacenter.org/presentData/geiadfrm.html
-# all indices are 1-based
-col.row.to.GEIA.grid.number <- function(
-  col.index,
-  row.index
-) {
-  (row.index * 1000) + col.index # return
-}
-
-# for metadata, see http://geiacenter.org/presentData/geiadfrm.html
-# all indices are 1-based
-# ASSERT: GEIA.grid.number always> 0
-GEIA.grid.number.to.lon.lat.vec <- function(
-  grid.n
-) {
-  # TODO: throw on grid.n <= 0
-  lon.lat.vec <- double(2) # [lon, lat]
-  lon.lat.vec[1] <- ((grid.n %% 1000) - 181) + 0.5 # longitude of cell center
-  lon.lat.vec[2] <- ((grid.n %/% 1000)- 91) + 0.5  # latitude of cell center
-  lon.lat.vec # return
-}
-
-# Calculate 1-d vector index from 2-d grid indices.
-# Note "grid index" != GEIA grid number !!!
-
-# GEIA numbers with latitude major:
-# > Grid Number = (j*1000) + i
-# > j   row number starting at 1 for 90S to 89S latitude,
-#                                    -90    -89
-# ...
-# > i   column number starting at 1 for 180W to 179W longitude,
-#                                       -180    -179
-
-# We can't do that, since, IN THIS CASE,
-# |latitudes| < |longitudes|: hashes will collide.
-# (If one has more lats than lons, do oppositely!)
-# Instead, smallest grid index == smallest lat index == 1
-
-# ASSERT: grid.index always> 0
-lon.lat.vec.to.grid.index <- function(
-  lon.lat.vec, # [lon, lat]
-  n.lon,       # number of longitudes
-  n.lat        # number of latitudes
-) {
-  lon <- lon.lat.vec[1] - 0.5 # lon NOT of cell center
-  lat <- lon.lat.vec[2] - 0.5 # lat NOT of cell center
-#  ((lat - 1) * n.lat) + lon # gotta deal with lat,lon < 0
-# lat=-90 -> row=1 , lon=-180 -> col=1
-#  ((lat + 90) * n.lat) + (lon + 181) # hashes collide
-  grid.index <- ((lon + 180) * n.lon) + (lat + 91)
-  if (grid.index <= 0) {
-    cat(sprintf(
-      'ERROR: %s: (lon=%.1f,lat=%.1f) -> global grid index=%.0f <= 0\n',
-      this.fn, lon, lat, global.emis.vec.index))
-    return # TODO: end execution of caller
-  }
-  grid.index
-}
-
-# Since we have way fewer records than gridcells (36143 << 64800),
-# we need to put NAs in the missing gridcells.
-# So
-# * calculate each grid index "in order"
-# * lookup value from GEIA.emis.mx
-# * write that value if found, or NA otherwise
-# TODO: make this code {more Rish, less procedural}
-create.global.emissions.vector <- function(
-  GEIA.emis.grids.dim, # total number of all gridcells that could contain values
-  GEIA.emis.mx.rows,   # total number of rows providing values
-  GEIA.emis.mx         # data structure: grid# -> value
-) {
-  global.emis.vec <-             # retval
-#    integer(GEIA.emis.grids.dim) # int vector for all gridcells
-    numeric(GEIA.emis.grids.dim)  # vector for all gridcells
-  global.emis.vec <-              # ... with all values=NA by default
-    # TODO: make size=(max(lat)*1000) + max(lon)
-    rep(NA, GEIA.emis.grids.dim)  # (actually need more space than that!)
-  for (i.row in 1:GEIA.emis.mx.rows) {
-    # get the GEIA grid# ...
-    GEIA.emis.grid.n <- GEIA.emis.mx[i.row,1]
-    # ... and the data for that gridcell ...
-    GEIA.emis.grid.val <- GEIA.emis.mx[i.row,2]
-    # ... and the corresponding index into OUR "grid vector"
-    lon.lat.vec <- # only for debugging
-      GEIA.grid.number.to.lon.lat.vec(GEIA.emis.grid.n)
-    global.emis.vec.index <-
-      lon.lat.vec.to.grid.index(
-        lon.lat.vec, n.lon=GEIA.emis.lon.dim, n.lat=GEIA.emis.lat.dim)
-    # ASSERT: never overwrite anything but NA!
-    if (is.na(global.emis.vec[global.emis.vec.index])) {
-# start debug
-#      cat(sprintf('debug: %s: writing value=%f from GEIA grid#=%.0f to global grid index=%.0f (lon=%.1f,lat=%.1f)\n',
-#        this.fn, GEIA.emis.grid.val, GEIA.emis.grid.n, 
-#        global.emis.vec.index, lon.lat.vec[1], lon.lat.vec[2]))
-#   end debug
-      global.emis.vec[global.emis.vec.index] <- GEIA.emis.grid.val
-    } else {
-# '%i' draws error? and '%d'??
-#      cat(sprintf('ERROR: %s: attempting to write value from GEIA grid#=%i to grid index=%i (lon=%d,lat=%d) -> %f',
-      cat(sprintf('ERROR: %s: attempting to write value=%f from GEIA grid#=%.0f to global grid index=%.0f (lon=%.1f,lat=%.1f)\n',
-        this.fn, GEIA.emis.grid.val, GEIA.emis.grid.n, 
-        global.emis.vec.index, lon.lat.vec[1], lon.lat.vec[2]))
-      return # TODO: end execution of caller
-    }
-  } # end for loop over GEIA data values
-  global.emis.vec # return
-} # end function create.global.emissions.vector
-
-# TODO: refactor! this plot.layer is copy/mod from ioapi::plotLayersForTimestep.r::plot.layer
-# changes to arguments:
-# * attrs.list -> x.centers, y.centers (since these are in degrees, not km)
-# * map not passed: I'm using global map via `map(add=TRUE)` TODO: fixme
-plot.layer <- function(
-  data,             # data to plot (required)
-  title,            # string for plot title (required?)
-                    # TODO: handle when null!
-  subtitle=NULL,    # string for plot subtitle
-  x.centers,        # centers of abscissa points
-  y.centers,        # centers of ordinate points
-  q.vec=NULL,       # for quantiles
-  colors
-) {
-  if (sum(!is.na(data)) && (!is.null(q.vec))) {
-    plot.list <- list(x=x.centers, y=y.centers, z=data)
-    quantiles <- quantile(c(data), q.vec, na.rm=TRUE)
-    quantiles.formatted <- format(as.numeric(quantiles), digits=3)
-# start debugging
-#      print(paste('Non-null image.plot for source layer==', i.layer, ', quantile bounds=='))
-#      print(quantiles)
-#   end debugging
-    if (is.null(subtitle)) {
-      image.plot(plot.list, xlab="", ylab="", axes=F, col=colors(100),
-        axis.args=list(at=quantiles, labels=quantiles.formatted),
-        main=title)
-    } else {
-      image.plot(plot.list, xlab="", ylab="", axes=F, col=colors(100),
-        axis.args=list(at=quantiles, labels=quantiles.formatted),
-        main=title, sub=subtitle)
-    }
-    map(add=TRUE)
-  } else {
-# debugging
-#      print(paste('Null image.plot for source layer=', i.layer))
-    if (is.null(subtitle)) {
-      plot(0, type="n", axes=F, xlab="", ylab="",
-        xlim=range(x.centers), ylim=range(y.centers),
-        main=title)
-    } else {
-      plot(0, type="n", axes=F, xlab="", ylab="",
-        xlim=range(x.centers), ylim=range(y.centers),
-        main=title, sub=subtitle)
-    }
-    map(add=TRUE)
-  } # end testing data
-} # end function plot.layer
-
-# TODO: refactor! subtitle.stats is copied from ioapi::plotLayersForTimestep.r
-subtitle.stats <- function(vec) {
-  return.str <- ""
-  # is it numeric, and not empty?
-  if (is.numeric(vec) && sum(!is.na(vec))) {
-#    unsparse.vec <- subset(vec, !is.na(vec)) # fail: intended for interactive use
-#    unsparse.vec <- na.omit(vec) # fail: omits all *rows* containing an NA!
-    grids <- length(vec)
-    grids.str <- sprintf('(of cells=%i)', grids)
-    unsparse.vec <- vec[!is.na(vec)]
-    obs <- length(unsparse.vec)
-    obs.str <- sprintf('obs=%i', obs)
-    # use constants defined above. TODO: compute these once!
-    max.str <- sprintf(max.str, max(unsparse.vec))
-    med.str <- sprintf(med.str, median(unsparse.vec))
-    min.str <- sprintf(min.str, min(unsparse.vec))
-    return.str <-
-      sprintf('%s %s: %s, %s, %s',
-              obs.str, grids.str, min.str, med.str, max.str)
-  } else {
-    return.str <-"no data"
-  }
-  return.str
-} # end function subtitle.stats
-
-# code----------------------------------------------------------------
-
-# process input
-library(maps)  # on tlrPanP5 as well as clusters
-
-# input file path
-GEIA.emis.txt.fp <- sprintf('%s/%s', GEIA.emis.txt.dir, GEIA.emis.txt.fn)
-# columns are grid#, mass
-GEIA.emis.mx  <-
-  as.matrix(read.table(GEIA.emis.txt.fp, skip=GEIA.emis.txt.n.header))
-# mask zeros? no, use NA for non-ocean areas
-# GEIA.emis.mx[GEIA.emis.mx == 0] <- NA
-# <simple input check>
-dim(GEIA.emis.mx) ## [1] 36143     2
-# start debug
-GEIA.emis.mx.rows <- dim(GEIA.emis.mx)[1]
-if (GEIA.emis.mx.rows > GEIA.emis.grids.dim) {
-  cat(sprintf('ERROR: %s: GEIA.emis.mx.rows=%.0d > GEIA.emis.grids.dim=%.0d\n',
-    this.fn, GEIA.emis.mx.rows, GEIA.emis.grids.dim))
-} else {
-  cat(sprintf('debug: %s: GEIA.emis.mx.rows < GEIA.emis.grids.dim\n',
-    this.fn))
-}
-#   end debug
-# </simple input check>
-
-global.emis.vec <-
-  create.global.emissions.vector(
-    GEIA.emis.grids.dim, GEIA.emis.mx.rows, GEIA.emis.mx)
-
-# <visual input check>
-# Need sorted lat and lon vectors: we know what those are a priori
-# Add 0.5 since grid centers
-lon.vec <- 0.5 +
-  seq(from=grid.lon.degree.start, by=grid.lon.degree.per.cell, length.out=GEIA.emis.lon.dim)
-lat.vec <- 0.5 +
-  seq(from=grid.lat.degree.start, by=grid.lat.degree.per.cell, length.out=GEIA.emis.lat.dim)
-
-# Create emissions matrix corresponding to those dimensional vectors
-# (i.e., global.emis.mx is the "projection" of global.emis.vec)
-# First, create empty global.emis.mx? No, fill from global.emis.vec.
-# NO: I cannot just fill global.emis.mx from global.emis.vec:
-# latter's/GEIA's grid numbering system ensures 1000 lons per lat!
-# Which overflows the "real-spatial" global.emis.mx :-(
-# So I need to fill global.emis.mx using a for loop to decode the grid indices :-(
-# (but at least I can fill in whatever direction I want :-)
-global.emis.mx <- matrix(
-  rep(NA, GEIA.emis.grids.dim), nrow=GEIA.emis.lat.dim, ncol=GEIA.emis.lon.dim)
-
-# 1: works if subsequently transposed: TODO: FIXME
-for (i.lon in 1:GEIA.emis.lon.dim) {
-  for (i.lat in 1:GEIA.emis.lat.dim) {
-# 2: fails with 'dimensions of z are not length(x)(-1) times length(y)(-1)'
-# for (i.lat in 1:GEIA.emis.lat.dim) {
-#   for (i.lon in 1:GEIA.emis.lon.dim) {
-# 3: fails with 'dimensions of z are not length(x)(-1) times length(y)(-1)'
-# for (i.lon in GEIA.emis.lon.dim:1) {
-#   for (i.lat in GEIA.emis.lat.dim:1) {
-# 4: fails with 'dimensions of z are not length(x)(-1) times length(y)(-1)'
-# for (i.lat in GEIA.emis.lat.dim:1) {
-#   for (i.lon in GEIA.emis.lon.dim:1) {
-    lon <- lon.vec[i.lon]
-    lat <- lat.vec[i.lat]
-    GEIA.emis.grid.val <-
-      global.emis.vec[
-        lon.lat.vec.to.grid.index(c(lon, lat),
-          n.lon=GEIA.emis.lon.dim, n.lat=GEIA.emis.lat.dim)]
-    if (!is.na(GEIA.emis.grid.val)) {
-      if (is.na(global.emis.mx[i.lat, i.lon])) {
-        global.emis.mx[i.lat, i.lon] <- GEIA.emis.grid.val
-# start debug
-#        cat(sprintf(
-#          'debug: %s: writing val=%f to global.emis.mx[%.0f,%.0f] for grid center=[%f,%f]\n',
-#          this.fn, GEIA.emis.grid.val, i.lon, i.lat, lon, lat))
-#   end debug
-      } else {
-        # error if overwriting presumably-previously-written non-NA!
-        cat(sprintf(
-          'ERROR: %s: overwriting val=%f with val=%f at global.emis.mx[%.0f,%.0f] for grid center=[%f,%f]\n',
-          this.fn, global.emis.mx[i.lat, i.lon], GEIA.emis.grid.val,
-          i.lon, i.lat, lon, lat))
-        return # TODO: abend
-      } # end testing target != NA (thus not previously written)
-    } # end testing source != NA (don't write if is.na(lookup)
-  } # end for loop over lats
-} # end for loop over lons
-
-
-# start debugging
-# pdf(file=pdf.fp, width=5.5, height=4.25)
-# # 1: TODO: FIXME: why do I need to transpose global.emis.mx?
-# image(lon.vec, lat.vec, t(global.emis.mx))
-# # 2,3,4: how it should work ?!?
-# # image(lon.vec, lat.vec, global.emis.mx)
-# map(add=TRUE)
-# # dev.off() # leave on for subsequent plot; this will be p1/2
-#   end debugging
-# </visual input check>
-
-# write output to netCDF
-library(ncdf4)
-
-# output file path
-netcdf.fp <- sprintf('%s/%s', netcdf.dir, netcdf.fn)
-
-# create dimensions and dimensional variables
-time.vec <- c(0) # annual value, corresponding to no specific year
-time.dim <- ncdim_def(
-  name=time.var.name,
-  units=time.var.units,
-  vals=time.vec,
-  unlim=TRUE,
-  create_dimvar=TRUE,
-  calendar=time.var.calendar,
-  longname=time.var.long_name)
-
-lon.dim <- ncdim_def(
-  name=lon.var.name,
-  units=lon.var.units,
-  vals=lon.vec,
-  unlim=FALSE,
-  create_dimvar=TRUE,
-  longname=lon.var.long_name)
-
-lat.dim <- ncdim_def(
-  name=lat.var.name,
-  units=lat.var.units,
-  vals=lat.vec,
-  unlim=FALSE,
-  create_dimvar=TRUE,
-  longname=lat.var.long_name)
-
-# create data variable (as container--can't put data until we have a file)
-emis.var <- ncvar_def(
-  name=emis.var.name,
-  units=emis.var.units,
-#  dim=c(time.dim, lat.dim, lon.dim),
-#  dim=list(time.dim, lat.dim, lon.dim),
-# note dim order desired for result=var(time, lat, lon)
-  dim=list(lon.dim, lat.dim, time.dim),
-  missval=as.double(emis.var._FillValue),
-  longname=emis.var.long_name,
-  prec="double")
-
-# get current time for creation_date
-# system(intern=TRUE) -> return char vector, one member per output line)
-netcdf.timestamp <- system('date', intern=TRUE)
-
-# create netCDF file
-netcdf.file <- nc_create(
-  filename=netcdf.fn,
-#  vars=list(emis.var),
-#  verbose=TRUE)
-  vars=list(emis.var))
-
-# Write data to data variable: gotta have file first.
-# Gotta convert 2d global.emis.mx[lat,lon] to 3d global.emis.arr[time,lat,lon]
-# Do this before adding _FillValue to prevent:
-# > Error in R_nc4_put_vara_double: NetCDF: Not a valid data type or _FillValue type mismatch
-## global.emis.arr <- global.emis.mx
-## dim(global.emis.arr) <- c(1, dim(global.emis.mx))
-## global.emis.arr[1,,] <- global.emis.mx
-
-# Note
-# * global.emis.mx[lat,lon]
-# * datavar needs [lon, lat, time] (with time *last*)
-
-ncvar_put(
-  nc=netcdf.file,
-  varid=emis.var,
-#  vals=global.emis.arr,
-  vals=t(global.emis.mx),
-#  start=rep.int(1, length(dim(global.emis.arr))),
-  start=c(1, 1, 1),
-#  count=dim(global.emis.arr))
-  count=c(-1,-1, 1)) # -1 -> all data
-
-# Write netCDF attributes
-# Note: can't pass *.dim as varid, even though these are coordinate vars :-(
-
-# add datavar attributes
-ncatt_put(
-  nc=netcdf.file,
-#  varid=lon.var,
-  varid=lon.var.name,
-  attname="comment",
-  attval=lon.var.comment,
-  prec="text")
-
-ncatt_put(
-  nc=netcdf.file,
-#  varid=lat.var,
-  varid=lat.var.name,
-  attname="comment",
-  attval=lat.var.comment,
-  prec="text")
-
-# put _FillValue after putting data!
-ncatt_put(
-  nc=netcdf.file,
-  varid=emis.var,
-  attname="_FillValue",
-  attval=emis.var._FillValue,
-  prec="float") # why is "emi_n2o:missing_value = -999."?
-
-# add global attributes (varid=0)
-ncatt_put(
-  nc=netcdf.file,
-  varid=0,
-  attname="creation_date",
-  attval=netcdf.timestamp,
-  prec="text")
-
-ncatt_put(
-  nc=netcdf.file,
-  varid=0,
-  attname="source_file",
-  attval=netcdf.source_file,
-  prec="text")
-
-ncatt_put(
-  nc=netcdf.file,
-  varid=0,
-  attname="Conventions",
-  attval=netcdf.Conventions,
-  prec="text")
-
-# flush to file (there may not be data on disk before this point)
-# nc_sync(netcdf.file) # so we don't hafta reopen the file, below
-# Nope: per David W. Pierce Mon, 27 Aug 2012 21:35:35 -0700, ncsync is not enough
-nc_close(netcdf.file)
-
-# NOTE: you must assign when you nc_open! (though not when you nc_close)
-# If you do
-## nc_open(netcdf.fp,
-##         write=FALSE,    # will only read below
-##         readunlim=TRUE) # it's a small file
-# you will get this error on your first attempt to ncvar_get
-# (though no error on ncvar_put, annoyingly enough):
-# > Error in if (nc$var[[li]]$hasAddOffset) addOffset = nc$var[[li]]$addOffset else addOffset = 0 :
-# > argument is of length zero
-# MERCI BEAUCOUP, Pascal Oettli!
-netcdf.file <- nc_open(
-  filename=netcdf.fp,
-  write=FALSE,    # will only read below
-  readunlim=TRUE) # it's a small file
-
-# <simple output check>
-system(sprintf('ls -alth %s', netcdf.fp))
-system(sprintf('ncdump -h %s', netcdf.fp))
-# TODO: call sample stats from netCDF.stats.to.stdout.r via function, not like this!
-# system('Rscript ./netCDF.stats.to.stdout.r netcdf.fp=./GEIA_N2O_oceanic.nc var.name=emi_n2o')
-source(netCDF.stats.fp)
-netCDF.stats.to.stdout(
-  netcdf.fp='./GEIA_N2O_oceanic.nc',
-  var.name='emi_n2o')
-# </simple output check>
-
-# <visual output check>
-# TODO: do plot-related refactoring! allow to work with projects={ioapi, this}
-# <copied from plotLayersForTimestep.r>
-library(fields)
-# double-sprintf-ing to set precision by constant: cool or brittle?
-stats.precision <- 3 # sigdigs to use for min, median, max of obs
-stat.str <- sprintf('%%.%ig', stats.precision)
-# use these in function=subtitle.stats as sprintf inputs
-max.str <- sprintf('max=%s', stat.str)
-med.str <- sprintf('med=%s', stat.str)
-min.str <- sprintf('min=%s', stat.str)
-# </copied from plotLayersForTimestep.r>
-
-# Get the data out of the datavar, to test reusability
-target.data <- ncvar_get(
-  nc=netcdf.file,
-  varid=emis.var,
-  # read all the data
-  start=rep(1, emis.var$ndims),
-  count=rep(-1, emis.var$ndims)) 
-
-# <simple output check/>
-dim(target.data) # n.lon, n.lat
-
-# <copy/mod from windowEmissions.r>
-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)
-probabilities.vec <- seq(0, 1, 1.0/(length(palette.vec) - 1))
-# </copy/mod from windowEmissions.r>
-
-pdf(file=pdf.fp, width=5.5, height=4.25)
-# <copy/mod from plotLayersForTimestep.r>
-plot.layer(target.data,
-  title=netcdf.title,
-  subtitle=subtitle.stats(target.data),
-  x.centers=lon.vec,
-  y.centers=lat.vec,
-  q.vec=probabilities.vec,
-  colors=colors)
-# </copy/mod from plotLayersForTimestep.r>
-map(add=TRUE)
-# </visual output check>
-
-# teardown
-dev.off()
-nc_close(netcdf.file)
-quit(save="no") # just exit

File GEIA_global_to_AQMEII_netCDF.sh

+#!/usr/bin/env bash
+
+# Requires R.
+# Driver for [GEIA_to_netCDF.r,global_to_AQMEII.r], converting
+# GEIA marine N2O emissions from global native format to netCDF over AQMEII.
+# Configure as needed for your platform.
+
+# ----------------------------------------------------------------------
+# constants with some simple manipulations
+# ----------------------------------------------------------------------
+
+THIS_FN="$(basename $0)"
+# will do the real work
+export CALL_REFORMAT_FN='GEIA_to_netCDF.r'
+export CALL_REGRID_FN='global_to_AQMEII.r'
+
+# 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_REFORMAT_FP="${WORK_DIR}/${CALL_REFORMAT_FN}"
+export CALL_REGRID_FP="${WORK_DIR}/${CALL_REGRID_FN}"
+
+# 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}"
+
+# Raw input data is GEIA marine N2O emissions in native format.
+# For details, see http://geiacenter.org/presentData/geiadfrm.html
+# Get from my repository (stable location, no unzip)
+export GEIA_RAW_SOURCE_FILE="http://geiacenter.org/data/n2ooc90y.1a.zip, 'GEIA Inventory n2ocg90yr1.1a  18 Dec 95'"
+export GEIA_RAW_CONVENTIONS='Contact: A.F. Bouwman, RIVM; Box 1, 3720 BA Bilthoven, Netherlands; tel.31-30-2743635; fax.31-30-2744417; email lex.bouwman@rivm.nl'
+GEIA_RAW_URI='https://github.com/downloads/TomRoche/GEIA_to_netCDF/N2OOC90Y.1A'
+GEIA_RAW_FN="$(basename ${GEIA_RAW_URI})"
+export GEIA_RAW_FP="${WORK_DIR}/${GEIA_RAW_FN}"
+
+export GEIA_RAW_HEADER_N='10'       # header.n lines to skip @ top
+# GEIA data is 1° x 1°, so ...
+export GEIA_RAW_LAT_DEGREE_START='-90.0'   # 90 S, scalar
+export GEIA_RAW_LAT_DEGREE_PER_CELL='1.0'
+export GEIA_RAW_LON_DEGREE_START='-180.0'  # 180 W, scalar
+export GEIA_RAW_LON_DEGREE_PER_CELL='1.0'
+
+GEIA_RAW_PDF_FN="geia_raw_$(date +${DATE_FORMAT}).pdf"
+export GEIA_RAW_PDF_FP="${WORK_DIR}/${GEIA_RAW_PDF_FN}" # path to PDF output
+
+# Intermediate product reformated to netCDF
+
+GEIA_REFORMAT_URI='https://bitbucket.org/tlroche/geia_regrid/downloads/GEIA_N2O_oceanic.nc'
+GEIA_REFORMAT_FN="$(basename ${GEIA_REFORMAT_URI})"
+export GEIA_REFORMAT_FP="${WORK_DIR}/${GEIA_REFORMAT_FN}"
+
+# TODO: automate getting this metadata
+export GEIA_REFORMAT_DATAVAR_NAME='emi_n2o'            # like EDGAR
+export GEIA_REFORMAT_DATAVAR_LONGNAME='N2O emissions'  # like CLM-CN
+export GEIA_REFORMAT_DATAVAR_UNIT='ton N2O-N/yr'       # value from header(GEIA.emis.txt.fn)
+export GEIA_REFORMAT_DATAVAR_NA='-999.0'               # like GFED
+export GEIA_REFORMAT_DATAVAR_BAND='3' # index of dim=TIME in emi_n2o(time, lat, lon)
+export GEIA_REFORMAT_DATAVAR_TOTAL='3.5959E+06'        # value from header(GEIA.emis.txt.fn)
+export GEIA_REFORMAT_TIME_VAR_NAME='time'              # like CLM-CN, EDGAR, GFED
+export GEIA_REFORMAT_TIME_VAR_LONG_NAME='time'         # like CLM-CN, EDGAR, GFED
+export GEIA_REFORMAT_TIME_VAR_UNITS='year'             # for now
+export GEIA_REFORMAT_TIME_VAR_CALENDAR='none: emissions attributed to no particular year'
+export GEIA_REFORMAT_LAT_VAR_NAME='lat'                # like CLM-CN, EDGAR, GFED
+export GEIA_REFORMAT_LAT_VAR_LONG_NAME='latitude'      # like CLM-CN, EDGAR, GFED
+export GEIA_REFORMAT_LAT_VAR_UNITS='degrees_north'     # like CLM-CN, EDGAR, GFED
+export GEIA_REFORMAT_LAT_VAR_COMMENT='center_of_cell'  # like EDGAR
+export GEIA_REFORMAT_LON_VAR_NAME='lon'                # like CLM-CN, EDGAR, GFED
+export GEIA_REFORMAT_LON_VAR_LONG_NAME='longitude'     # like CLM-CN, EDGAR, GFED
+export GEIA_REFORMAT_LON_VAR_UNITS='degrees_east'      # like CLM-CN, EDGAR, GFED
+export GEIA_REFORMAT_LON_VAR_COMMENT='center_of_cell'  # like EDGAR
+
+GEIA_REFORMAT_PDF_FN="geia_reformat_$(date +${DATE_FORMAT}).pdf"
+export GEIA_REFORMAT_PDF_FP="${WORK_DIR}/${GEIA_REFORMAT_PDF_FN}"
+export GEIA_REFORMAT_PDF_TITLE='GEIA annual oceanic N2O emissions'
+
+# 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!
+export TEMPLATE_VAR_NAME='emi_n2o'
+
+TEMPLATE_INPUT_URI='https://bitbucket.org/tlroche/geia_regrid/downloads/emis_mole_all_20080101_12US1_cmaq_cb05_soa_2008ab_08c.EXTENTS_INPUT.nc'
+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
+
+# Output is regridded netCDF.
+
+GEIA_REGRID_URI='https://bitbucket.org/tlroche/geia_regrid/downloads/GEIA_N2O_oceanic_regrid.nc'
+GEIA_REGRID_FN="$(basename ${GEIA_REGRID_URI})"
+export GEIA_REGRID_FP="${WORK_DIR}/${GEIA_REGRID_FN}"
+# TODO: automate getting this metadata
+export GEIA_REGRID_X_VAR_NAME='COL'
+export GEIA_REGRID_Y_VAR_NAME='ROW'
+export GEIA_REGRID_DATAVAR_UNIT="${GEIA_REFORMAT_DATAVAR_UNIT}"
+# export GEIA_REGRID_PDF_TITLE="GEIA annual oceanic N2O emissions regridded to AQMEII-NA\n${GEIA_REGRID_DATAVAR_UNIT}"
+export GEIA_REGRID_PDF_TITLE='GEIA annual oceanic N2O emissions regridded to AQMEII-NA'
+
+GEIA_REGRID_PDF_FN="geia_regrid_$(date +${DATE_FORMAT}).pdf"
+export GEIA_REGRID_PDF_FP="${WORK_DIR}/${GEIA_REGRID_PDF_FN}"
+
+# map*.dat is used to create a map for image.plot-ing
+# use one of following, depending on units
+
+# units=km
+# MAP_TABLE_URI='https://bitbucket.org/tlroche/geia_regrid/downloads/map.CMAQkm.world.dat'
+# units=m
+MAP_TABLE_URI='https://bitbucket.org/tlroche/geia_regrid/downloads/map.CMAQm.world.dat'
+MAP_TABLE_FN="$(basename ${MAP_TABLE_URI})"
+export MAP_TABLE_FP="${WORK_DIR}/${MAP_TABLE_FN}"
+
+# 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 "${MAP_TABLE_FP}" ]] ; then
+  echo -e "${THIS_FN}: ERROR: MAP_TABLE_FP not defined"
+  exit 1
+fi
+if [[ ! -r "${MAP_TABLE_FP}" ]] ; then
+  for CMD in \
+    "wget --no-check-certificate -c -O ${MAP_TABLE_FP} ${MAP_TABLE_URI}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+if [[ -z "${STAT_SCRIPT_FP}" ]] ; then
+  echo -e "${THIS_FN}: ERROR: STAT_SCRIPT_FP not defined"
+  exit 2
+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
+
+if [[ -z "${PLOT_SCRIPT_FP}" ]] ; then
+  echo -e "${THIS_FN}: ERROR: PLOT_SCRIPT_FP not defined"
+  exit 3
+fi
+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 [[ -z "${GEIA_RAW_FP}" ]] ; then
+  echo -e "${THIS_FN}: ERROR: GEIA_RAW_FP not defined"
+  exit 4
+fi
+if [[ ! -r "${GEIA_RAW_FP}" ]] ; then
+  for CMD in \
+    "wget --no-check-certificate -c -O ${GEIA_RAW_FP} ${GEIA_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 "${GEIA_REFORMAT_FP}" ]] ; then
+  echo -e "${THIS_FN}: ERROR: GEIA_REFORMAT_FP not defined"
+  exit 6
+fi
+# if [[ ! -r "${GEIA_REFORMAT_FP}" ]] ; then
+#   for CMD in \
+#     "wget --no-check-certificate -c -O ${GEIA_REFORMAT_FP} ${GEIA_REFORMAT_URI}" \
+#   ; do
+#     echo -e "$ ${CMD}"
+#     eval "${CMD}"
+#   done
+# fi
+
+# only for testing. TODO: flag control
+if [[ -z "${GEIA_REGRID_FP}" ]] ; then
+  echo -e "${THIS_FN}: ERROR: GEIA_REGRID_FP not defined"
+  exit 7
+fi
+# if [[ ! -r "${GEIA_REGRID_FP}" ]] ; then
+#   for CMD in \
+#     "wget --no-check-certificate -c -O ${GEIA_REGRID_FP} ${GEIA_REGRID_URI}" \
+#   ; do
+#     echo -e "$ ${CMD}"
+#     eval "${CMD}"
+#   done
+# fi
+
+# ----------------------------------------------------------------------
+# payload
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# convert from native format to netCDF
+# ----------------------------------------------------------------------
+
+# If this fails ...
+"${RSCRIPT}" "${CALL_REFORMAT_FP}"
+# ... just start R ...
+# "${R}"
+# ... and source ${CALL_REFORMAT_FP}, e.g.,
+# source('....r')
+
+# After exiting R, show cwd and display output PDF.
+for CMD in \
+  "ls -alht ${WORK_DIR}" \
+  "${PDF_VIEWER} ${GEIA_REFORMAT_PDF_FP} &" \
+; do
+  echo -e "$ ${CMD}"
+  eval "${CMD}"
+done
+
+# ----------------------------------------------------------------------
+# 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 output PDF.
+for CMD in \
+  "ls -alht ${WORK_DIR}" \
+  "${PDF_VIEWER} ${GEIA_REGRID_PDF_FP} &" \
+; do
+  echo -e "$ ${CMD}"
+  eval "${CMD}"
+done

File GEIA_to_netCDF.r

+# R code to write GEIA emissions from the distributed text format
+# (space separated values, 10-line header, no comments)
+# to netCDF format. Unfortunately all hardcoded.
+
+# ----------------------------------------------------------------------
+# constants
+# ----------------------------------------------------------------------
+
+# set in driver script; MUST exist
+this.fn <- Sys.getenv('CALL_REFORMAT_FN')
+
+work.dir <- Sys.getenv('WORK_DIR') 
+pdf.dir <- Sys.getenv('PDF_DIR')   # set in driver script; MUST exist
+pdf.fp <- Sys.getenv('GEIA_REFORMAT_PDF_FP')     # set in driver script
+
+# input metadata
+# input file path
+GEIA.emis.txt.fp <- Sys.getenv('GEIA_RAW_FP') # ... text file from GEIA ...
+# GEIA.emis.txt.fn <- 'N2OOC90Y.1A.short'     # ... text file from GEIA ...
+# GEIA.emis.txt.fn <- 'N2OOC90Y.1A.sorted'    # ... text file from GEIA ...
+GEIA.emis.txt.n.header <- # n lines to skip @ top
+  as.numeric( Sys.getenv('GEIA_RAW_HEADER_N'))
+
+grid.lat.degree.start <- as.numeric(Sys.getenv('GEIA_RAW_LAT_DEGREE_START'))
+grid.lat.degree.per.cell <- as.numeric(Sys.getenv('GEIA_RAW_LAT_DEGREE_PER_CELL'))
+grid.lon.degree.start <- as.numeric(Sys.getenv('GEIA_RAW_LON_DEGREE_START'))
+grid.lon.degree.per.cell <- as.numeric(Sys.getenv('GEIA_RAW_LON_DEGREE_PER_CELL'))
+GEIA.emis.lat.dim <- 180.0 / grid.lat.degree.per.cell # total lats, scalar
+GEIA.emis.lon.dim <- 360.0 / grid.lon.degree.per.cell # total lons, scalar
+GEIA.emis.grids.dim <-
+  GEIA.emis.lat.dim * GEIA.emis.lon.dim # total gridcells, scalar
+
+# output metadata
+netcdf.fp <- Sys.getenv('GEIA_REFORMAT_FP')
+netcdf.title <- Sys.getenv('GEIA_REFORMAT_PDF_TITLE')
+netcdf.source_file <- Sys.getenv('GEIA_RAW_SOURCE_FILE')
+netcdf.Conventions <- Sys.getenv('GEIA_RAW_CONVENTIONS')
+time.var.name <- Sys.getenv('GEIA_REFORMAT_TIME_VAR_NAME')
+time.var.long_name <- Sys.getenv('GEIA_REFORMAT_TIME_VAR_LONG_NAME')
+time.var.units <- Sys.getenv('GEIA_REFORMAT_TIME_VAR_UNITS')
+time.var.calendar <- Sys.getenv('GEIA_REFORMAT_TIME_VAR_CALENDAR')
+lat.var.name <- Sys.getenv('GEIA_REFORMAT_LAT_VAR_NAME')
+lat.var.long_name <- Sys.getenv('GEIA_REFORMAT_LAT_VAR_LONG_NAME')
+lat.var.units <- Sys.getenv('GEIA_REFORMAT_LAT_VAR_UNITS')
+lat.var.comment <- Sys.getenv('GEIA_REFORMAT_LAT_VAR_COMMENT')
+lon.var.name <- Sys.getenv('GEIA_REFORMAT_LON_VAR_NAME')
+lon.var.long_name <- Sys.getenv('GEIA_REFORMAT_LON_VAR_LONG_NAME')
+lon.var.units <- Sys.getenv('GEIA_REFORMAT_LON_VAR_UNITS')
+lon.var.comment <- Sys.getenv('GEIA_REFORMAT_LON_VAR_COMMENT')
+emis.var.name <- Sys.getenv('GEIA_REFORMAT_DATAVAR_NAME')
+emis.var.long_name <- Sys.getenv('GEIA_REFORMAT_DATAVAR_LONGNAME')
+emis.var.units <- Sys.getenv('GEIA_REFORMAT_DATAVAR_UNIT')
+emis.var._FillValue <- as.numeric(Sys.getenv('GEIA_REFORMAT_DATAVAR_NA'))
+
+# utilities
+stat.script.fp <- Sys.getenv('STAT_SCRIPT_FP')
+plot.script.fp <- Sys.getenv('PLOT_SCRIPT_FP')
+
+# ----------------------------------------------------------------------
+# functions
+# ----------------------------------------------------------------------
+
+# for metadata, see http://geiacenter.org/presentData/geiadfrm.html
+# all indices are 1-based
+col.row.to.GEIA.grid.number <- function(
+  col.index,
+  row.index
+) {
+  (row.index * 1000) + col.index # return
+}
+
+# for metadata, see http://geiacenter.org/presentData/geiadfrm.html
+# all indices are 1-based
+# ASSERT: GEIA.grid.number always> 0
+GEIA.grid.number.to.lon.lat.vec <- function(
+  grid.n
+) {
+  # TODO: throw on grid.n <= 0
+  lon.lat.vec <- double(2) # [lon, lat]
+  lon.lat.vec[1] <- ((grid.n %% 1000) - 181) + 0.5 # longitude of cell center
+  lon.lat.vec[2] <- ((grid.n %/% 1000)- 91) + 0.5  # latitude of cell center
+  lon.lat.vec # return
+}
+
+# Calculate 1-d vector index from 2-d grid indices.
+# Note "grid index" != GEIA grid number !!!
+
+# GEIA numbers with latitude major:
+# > Grid Number = (j*1000) + i
+# > j   row number starting at 1 for 90S to 89S latitude,
+#                                    -90    -89
+# ...
+# > i   column number starting at 1 for 180W to 179W longitude,
+#                                       -180    -179
+
+# We can't do that, since, IN THIS CASE,
+# |latitudes| < |longitudes|: hashes will collide.
+# (If one has more lats than lons, do oppositely!)
+# Instead, smallest grid index == smallest lat index == 1
+
+# ASSERT: grid.index always> 0
+lon.lat.vec.to.grid.index <- function(
+  lon.lat.vec, # [lon, lat]
+  n.lon,       # number of longitudes
+  n.lat        # number of latitudes
+) {
+  lon <- lon.lat.vec[1] - 0.5 # lon NOT of cell center
+  lat <- lon.lat.vec[2] - 0.5 # lat NOT of cell center
+#  ((lat - 1) * n.lat) + lon # gotta deal with lat,lon < 0
+# lat=-90 -> row=1 , lon=-180 -> col=1
+#  ((lat + 90) * n.lat) + (lon + 181) # hashes collide
+  grid.index <- ((lon + 180) * n.lon) + (lat + 91)
+  if (grid.index <= 0) {
+    cat(sprintf(
+      '%s: ERROR: (lon=%.1f,lat=%.1f) -> global grid index=%.0f <= 0\n',
+      this.fn, lon, lat, global.emis.vec.index))
+    return # TODO: end execution of caller
+  }
+  grid.index
+}
+
+# Since we have way fewer records than gridcells (36143 << 64800),
+# we need to put NAs in the missing gridcells.
+# So
+# * calculate each grid index "in order"
+# * lookup value from GEIA.emis.mx
+# * write that value if found, or NA otherwise
+# TODO: make this code {more Rish, less procedural}
+create.global.emissions.vector <- function(
+  GEIA.emis.grids.dim, # total number of all gridcells that could contain values
+  GEIA.emis.mx.rows,   # total number of rows providing values
+  GEIA.emis.mx         # data structure: grid# -> value
+) {
+  global.emis.vec <-             # retval
+#    integer(GEIA.emis.grids.dim) # int vector for all gridcells
+    numeric(GEIA.emis.grids.dim)  # vector for all gridcells
+  global.emis.vec <-              # ... with all values=NA by default
+    # TODO: make size=(max(lat)*1000) + max(lon)
+    rep(NA, GEIA.emis.grids.dim)  # (actually need more space than that!)
+  for (i.row in 1:GEIA.emis.mx.rows) {
+    # get the GEIA grid# ...
+    GEIA.emis.grid.n <- GEIA.emis.mx[i.row,1]
+    # ... and the data for that gridcell ...
+    GEIA.emis.grid.val <- GEIA.emis.mx[i.row,2]
+    # ... and the corresponding index into OUR "grid vector"
+    lon.lat.vec <- # only for debugging
+      GEIA.grid.number.to.lon.lat.vec(GEIA.emis.grid.n)
+    global.emis.vec.index <-
+      lon.lat.vec.to.grid.index(
+        lon.lat.vec, n.lon=GEIA.emis.lon.dim, n.lat=GEIA.emis.lat.dim)
+    # ASSERT: never overwrite anything but NA!
+    if (is.na(global.emis.vec[global.emis.vec.index])) {
+# start debug
+#      cat(sprintf('debug: %s: writing value=%f from GEIA grid#=%.0f to global grid index=%.0f (lon=%.1f,lat=%.1f)\n',
+#        this.fn, GEIA.emis.grid.val, GEIA.emis.grid.n, 
+#        global.emis.vec.index, lon.lat.vec[1], lon.lat.vec[2]))
+#   end debug
+      global.emis.vec[global.emis.vec.index] <- GEIA.emis.grid.val
+    } else {
+# '%i' draws error? and '%d'??
+#      cat(sprintf('%s: ERROR: attempting to write value from GEIA grid#=%i to grid index=%i (lon=%d,lat=%d) -> %f',
+      cat(sprintf('%s: ERROR: attempting to write value=%f from GEIA grid#=%.0f to global grid index=%.0f (lon=%.1f,lat=%.1f)\n',
+        this.fn, GEIA.emis.grid.val, GEIA.emis.grid.n, 
+        global.emis.vec.index, lon.lat.vec[1], lon.lat.vec[2]))
+      return # TODO: end execution of caller
+    }
+  } # end for loop over GEIA data values
+  global.emis.vec # return
+} # end function create.global.emissions.vector
+
+# TODO: refactor! this plot.layer is copy/mod from ioapi::plotLayersForTimestep.r::plot.layer
+# changes to arguments:
+# * attrs.list -> x.centers, y.centers (since these are in degrees, not km)
+# * map not passed: I'm using global map via `map(add=TRUE)` TODO: fixme
+plot.layer <- function(
+  data,             # data to plot (required)
+  title,            # string for plot title (required?)
+                    # TODO: handle when null!
+  subtitle=NULL,    # string for plot subtitle
+  x.centers,        # centers of abscissa points
+  y.centers,        # centers of ordinate points
+  q.vec=NULL,       # for quantiles
+  colors
+) {
+  library(fields)
+  library(maps)
+  if (sum(!is.na(data)) && (!is.null(q.vec))) {
+    plot.list <- list(x=x.centers, y=y.centers, z=data)
+    quantiles <- quantile(c(data), q.vec, na.rm=TRUE)
+    quantiles.formatted <- format(as.numeric(quantiles), digits=3)
+# start debugging
+#      print(paste('Non-null image.plot for source layer==', i.layer, ', quantile bounds=='))
+#      print(quantiles)
+#   end debugging
+    if (is.null(subtitle)) {
+      image.plot(plot.list, xlab="", ylab="", axes=F, col=colors(100),
+        axis.args=list(at=quantiles, labels=quantiles.formatted),
+        main=title)
+    } else {
+      image.plot(plot.list, xlab="", ylab="", axes=F, col=colors(100),
+        axis.args=list(at=quantiles, labels=quantiles.formatted),
+        main=title, sub=subtitle)
+    }
+    map(add=TRUE)
+  } else {
+# debugging
+#      print(paste('Null image.plot for source layer=', i.layer))
+    if (is.null(subtitle)) {
+      plot(0, type="n", axes=F, xlab="", ylab="",
+        xlim=range(x.centers), ylim=range(y.centers),
+        main=title)
+    } else {
+      plot(0, type="n", axes=F, xlab="", ylab="",
+        xlim=range(x.centers), ylim=range(y.centers),
+        main=title, sub=subtitle)
+    }
+    map(add=TRUE)
+  } # end testing data
+} # end function plot.layer
+
+# TODO: refactor! subtitle.stats is copied from ioapi::plotLayersForTimestep.r
+subtitle.stats <- function(vec) {
+  return.str <- ""
+  # is it numeric, and not empty?
+  if (is.numeric(vec) && sum(!is.na(vec))) {
+#    unsparse.vec <- subset(vec, !is.na(vec)) # fail: intended for interactive use
+#    unsparse.vec <- na.omit(vec) # fail: omits all *rows* containing an NA!
+    grids <- length(vec)
+    grids.str <- sprintf('(of cells=%i)', grids)
+    unsparse.vec <- vec[!is.na(vec)]
+    obs <- length(unsparse.vec)
+    obs.str <- sprintf('obs=%i', obs)
+    # use constants defined above. TODO: compute these once!
+    max.str <- sprintf(max.str, max(unsparse.vec))
+    med.str <- sprintf(med.str, median(unsparse.vec))
+    min.str <- sprintf(min.str, min(unsparse.vec))
+    return.str <-
+      sprintf('%s %s: %s, %s, %s',
+              obs.str, grids.str, min.str, med.str, max.str)
+  } else {
+    return.str <-"no data"
+  }
+  return.str
+} # end function subtitle.stats
+
+
+# ----------------------------------------------------------------------
+# payload
+# ----------------------------------------------------------------------
+
+# process input
+
+# columns are grid#, mass
+GEIA.emis.mx  <-
+  as.matrix(read.table(GEIA.emis.txt.fp, skip=GEIA.emis.txt.n.header))
+# mask zeros? no, use NA for non-ocean areas
+# GEIA.emis.mx[GEIA.emis.mx == 0] <- NA
+# <simple input check>
+dim(GEIA.emis.mx) ## [1] 36143     2
+# start debug
+GEIA.emis.mx.rows <- dim(GEIA.emis.mx)[1]
+if (GEIA.emis.mx.rows > GEIA.emis.grids.dim) {
+  cat(sprintf('%s: ERROR: GEIA.emis.mx.rows=%.0d > GEIA.emis.grids.dim=%.0d\n',
+    this.fn, GEIA.emis.mx.rows, GEIA.emis.grids.dim))
+} else {
+  cat(sprintf('%s: debug: GEIA.emis.mx.rows < GEIA.emis.grids.dim\n',
+    this.fn))
+}
+#   end debug
+# </simple input check>
+
+global.emis.vec <-
+  create.global.emissions.vector(
+    GEIA.emis.grids.dim, GEIA.emis.mx.rows, GEIA.emis.mx)
+
+# <visual input check>
+# Need sorted lat and lon vectors: we know what those are a priori
+# Add 0.5 since grid centers
+lon.vec <- 0.5 +
+  seq(from=grid.lon.degree.start, by=grid.lon.degree.per.cell, length.out=GEIA.emis.lon.dim)
+lat.vec <- 0.5 +
+  seq(from=grid.lat.degree.start, by=grid.lat.degree.per.cell, length.out=GEIA.emis.lat.dim)
+
+# Create emissions matrix corresponding to those dimensional vectors
+# (i.e., global.emis.mx is the "projection" of global.emis.vec)
+# First, create empty global.emis.mx? No, fill from global.emis.vec.
+# NO: I cannot just fill global.emis.mx from global.emis.vec:
+# latter's/GEIA's grid numbering system ensures 1000 lons per lat!
+# Which overflows the "real-spatial" global.emis.mx :-(
+# So I need to fill global.emis.mx using a for loop to decode the grid indices :-(
+# (but at least I can fill in whatever direction I want :-)
+global.emis.mx <- matrix(
+  rep(NA, GEIA.emis.grids.dim), nrow=GEIA.emis.lat.dim, ncol=GEIA.emis.lon.dim)
+
+# 1: works if subsequently transposed: TODO: FIXME
+for (i.lon in 1:GEIA.emis.lon.dim) {
+  for (i.lat in 1:GEIA.emis.lat.dim) {
+# 2: fails with 'dimensions of z are not length(x)(-1) times length(y)(-1)'
+# for (i.lat in 1:GEIA.emis.lat.dim) {
+#   for (i.lon in 1:GEIA.emis.lon.dim) {
+# 3: fails with 'dimensions of z are not length(x)(-1) times length(y)(-1)'
+# for (i.lon in GEIA.emis.lon.dim:1) {
+#   for (i.lat in GEIA.emis.lat.dim:1) {
+# 4: fails with 'dimensions of z are not length(x)(-1) times length(y)(-1)'
+# for (i.lat in GEIA.emis.lat.dim:1) {
+#   for (i.lon in GEIA.emis.lon.dim:1) {
+    lon <- lon.vec[i.lon]
+    lat <- lat.vec[i.lat]
+    GEIA.emis.grid.val <-
+      global.emis.vec[
+        lon.lat.vec.to.grid.index(c(lon, lat),
+          n.lon=GEIA.emis.lon.dim, n.lat=GEIA.emis.lat.dim)]
+    if (!is.na(GEIA.emis.grid.val)) {
+      if (is.na(global.emis.mx[i.lat, i.lon])) {
+        global.emis.mx[i.lat, i.lon] <- GEIA.emis.grid.val
+# start debug
+#        cat(sprintf(
+#          'debug: %s: writing val=%f to global.emis.mx[%.0f,%.0f] for grid center=[%f,%f]\n',
+#          this.fn, GEIA.emis.grid.val, i.lon, i.lat, lon, lat))
+#   end debug
+      } else {
+        # error if overwriting presumably-previously-written non-NA!
+        cat(sprintf(
+          '%s: ERROR: overwriting val=%f with val=%f at global.emis.mx[%.0f,%.0f] for grid center=[%f,%f]\n',
+          this.fn, global.emis.mx[i.lat, i.lon], GEIA.emis.grid.val,
+          i.lon, i.lat, lon, lat))
+        return # TODO: abend
+      } # end testing target != NA (thus not previously written)
+    } # end testing source != NA (don't write if is.na(lookup)
+  } # end for loop over lats
+} # end for loop over lons
+
+# start debugging
+# pdf(file=pdf.fp, width=5.5, height=4.25)
+# # 1: TODO: FIXME: why do I need to transpose global.emis.mx?
+# image(lon.vec, lat.vec, t(global.emis.mx))
+# # 2,3,4: how it should work ?!?
+# # image(lon.vec, lat.vec, global.emis.mx)
+# map(add=TRUE)
+# # dev.off() # leave on for subsequent plot; this will be p1/2
+#   end debugging
+# </visual input check>
+
+# write output to netCDF
+
+# create dimensions and dimensional variables
+time.vec <- c(0) # annual value, corresponding to no specific year
+time.dim <- ncdf4::ncdim_def(
+  name=time.var.name,
+  units=time.var.units,
+  vals=time.vec,
+  unlim=TRUE,
+  create_dimvar=TRUE,
+  calendar=time.var.calendar,
+  longname=time.var.long_name)
+
+lon.dim <- ncdf4::ncdim_def(
+  name=lon.var.name,
+  units=lon.var.units,
+  vals=lon.vec,
+  unlim=FALSE,
+  create_dimvar=TRUE,
+  longname=lon.var.long_name)
+
+lat.dim <- ncdf4::ncdim_def(
+  name=lat.var.name,
+  units=lat.var.units,
+  vals=lat.vec,
+  unlim=FALSE,
+  create_dimvar=TRUE,
+  longname=lat.var.long_name)
+
+# create data variable (as container--can't put data until we have a file)
+emis.var <- ncdf4::ncvar_def(
+  name=emis.var.name,
+  units=emis.var.units,
+#  dim=c(time.dim, lat.dim, lon.dim),
+#  dim=list(time.dim, lat.dim, lon.dim),
+# note dim order desired for result=var(time, lat, lon)
+  dim=list(lon.dim, lat.dim, time.dim),
+  missval=as.double(emis.var._FillValue),
+  longname=emis.var.long_name,
+  prec="double")
+
+# get current time for creation_date
+# system(intern=TRUE) -> return char vector, one member per output line)
+netcdf.timestamp <- system('date', intern=TRUE)
+
+# create netCDF file
+netcdf.file <- ncdf4::nc_create(
+  filename=netcdf.fp,
+#  vars=list(emis.var),
+#  verbose=TRUE)
+  vars=list(emis.var))
+
+# Write data to data variable: gotta have file first.
+# Gotta convert 2d global.emis.mx[lat,lon] to 3d global.emis.arr[time,lat,lon]
+# Do this before adding _FillValue to prevent:
+# > Error in R_nc4_put_vara_double: NetCDF: Not a valid data type or _FillValue type mismatch
+## global.emis.arr <- global.emis.mx
+## dim(global.emis.arr) <- c(1, dim(global.emis.mx))
+## global.emis.arr[1,,] <- global.emis.mx
+
+# Note
+# * global.emis.mx[lat,lon]
+# * datavar needs [lon, lat, time] (with time *last*)
+
+ncdf4::ncvar_put(
+  nc=netcdf.file,
+  varid=emis.var,
+#  vals=global.emis.arr,
+  vals=t(global.emis.mx),
+#  start=rep.int(1, length(dim(global.emis.arr))),
+  start=c(1, 1, 1),
+#  count=dim(global.emis.arr))
+  count=c(-1,-1, 1)) # -1 -> all data
+
+# Write netCDF attributes
+# Note: can't pass *.dim as varid, even though these are coordinate vars :-(
+
+# add datavar attributes
+ncdf4::ncatt_put(
+  nc=netcdf.file,
+#  varid=lon.var,
+  varid=lon.var.name,
+  attname="comment",
+  attval=lon.var.comment,
+  prec="text")
+
+ncdf4::ncatt_put(
+  nc=netcdf.file,
+#  varid=lat.var,
+  varid=lat.var.name,
+  attname="comment",
+  attval=lat.var.comment,
+  prec="text")
+
+# put _FillValue after putting data!
+ncdf4::ncatt_put(
+  nc=netcdf.file,
+  varid=emis.var,
+  attname="_FillValue",
+  attval=emis.var._FillValue,
+  prec="float") # why is "emi_n2o:missing_value = -999."?
+
+# add global attributes (varid=0)
+ncdf4::ncatt_put(
+  nc=netcdf.file,
+  varid=0,
+  attname="creation_date",
+  attval=netcdf.timestamp,
+  prec="text")
+
+ncdf4::ncatt_put(
+  nc=netcdf.file,
+  varid=0,
+  attname="source_file",
+  attval=netcdf.source_file,
+  prec="text")
+
+ncdf4::ncatt_put(
+  nc=netcdf.file,
+  varid=0,
+  attname="Conventions",
+  attval=netcdf.Conventions,
+  prec="text")
+
+# flush to file (there may not be data on disk before this point)
+# nc_sync(netcdf.file) # so we don't hafta reopen the file, below
+# Nope: per David W. Pierce Mon, 27 Aug 2012 21:35:35 -0700, ncsync is not enough
+ncdf4::nc_close(netcdf.file)
+
+# NOTE: you must assign when you nc_open! (though not when you nc_close)
+# If you do
+## nc_open(netcdf.fp,
+##         write=FALSE,    # will only read below
+##         readunlim=TRUE) # it's a small file
+# you will get this error on your first attempt to ncvar_get
+# (though no error on ncvar_put, annoyingly enough):
+# > Error in if (nc$var[[li]]$hasAddOffset) addOffset = nc$var[[li]]$addOffset else addOffset = 0 :
+# > argument is of length zero
+# MERCI BEAUCOUP, Pascal Oettli!
+netcdf.file <- ncdf4::nc_open(
+  filename=netcdf.fp,
+  write=FALSE,    # will only read below
+  readunlim=TRUE) # it's a small file
+
+# <simple output check>
+system(sprintf('ls -alth %s', netcdf.fp))
+system(sprintf('ncdump -h %s', netcdf.fp))
+# TODO: call sample stats from netCDF.stats.to.stdout.r via function, not like this!
+# system('Rscript ./netCDF.stats.to.stdout.r netcdf.fp=./GEIA_N2O_oceanic.nc var.name=emi_n2o')
+source(stat.script.fp)
+netCDF.stats.to.stdout(
+  netcdf.fp='./GEIA_N2O_oceanic.nc',
+  var.name='emi_n2o')
+# </simple output check>
+
+# <visual output check>
+# TODO: do plot-related refactoring! allow to work with projects={ioapi, this}
+# <copied from plotLayersForTimestep.r>
+# double-sprintf-ing to set precision by constant: cool or brittle?
+stats.precision <- 3 # sigdigs to use for min, median, max of obs
+stat.str <- sprintf('%%.%ig', stats.precision)
+# use these in function=subtitle.stats as sprintf inputs
+max.str <- sprintf('max=%s', stat.str)
+med.str <- sprintf('med=%s', stat.str)
+min.str <- sprintf('min=%s', stat.str)
+# </copied from plotLayersForTimestep.r>
+
+# Get the data out of the datavar, to test reusability
+target.data <- ncdf4::ncvar_get(
+  nc=netcdf.file,
+  varid=emis.var,
+  # read all the data
+  start=rep(1, emis.var$ndims),
+  count=rep(-1, emis.var$ndims)) 
+
+# <simple output check/>
+dim(target.data) # n.lon, n.lat
+
+# <copy/mod from windowEmissions.r>
+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)
+probabilities.vec <- seq(0, 1, 1.0/(length(palette.vec) - 1))
+# </copy/mod from windowEmissions.r>
+
+pdf(file=pdf.fp, width=5.5, height=4.25)
+# <copy/mod from plotLayersForTimestep.r>
+plot.layer(target.data,
+  title=netcdf.title,
+  subtitle=subtitle.stats(target.data),
+  x.centers=lon.vec,
+  y.centers=lat.vec,
+  q.vec=probabilities.vec,
+  colors=colors)
+# </copy/mod from plotLayersForTimestep.r>
+library(maps)
+maps::map(add=TRUE)
+# </visual output check>
+
+# teardown
+dev.off()
+ncdf4::nc_close(netcdf.file)
+quit(save="no") # just exit

File GEIA_to_netCDF.sh

-#!/usr/bin/env bash
-# Driver for GEIA.to.netCDF.r--configure as needed for your platform.
-# NOTE regarding R and source-ing below!
-
-# 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
-# for R on EMVL: don't ask :-(
-R_DIR='/usr/local/R-2.15.0/bin'
-R="${R_DIR}/R"
-RSCRIPT="${R_DIR}/Rscript"
-
-export TEST_DIR='.' # keep it simple for now: same dir as top of git repo
-mkdir -p ${TEST_DIR}
-
-# 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}"
-PDF_FN="test_$(date +${DATE_FORMAT}).pdf"
-export PDF_FP="${TEST_DIR}/${PDF_FN}" # path to PDF output
-
-# this does the converting and plotting
-CONVERT_PLOT_R="${TEST_DIR}/GEIA.to.netCDF.r" 
-
-# get GEIA input text from my repository (stable location, no unzip)
-GEIA_INPUT_URI='https://github.com/downloads/TomRoche/GEIA_to_netCDF/N2OOC90Y.1A'
-GEIA_INPUT_FN="$(basename ${GEIA_INPUT_URI})"
-export GEIA_INPUT_FP="${TEST_DIR}/${GEIA_INPUT_FN}"
-if [[ ! -r "${GEIA_INPUT_FP}" ]] ; then
-  for CMD in \
-    "wget --no-check-certificate -c -O ${GEIA_INPUT_FP} ${GEIA_INPUT_URI}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# If this fails ...
-"${RSCRIPT}" "${CONVERT_PLOT_R}"
-# ... just start R ...
-# "${R}"
-# ... and source ${CONVERT_PLOT_R}, e.g.,
-# source('./GEIA.to.netCDF.r')
-
-# After exiting R, show output files and display output PDF.
-for CMD in \
-  "ls -alht ${TEST_DIR}" \
-  "${PDF_VIEWER} ${PDF_FP}" \
-; do
-  echo -e "$ ${CMD}"
-  eval "${CMD}"
-done

File global_to_AQMEII.r

+# R code to 2D-regrid GEIA global/unprojected netCDF data to AQMEII-NA, an LCC-projected subdomain.
+# If running manually in R console, remember to run setup actions: `source ./regrid_GEIA_netCDF.sh`
+
+## library(ncdf4)
+## library(fields)
+## # run if rgeos not available
+## # gpclibPermit()
+## 
+
+# ----------------------------------------------------------------------
+# constants
+# ----------------------------------------------------------------------
+
+this.fn <- Sys.getenv('CALL_REGRID_FN')
+
+# all the following env vars must be set and exported in driver script
+work.dir <- Sys.getenv('WORK_DIR')
+pdf.fp <- Sys.getenv('GEIA_REGRID_PDF_FP')
+pdf.er <- Sys.getenv('PDF_VIEWER')
+in.fp <- Sys.getenv('GEIA_REFORMAT_FP')
+in.band <- Sys.getenv('GEIA_REFORMAT_BAND')
+data.var.name <- Sys.getenv('GEIA_REFORMAT_DATAVAR_NAME')
+data.var.longname <- Sys.getenv('GEIA_REFORMAT_DATAVAR_LONGNAME')
+data.var.unit <- Sys.getenv('GEIA_REFORMAT_DATAVAR_UNIT')
+data.var.na <- as.numeric(Sys.getenv('GEIA_REFORMAT_DATAVAR_NA')) # must convert from string
+template.var.name <- Sys.getenv('TEMPLATE_VAR_NAME')
+template.in.fp <- Sys.getenv('TEMPLATE_INPUT_FP')
+template.band <- Sys.getenv('TEMPLATE_INPUT_BAND')
+out.fp <- Sys.getenv('GEIA_REGRID_FP')
+x.var.name <- Sys.getenv('GEIA_REGRID_X_VAR_NAME')
+y.var.name <- Sys.getenv('GEIA_REGRID_Y_VAR_NAME')
+
+# ----------------------------------------------------------------------
+# functions
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# payload
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# setup
+# ----------------------------------------------------------------------
+
+# 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
+
+stat.script.fp <- Sys.getenv('STAT_SCRIPT_FP')
+source(stat.script.fp) # produces errant error=
+#> netCDF.stats.to.stdout.r: no arguments supplied, exiting
+plot.script.fp <- Sys.getenv('PLOT_SCRIPT_FP')
+source(plot.script.fp)
+
+# if any of the above failed, check how you ran (e.g., `source`d) driver script
+# plot constants below
+
+# check input
+
+system(sprintf('ncdump -h %s', in.fp))
+# netcdf GEIA_N2O_oceanic {
+# dimensions:
+#   lon = 360 ;
+#   lat = 180 ;
+#   time = UNLIMITED ; // (1 currently)
+# variables:
+#   double lon(lon) ;
+# ...
+#   double lat(lat) ;
+# ...
+#   double time(time) ;
+# ...
+#   double emi_n2o(time, lat, lon) ;
+#     emi_n2o:units = "ton N2O-N/yr" ;
+#     emi_n2o:missing_value = -999. ;
+#     emi_n2o:long_name = "N2O emissions" ;
+#     emi_n2o:_FillValue = -999.f ;
+
+# Note missing/fill values: more below.
+
+netCDF.stats.to.stdout(netcdf.fp=in.fp, var.name=data.var.name)
+# For /tmp/projectRasterTest/GEIA_N2O_oceanic.nc var=emi_n2o
+#       cells=64800
+#       obs=36143
+#       min=5.96e-08
+#       q1=30.4
+#       med=67.7
+#       mean=99.5
+#       q3=140
+#       max=1.17e+03
+# Note min and max of the input, and that min > 0.
+
+# make input raster
+
+library(raster)
+in.raster <- raster::raster(in.fp, varname=data.var.name, band=in.band)
+in.raster
+# class       : RasterLayer 
+# dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
+# resolution  : 1, 1  (x, y)
+# extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
+# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
+# values      : /tmp/projectRasterTest/GEIA_N2O_oceanic.nc 
+# layer name  : N2O.emissions 
+# z-value     : 0 
+# zvar        : emi_n2o 
+
+# make template raster to format output.
+# Note below that one should be able to projectRaster(...) with a resolution (i.e.,
+# from=in.raster, res=c(12e3, 12e3), crs=out.proj4,
+# ) as well as a template (i.e.,
+# from=in.raster, to=template.raster, crs=out.proj4,
+# ). So why create a template? Because just giving a resolution hangs.
+
+# look @ input (one of the PHASE AQ inputs) to template/extents raster:
+# truncate 'ncdump' to avoid IOAPI cruft
+system(sprintf('ncdump -h %s | head -n 13', template.in.fp))
+# netcdf emis_mole_all_20080101_12US1_cmaq_cb05_soa_2008ab_08c.EXTENTS_INPUT {
+# dimensions:
+#   TSTEP = UNLIMITED ; // (25 currently)
+#   LAY = 1 ;
+#   ROW = 299 ;
+#   COL = 459 ;
+# variables:
+#   float emi_n2o(TSTEP, LAY, ROW, COL) ;
+#     emi_n2o:long_name = "XYL             " ;
+#     emi_n2o:units = "moles/s         " ;
+#     emi_n2o:var_desc = "Model species XYL                                                               " ;
+
+# 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
+
+template.in.raster <-
+  raster::raster(template.in.fp, varname=template.var.name, band=template.band)
+template.raster <- raster::projectExtent(template.in.raster, crs=out.proj4)
+#> 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
+template.raster@extent <- template.extents
+# should resemble the domain specification @
+# https://github.com/TomRoche/cornbeltN2O/wiki/AQMEII-North-American-domain#wiki-EPA
+template.raster
+# 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 
+
+# start debug---------------------------------------------------------
+# regrid.start.time <- system('date', intern=TRUE)
+# regrid.start.str <- sprintf('start regrid @ %s', regrid.start.time)
+# cat(sprintf('%s\n', regrid.start.str))
+#   end debug---------------------------------------------------------
+
+# 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.proj4,
+    # 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=data.var.na,
+    varname=data.var.name, 
+    varunit=data.var.unit,
+    longname=data.var.longname,
+    xname=x.var.name,
+    yname=y.var.name,
+    filename=out.fp)
+# above fails to set CRS, so
+library(sp)
+out.raster@crs <- sp::CRS(out.proj4)
+out.raster
+# 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 
+# data source : /tmp/projectRasterTest/GEIA_N2O_oceanic_regrid.nc 
+# names       : N2O.emissions 
+# zvar        : emi_n2o 
+
+# start debug---------------------------------------------------------
+# regrid.end.time <- system('date', intern=TRUE)
+# cat(sprintf('  end regrid @ %s\n', regrid.end.time))
+# cat(sprintf('%s\n', regrid.start.time))
+#   end debug---------------------------------------------------------
+
+system(sprintf('ls -alht %s', work.dir))
+system(sprintf('ncdump -h %s', out.fp))
+# netcdf GEIA_N2O_oceanic_regrid {
+# dimensions:
+#   COL = 459 ;
+#   ROW = 299 ;
+# variables:
+#   double COL(COL) ;
+#     COL:units = "meter" ;
+#     COL:long_name = "COL" ;
+#   double ROW(ROW) ;
+#     ROW:units = "meter" ;
+#     ROW:long_name = "ROW" ;
+#   float emi_n2o(ROW, COL) ;
+#     emi_n2o:units = "ton N2O-N/yr" ;
+#     emi_n2o:_FillValue = -999. ;
+#     emi_n2o:missing_value = -999. ;
+#     emi_n2o:long_name = "N2O emissions" ;
+#     emi_n2o:projection = "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000" ;
+#     emi_n2o:projection_format = "PROJ.4" ;
+#     emi_n2o:min = 0.83728 ;
+#     emi_n2o:max = 522.693774638276 ;
+# ...
+
+netCDF.stats.to.stdout(netcdf.fp=out.fp, var.name=data.var.name)
+# For ./GEIA_N2O_oceanic_regrid.nc var=emi_n2o
+#   cells=137241
+#   obs=46473
+#   min=0.837
+#   q1=11.9
+#   med=54.1
+#   mean=71.9
+#   q3=85.5
+#   max=523
+
+# plot to PDF
+
+# plot constants
+
+library(maptools)
+data(wrld_simpl) # from maptools
+# map.us.unproj <- wrld_simpl[wrld_simpl$ISO3 == 'USA', ]  # unprojected
+library(rgdal)
+# map.us.proj <- rgdal::spTransform(map.us.unproj, sp::CRS(out.proj4))  # projected
+# do North America instead
+map.us.unproj <- wrld_simpl[wrld_simpl$ISO3 %in% c('CAN', 'MEX', 'USA'),]
+map.us.proj <-
+  rgdal::spTransform(map.us.unproj, sp::CRS(out.proj4)) # projected
+
+title <- sprintf('%s\n%s',
+  Sys.getenv('GEIA_REGRID_PDF_TITLE'), Sys.getenv('GEIA_REGRID_DATAVAR_UNIT'))
+
+out.data.vec <- getValues(out.raster) # raster data as vector
+subtitle <- subtitle.stats(out.data.vec)
+
+# package=M3 map for fields::image.plot
+map.table.fp <- 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
+
+x.centers <- raster.centers.x(out.raster)
+y.centers <- raster.centers.y(out.raster)
+quantiles <- quantile(out.data.vec, probabilities.vec, na.rm=TRUE)
+quantiles.formatted <- format(as.numeric(quantiles), digits=3)
+
+# repeat for each plot file
+pdf(file=pdf.fp, width=5.5, height=4.25)
+
+# plot page 1: raster::plot-------------------------------------------
+# plot(out.raster, main=title, sub=subtitle) # works
+plot(out.raster, # remaining args from image.plot
+  main=title, sub=subtitle,
+  xlab='', ylab='', axes=F, col=colors(100),
+  axis.args=list(at=quantiles, labels=quantiles.formatted))
+# add a projected CONUS map
+plot(map.us.proj, add=TRUE)
+
+# plot page 2: fields::image.plot-------------------------------------
+
+plot.raster(
+  raster=out.raster,
+  title=title, 
+  subtitle=subtitle,
+  q.vec=probabilities.vec,
+  colors,
+  map.cmaq
+)
+
+# step through plot.raster:
+# package=fields needs data as matrix, not vector
+# out.data.mx <- out.data.vec
+# dim(out.data.mx) <- c(length(x.centers), length(y.centers)) # cols, rows
+# plot.data(out.data.mx, title, subtitle, x.centers, y.centers, probabilities.vec, colors, map.cmaq)
+
+# step through plot.data
+# plot.list <- list(x=x.centers, y=y.centers, z=out.data.mx)
+# image.plot(plot.list, xlab='', ylab='', axes=F, col=colors(100),
+#   axis.args=list(at=quantiles, labels=quantiles.formatted),
+#   main=title, sub=subtitle)
+# lines(map.cmaq)
+
+#   end image.plot----------------------------------------------------
+
+# flush the plot device
+dev.off()

File regrid.global.to.AQMEII.r

-# R code to 2D-regrid GEIA global/unprojected netCDF data to AQMEII-NA, an LCC-projected subdomain.
-# If running manually in R console, remember to run setup actions: `source ./regrid_GEIA_netCDF.sh`
-
-library(ncdf4)
-library(fields)
-library(maptools)
-# run if rgeos not available
-# gpclibPermit()
-library(rgdal)
-library(raster)
-library(M3) # for extents calculation
-data(wrld_simpl) # from maptools
-
-# constants-----------------------------------------------------------
-
-this.fn <- 'regrid.global.to.AQMEII.r'  # TODO: get like $0
-
-# all the following env vars must be set and exported in driver script
-test.dir <- Sys.getenv('TEST_DIR')
-pdf.dir <- Sys.getenv('PDF_DIR')
-pdf.fp <- Sys.getenv('PDF_FP')
-pdf.er <- Sys.getenv('PDF_VIEWER')
-in.fp <- Sys.getenv('DATA_INPUT_FP')
-in.band <- Sys.getenv('DATA_INPUT_BAND')
-data.var.name <- Sys.getenv('DATA_VAR_NAME')
-data.var.longname <- Sys.getenv('DATA_VAR_LONGNAME')
-data.var.unit <- Sys.getenv('DATA_VAR_UNIT')
-data.var.na <- as.numeric(Sys.getenv('DATA_VAR_NA')) # must convert from string
-template.var.name <- Sys.getenv('TEMPLATE_VAR_NAME')
-template.in.fp <- Sys.getenv('TEMPLATE_INPUT_FP')
-template.band <- Sys.getenv('TEMPLATE_INPUT_BAND')
-out.fp <- Sys.getenv('DATA_OUTPUT_FP')
-
-# coordinate reference system:
-# use package=M3 to get CRS from template file
-out.crs <- get.proj.info.M3(template.in.fp)
-cat(sprintf('out.crs=%s\n', out.crs)) # debugging
-# out.crs=+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000
-
-stat.script.fp <- Sys.getenv('STAT_SCRIPT_FP')
-source(stat.script.fp) # produces errant error=
-#> netCDF.stats.to.stdout.r: no arguments supplied, exiting
-plot.script.fp <- Sys.getenv('PLOT_SCRIPT_FP')
-source(plot.script.fp)
-
-# if any of the above failed, check how you ran (e.g., `source`d) driver script
-# plot constants below
-
-# payload-------------------------------------------------------------
-
-# check input
-
-system(sprintf('ncdump -h %s', in.fp))
-# netcdf GEIA_N2O_oceanic {
-# dimensions:
-#   lon = 360 ;
-#   lat = 180 ;
-#   time = UNLIMITED ; // (1 currently)
-# variables:
-#   double lon(lon) ;
-# ...
-#   double lat(lat) ;
-# ...
-#   double time(time) ;
-# ...
-#   double emi_n2o(time, lat, lon) ;
-#     emi_n2o:units = "ton N2O-N/yr" ;
-#     emi_n2o:missing_value = -999. ;
-#     emi_n2o:long_name = "N2O emissions" ;
-#     emi_n2o:_FillValue = -999.f ;
-
-# Note missing/fill values: more below.
-
-netCDF.stats.to.stdout(netcdf.fp=in.fp, var.name=data.var.name)
-# For /tmp/projectRasterTest/GEIA_N2O_oceanic.nc var=emi_n2o
-#       cells=64800
-#       obs=36143
-#       min=5.96e-08
-#       q1=30.4
-#       med=67.7
-#       mean=99.5
-#       q3=140
-#       max=1.17e+03
-# Note min and max of the input, and that min > 0.
-
-# make input raster
-
-in.raster <- raster(in.fp, varname=data.var.name, band=in.band)
-in.raster
-# class       : RasterLayer 
-# dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
-# resolution  : 1, 1  (x, y)
-# extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
-# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
-# values      : /tmp/projectRasterTest/GEIA_N2O_oceanic.nc 
-# layer name  : N2O.emissions 
-# z-value     : 0 
-# zvar        : emi_n2o 
-
-# make template raster to format output.
-# Note below that one should be able to projectRaster(...) with a resolution (i.e.,
-# from=in.raster, res=c(12e3, 12e3), crs=out.crs,
-# ) as well as a template (i.e.,
-# from=in.raster, to=template.raster, crs=out.crs,
-# ). So why create a template? Because just giving a resolution hangs.
-
-# look @ input (one of the PHASE AQ inputs) to template/extents raster:
-# truncate 'ncdump' to avoid IOAPI cruft
-system(sprintf('ncdump -h %s | head -n 13', template.in.fp))
-# netcdf emis_mole_all_20080101_12US1_cmaq_cb05_soa_2008ab_08c.EXTENTS_INPUT {
-# dimensions:
-#   TSTEP = UNLIMITED ; // (25 currently)
-#   LAY = 1 ;
-#   ROW = 299 ;
-#   COL = 459 ;
-# variables:
-#   float emi_n2o(TSTEP, LAY, ROW, COL) ;
-#     emi_n2o:long_name = "XYL             " ;
-#     emi_n2o:units = "moles/s         " ;
-#     emi_n2o:var_desc = "Model species XYL                                                               " ;
-
-# use M3 to get extents from template file (thanks CGN!)
-extents.info <- get.grid.info.M3(template.in.fp)
-extents.xmin <- extents.info$x.orig
-extents.xmax <- max(
-  get.coord.for.dimension(
-    file=template.in.fp, dimension="col", position="upper", units="m")$coords)
-extents.ymin <- extents.info$y.orig
-extents.ymax <- max(
-  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 <-
-  extent(extents.xmin, extents.xmax, extents.ymin, extents.ymax)
-template.extents
-
-template.in.raster <- raster(template.in.fp, varname=template.var.name, band=template.band)
-template.raster <- projectExtent(template.in.raster, crs=out.crs)
-#> Warning message:
-#> In projectExtent(template.in.raster, out.crs) :
-#>   158 projected point(s) not finite
-# is that "projected point(s) not finite" warning important? Probably not, per Hijmans
-template.raster@extent <- template.extents
-# should resemble the domain specification @
-# https://github.com/TomRoche/cornbeltN2O/wiki/AQMEII-North-American-domain#wiki-EPA
-template.raster
-# 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 
-
-# start debug---------------------------------------------------------
-# regrid.start.time <- system('date', intern=TRUE)
-# regrid.start.str <- sprintf('start regrid @ %s', regrid.start.time)
-# cat(sprintf('%s\n', regrid.start.str))
-#   end debug---------------------------------------------------------
-
-# at last: do the regridding
-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
-    NAflag=data.var.na,
-    varname=data.var.name, 
-    varunit=data.var.unit,
-    longname=data.var.longname,
-    xname='COL',
-    yname='ROW',
-    filename=out.fp)
-# above fails to set CRS, so
-out.raster@crs <- CRS(out.crs)
-out.raster
-# 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 
-# data source : /tmp/projectRasterTest/GEIA_N2O_oceanic_regrid.nc 
-# names       : N2O.emissions 
-# zvar        : emi_n2o 
-
-# start debug---------------------------------------------------------
-# regrid.end.time <- system('date', intern=TRUE)
-# cat(sprintf('  end regrid @ %s\n', regrid.end.time))
-# cat(sprintf('%s\n', regrid.start.time))
-#   end debug---------------------------------------------------------
-
-system(sprintf('ls -alht %s', test.dir))
-system(sprintf('ncdump -h %s', out.fp))
-# netcdf GEIA_N2O_oceanic_regrid {
-# dimensions:
-#   COL = 459 ;
-#   ROW = 299 ;
-# variables:
-#   double COL(COL) ;
-#     COL:units = "meter" ;
-#     COL:long_name = "COL" ;
-#   double ROW(ROW) ;
-#     ROW:units = "meter" ;
-#     ROW:long_name = "ROW" ;
-#   float emi_n2o(ROW, COL) ;
-#     emi_n2o:units = "ton N2O-N/yr" ;
-#     emi_n2o:_FillValue = -999. ;
-#     emi_n2o:missing_value = -999. ;
-#     emi_n2o:long_name = "N2O emissions" ;
-#     emi_n2o:projection = "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000" ;
-#     emi_n2o:projection_format = "PROJ.4" ;
-#     emi_n2o:min = 0.83728 ;
-#     emi_n2o:max = 522.693774638276 ;
-# ...
-
-netCDF.stats.to.stdout(netcdf.fp=out.fp, var.name=data.var.name)
-# For ./GEIA_N2O_oceanic_regrid.nc var=emi_n2o
-#   cells=137241
-#   obs=46473
-#   min=0.837
-#   q1=11.9
-#   med=54.1
-#   mean=71.9
-#   q3=85.5
-#   max=523
-
-# plot to PDF---------------------------------------------------------
-
-# plot constants
-
-# maps from package=raster
-# map.us.unproj <- wrld_simpl[wrld_simpl$ISO3 == 'USA', ]  # unprojected
-# map.us.proj <- spTransform(map.us.unproj, CRS(out.crs))  # projected
-# do North America instead
-map.us.unproj <- wrld_simpl[wrld_simpl$ISO3 %in% c('CAN', 'MEX', 'USA'),]
-map.us.proj <-
-  spTransform(map.us.unproj, CRS(out.crs)) # projected
-
-title <- 'GEIA annual oceanic N2O regridded to AQMEII-NA'
-units <- '(ton N2O-N/yr)' # TODO: get this from netCDF file
-# combine them. TODO: how to make units smaller, italicized?
-title <- sprintf('%s\n%s', title, units)
-
-out.data.vec <- getValues(out.raster) # raster data as vector
-subtitle <- subtitle.stats(out.data.vec)
-
-# package=M3 map for fields::image.plot
-map.table.fp <- 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
-
-x.centers <- raster.centers.x(out.raster)
-y.centers <- raster.centers.y(out.raster)
-quantiles <- quantile(out.data.vec, probabilities.vec, na.rm=TRUE)
-quantiles.formatted <- format(as.numeric(quantiles), digits=3)
-
-# repeat for each plot file
-pdf(file=pdf.fp, width=5.5, height=4.25)
-
-# plot page 1: raster::plot-------------------------------------------
-# plot(out.raster, main=title, sub=subtitle) # works
-plot(out.raster, # remaining args from image.plot
-  main=title, sub=subtitle,
-  xlab='', ylab='', axes=F, col=colors(100),
-  axis.args=list(at=quantiles, labels=quantiles.formatted))
-# add a projected CONUS map
-plot(map.us.proj, add=TRUE)
-
-# plot page 2: fields::image.plot-------------------------------------
-
-plot.raster(
-  raster=out.raster,
-  title=title, 
-  subtitle=subtitle,
-  q.vec=probabilities.vec,
-  colors,
-  map.cmaq
-)
-
-# step through plot.raster:
-# package=fields needs data as matrix, not vector
-# out.data.mx <- out.data.vec
-# dim(out.data.mx) <- c(length(x.centers), length(y.centers)) # cols, rows
-# plot.data(out.data.mx, title, subtitle, x.centers, y.centers, probabilities.vec, colors, map.cmaq)
-
-# step through plot.data
-# plot.list <- list(x=x.centers, y=y.centers, z=out.data.mx)
-# image.plot(plot.list, xlab='', ylab='', axes=F, col=colors(100),
-#   axis.args=list(at=quantiles, labels=quantiles.formatted),
-#   main=title, sub=subtitle)
-# lines(map.cmaq)
-
-#   end image.plot----------------------------------------------------
-
-# flush the plot device
-dev.off()

File regrid_GEIA_netCDF.sh

-#!/usr/bin/env bash
-# Driver for regrid.global.to.AQMEII.r--configure as needed for your platform.
-
-# start setup code copied from GEIA_to_netCDF.sh----------------------
-# TODO: refactor
-
-# 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
-# for R on EMVL: don't ask :-(
-R_DIR='/usr/local/R-2.15.0/bin'
-R="${R_DIR}/R"
-RSCRIPT="${R_DIR}/Rscript"
-
-export TEST_DIR='.' # keep it simple for now: same dir as top of git repo
-mkdir -p ${TEST_DIR}
-