Source

MOZART-global to AQMEII-NA / 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)