Source

MOZART-global to AQMEII-NA / levelplot.netCDF.global.r

Full commit
# Code to rasterVis::levelplot global data in a netCDF file.
# Does not manipulate the data, just displays.
# TODO: make generic, pass parameters

library(ncdf4)        # on clusters only
library(maps)         # on tlrPanP5 as well as clusters

# library(latticeExtra) # on tlrPanP5 as well as clusters
# library(raster)
library(rasterVis) # loads both

# library(M3)           # for extents calculation
library(maptools)     # on tlrPanP5 as well as clusters
data(wrld_simpl)      # from maptools

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

# all the following env vars must be set and exported in driver script (TODO: fail if not)
this.fn <- Sys.getenv('CALL_R_FN') 
test.dir <- Sys.getenv('TEST_DIR')
pdf.dir <- Sys.getenv('PDF_DIR')
pdf.fp <- Sys.getenv('PDF_FP')
pdf.er <- Sys.getenv('PDF_VIEWER')
sigdigs <- as.numeric(Sys.getenv('SIGDIGS')) # significant digits for lattice strips
in.fp <- Sys.getenv('INPUT_DATA_FP')
#in.band <- Sys.getenv('DATA_INPUT_BAND')
in.band <- as.numeric(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

# coordinate reference system: here, just global lon-lat
global.proj <- CRS('+proj=longlat +ellps=WGS84')

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)

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

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

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

# 1. Load the input data

system(sprintf('ncdump -h %s', in.fp))
# netcdf \2008N2O_restart {
# dimensions:
#         lev = 56 ;
#         lat = 96 ;
#         lon = 144 ;
#         ilev = 57 ;
# variables:
#         double N2O(lev, lat, lon) ;
#                 N2O:long_name = "N2O" ;
#                 N2O:units = "ppmV" ;
#         double lon(lon) ;
#                 lon:long_name = "longitude" ;
#                 lon:units = "degrees_east" ;
#         double lat(lat) ;
#                 lat:long_name = "latitude" ;
#                 lat:units = "degrees_north" ;
#         double lev(lev) ;
#                 lev:long_name = "hybrid level at layer midpoints (1000*(A+B))" ;
#                 lev:units = "hybrid_sigma_pressure" ;

netCDF.stats.to.stdout(netcdf.fp=in.fp, var.name=data.var.name)
# For /home/rtd/code/regridding/MOZART_global_to_AQMEII-NA/2008N2O_restart.nc var=N2O
#         cells=774144
#         obs=774144
#         min=0.0146
#         q1=0.224
#         med=0.483
#         mean=0.366
#         q3=0.485
#         max=0.502
# Note min and max of the input, and that min > 0.

# make global input RasterBrick. note this encapsulates ncdf4 calls!
global.raster <- brick(in.fp, varname=data.var.name, crs=global.proj)
global.raster
# 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 +ellps=WGS84 +towgs84=0,0,0 
# data source : /home/rtd/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 
# get number of layers, for `levelplot`
layers.n <- nbands(global.raster)

# make input longitudes NOT zero-based (which causes problems for wrld_simpl, and probably other data)
global.raster <- rotate(global.raster,
                        overwrite=TRUE) # else levelplot does one layer per page?
global.raster
# 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 +ellps=WGS84 +towgs84=0,0,0 
# data source : in memory
# 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, ... 
# min values  :        0.01459371,        0.01522770,        0.01654624,        0.01816243,        0.01947412,        0.02304307,        0.02800324,       0.03213281,        0.03435659,        0.04158530,        0.04842525,        0.05333659,       0.06187457,        0.07327773,        0.09040873, ... 
# max values  :        0.06437672,        0.08462274,        0.10636217,        0.12670186,        0.14770647,        0.16961053,        0.19968871,       0.22470128,        0.23804380,        0.26975339,        0.29912971,        0.32279733,       0.34366432,        0.37546554,        0.40817063, ... 
# hybrid_sigma_pressure: 1.86787999700755, 992.500010610456 (min, max)

# get a data.frame for use with `levelplot`
# TODO: if using rasterVis::levelplot, get names with `layerNames`
global.df <- as.data.frame(global.raster)
# start debugging
#head(names(global.df))
#tail(names(global.df))
#   end debugging
#library('stringr')
#names(global.df) <- formatC(as.numeric(str_replace(names(global.df), "X", fixed(""))), format="g", digits=4)
# just use package=base for now
names(global.df) <- formatC(as.numeric(sub(x=names(global.df), pattern="X", replacement="", fixed=TRUE)), format="g", digits=sigdigs)
# debugging
head(names(global.df))
tail(names(global.df))
summary(unlist(global.df))
# [1] "1.868" "2.353" "2.948" "3.677" "4.562" "5.632"
# [1] "917.5" "932.5" "947.5" "962.5" "977.5" "992.5"
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
# 0.01459 0.22440 0.48350 0.36550 0.48530 0.50210

# get a global map
global.map <- map('world', plot=FALSE)
global.map.shp <- map2SpatialLines(global.map, proj4string=global.proj)

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

# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
#  layers,
  margin=FALSE,
#  FUN.margin=mean,
#  maxpixels=1e5,
#  par.settings=rasterTheme,
#  between=list(x=0.5, y=0.2),
#  as.table=TRUE,
#  xlab='', ylab='', main='', # arguments for levelplot
  names.attr=names(global.df),
#  scales=list(), scales.margin=NULL,
#  xscale.components=xscale.raster,
#  yscale.components=yscale.raster,
#  zscaleLog=NULL,
#  colorkey=list(space='right'),
#  panel=panel.levelplot,
#  contour=FALSE, region=TRUE, labels=FALSE,
  layout=c(1,layers.n), # all layers in one "atmospheric" column
#  strip.left=strip.custom(  # ... draw strip to left of panel
#  strip.right=strip.custom(  # ... draw strip to right of panel
#    factor.levels=names(global.df),
#    strip.levels=TRUE,
#    horizontal=TRUE,
#    strip.names=FALSE,
#    # gotta shrink strip text size to fit strip width: more on that separately
#    par.strip.text=list(cex=0.5)
#  )
# ...
) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))

# teardown
dev.off()
# if running from console
# quit(save="no") # just exit