MOZART-global to AQMEII-NA / regrid_surface_elevation__MOZART.r

# R code to regrid pre-cropped surface elevation netCDF data to match other data over the same extents (but diferently gridded).
# TODO: make generic, pass parameters

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

cat('Setting global constants ...\n') # debugging

# global constants, mostly from environment---------------------------

cat('Loading package=raster ...\n') # debugging
library(raster)
# coordinate reference system: here, just global lon-lat
global.proj <- CRS('+proj=longlat +ellps=WGS84')

# all the following env vars must be set and exported in driver script (TODO: fail if not)
this.fn <- Sys.getenv('CALL_R_FN')  # FIXME: gets overidden by `source` below
cat(sprintf('this.fn=%s\n', this.fn)) # debugging
work.dir <- Sys.getenv('WORK_DIR')
cat(sprintf('work.dir=%s\n', work.dir)) # debugging

sigdigs <- as.numeric(Sys.getenv('SIGDIGS')) # significant digits for lattice strips

stat.script.fp <- Sys.getenv('STAT_SCRIPT_FP')
source(stat.script.fp) # produces errant error=
#> netCDF.stats.to.stdout.r: no arguments supplied, exiting
# unnecessary here
# plot.script.fp <- Sys.getenv('PLOT_SCRIPT_FP')
# source(plot.script.fp)

# local constants, from environment-----------------------------------
# taking from Rscript call failed, annoyingly and mysteriously :-(

cat('Setting local constants ...\n') # debugging

in.fp <- Sys.getenv('IN_FP')
in.datavar.name <- Sys.getenv('IN_DATAVAR_NAME')
in.datavar.band <- Sys.getenv('IN_DATAVAR_BAND')
in.crs <- global.proj # TODO: get from file!

template.fp <- Sys.getenv('TEMPLATE_FP')
template.datavar.name <- Sys.getenv('TEMPLATE_DATAVAR_NAME')

raster.brick <- Sys.getenv('RASTER_BRICK')
raster.rotate <- Sys.getenv('RASTER_ROTATE')
raster.plot <- Sys.getenv('RASTER_PLOT')

out.fp <- Sys.getenv('OUT_FP')
out.crs <- in.crs
out.datavar.name <- Sys.getenv('OUT_DATAVAR_NAME')
out.datavar.type <- Sys.getenv('OUT_DATAVAR_TYPE')
out.datavar.band <- Sys.getenv('OUT_DATAVAR_BAND')
out.datavar.na <- Sys.getenv('OUT_DATAVAR_NA')
out.datavar.unit <- Sys.getenv('OUT_DATAVAR_UNIT')
out.datavar.longname <- Sys.getenv('OUT_DATAVAR_LONGNAME')
out.datavar.coord.x.name <- Sys.getenv('OUT_DATAVAR_COORD_X_NAME')
out.datavar.coord.y.name <- Sys.getenv('OUT_DATAVAR_COORD_Y_NAME')
# out.datavar.coord.z.name <- Sys.getenv('OUT_DATAVAR_COORD_Z_NAME')
# out.datavar.coord.z.unit <- Sys.getenv('OUT_DATAVAR_COORD_Z_UNIT')

out.pdf.fp <- Sys.getenv('OUT_PDF_FP')
out.pdf.height <- Sys.getenv('OUT_PDF_HEIGHT')
out.pdf.width <- Sys.getenv('OUT_PDF_WIDTH')

# # start debugging
# # TODO: check paths
# cat(sprintf('in.fp==%s\n', in.fp))
# cat(sprintf('in.datavar.name==%s\n', in.datavar.name))
# cat(sprintf('in.datavar.band==%s\n', in.datavar.band))

# cat(sprintf('raster.brick==%s\n', raster.brick))
# cat(sprintf('raster.rotate==%s\n', raster.rotate))
# cat(sprintf('raster.plot==%s\n', raster.plot))

# cat(sprintf('template.fp==%s\n', template.fp))
# cat(sprintf('template.datavar.name==%s\n', template.datavar.name))

# cat(sprintf('out.fp==%s\n', out.fp))
# cat(sprintf('out.datavar.name==%s\n', out.datavar.name))
# cat(sprintf('out.datavar.type==%s\n', out.datavar.type))
# cat(sprintf('out.datavar.band==%s\n', out.datavar.band))
# cat(sprintf('out.datavar.na==%s\n', out.datavar.na))
# cat(sprintf('out.datavar.unit==%s\n', out.datavar.unit))
# cat(sprintf('out.datavar.longname==%s\n', out.datavar.longname))
# cat(sprintf('out.datavar.coord.x.name==%s\n', out.datavar.coord.x.name))
# cat(sprintf('out.datavar.coord.y.name==%s\n', out.datavar.coord.y.name))
# # cat(sprintf('out.datavar.coord.z.name==%s\n', out.datavar.coord.z.name))
# # cat(sprintf('out.datavar.coord.z.unit==%s\n', out.datavar.coord.z.unit))
# cat(sprintf('out.pdf.fp==%s\n', out.pdf.fp))
# cat(sprintf('out.pdf.height==%s\n', out.pdf.height))
# cat(sprintf('out.pdf.width==%s\n', out.pdf.width))
# #   end debugging

raster.brick <- as.logical(raster.brick)
raster.rotate <- as.logical(raster.rotate)
raster.plot <- as.logical(raster.plot)
in.datavar.band <- as.numeric(in.datavar.band)
out.datavar.na <- as.numeric(out.datavar.na)
out.datavar.band <- as.numeric(out.datavar.band)
out.pdf.height <- as.numeric(out.pdf.height)
out.pdf.width <- as.numeric(out.pdf.width)

# ensure global scope, so raster.aqmeii can be accessed inside raster.plot conditional? nope, fail
# raster.crop <- TRUE # kludge

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

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

# 1. Load the input datavar

if (raster.brick) { # multiple vertical layers
  in.raster <- brick(in.fp, varname=in.datavar.name, crs=global.proj)
  layers.n <- nbands(in.raster)
} else {
  in.raster <- raster(in.fp, varname=in.datavar.name, crs=global.proj)
  layers.n <- 1
}
# start debugging
in.raster 
# e.g.
# class       : RasterLayer 
# dimensions  : 2221, 4351, 9663571  (nrow, ncol, ncell)
# resolution  : 0.01666667, 0.01666667  (x, y)
# extent      : -131.0083, -58.49167, 22.49167, 59.50833  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 
# data source : ./ETOPO1_Ice_g_gmt4.grd_region_zS.nc 
# names       : surface.elevation 
# zvar        : zS 

netCDF.stats.to.stdout(netcdf.fp=in.fp, var.name=in.datavar.name)
# e.g.
# For ./ETOPO1_Ice_g_gmt4.grd_region_zS.nc var=zS
#       cells=9663571
#       obs=9663571
#       min=-6.63e+03
#       q1=-1.36e+03
#       med=203
#       mean=-659
#       q3=569
#       max=4.16e+03
#   end debugging

# 2. Rotate the input global datavar if requested:
#    zero-based longitudes cause problems for wrld_simpl (and probably other data)
if (raster.rotate) {
  in.raster <- rotate(in.raster,
    overwrite=TRUE) # else levelplot does one layer per page?
  in.raster # debugging
}

# 3. Create template raster.
template.raster <-projectExtent(
  raster(template.fp, varname=template.datavar.name), crs=out.crs)
# template.raster <- projectExtent(template.in.raster, crs=out.crs)
# start debugging
template.raster 
# e.g.
# class       : RasterLayer 
# dimensions  : 21, 30, 630  (nrow, ncol, ncell)
# resolution  : 2.5, 1.804511  (x, y)
# extent      : -131.25, -56.25, 22.73684, 60.63158  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +ellps=WGS84 

netCDF.stats.to.stdout(netcdf.fp=template.fp, var.name=template.datavar.name)
# e.g.
# For ./2008N2O_restart_region_PS.nc var=PS
#       cells=630
#       obs=630
#       min=7.43e+04
#       q1=9.59e+04
#       med=9.93e+04
#       mean=9.74e+04
#       q3=1.02e+05
#       max=1.03e+05
#   end debugging

# 4. Regrid input using template.
out.raster <-
  projectRaster(
    # give a template with extents--fast, but gotta calculate extents
    from=in.raster, to=template.raster, crs=out.crs,
    # give a resolution instead of a template? no, that hangs
#    from=in.raster, res=grid.res, crs=out.crs,
    method='bilinear', overwrite=TRUE, format='CDF',
    # args from writeRaster
    # datavar type will truncate unless set?
    datatype=out.datavar.type,
    NAflag=out.datavar.na,
    varname=out.datavar.name, 
    # netCDF-specific arguments
    varunit=out.datavar.unit,
    longname=out.datavar.longname,
    xname=out.datavar.coord.x.name,
    yname=out.datavar.coord.y.name,
#    zname=out.datavar.coord.z.name,
#    zunit=out.datavar.coord.z.unit
    filename=out.fp)
# if above fails to set CRS
# out.raster@crs <- CRS(out.crs)

# start debugging
# regrid.end.time <- system('date', intern=TRUE)
# cat(sprintf('  end regrid @ %s\n', regrid.end.time))
# cat(sprintf('%s\n', regrid.start.time))

out.raster 
# e.g.
# class       : RasterLayer 
# dimensions  : 21, 30, 630  (nrow, ncol, ncell)
# resolution  : 2.5, 1.804511  (x, y)
# extent      : -131.25, -56.25, 22.73684, 60.63158  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +ellps=WGS84 
# data source : ./ETOPO1_Ice_g_gmt4.grd_region_zS_regrid.nc 
# names       : surface.elevation 
# zvar        : zS 

netCDF.stats.to.stdout(netcdf.fp=out.fp, var.name=out.datavar.name)
# e.g.
# For ./ETOPO1_Ice_g_gmt4.grd_region_zS_regrid.nc var=zS
#       cells=630
#       obs=580
#       min=-6.16e+03
#       q1=-1.54e+03
#       med=200
#       mean=-658
#       q3=565
#       max=3.82e+03
# problem! significant decrease from orig max=4.16e+03, but this just regrids over same region ...

system(sprintf('ls -alht %s | head', work.dir))
system(sprintf('ncdump -h %s', out.fp))
#   end debugging

# 5. Plot data if requested

# 5.0 Massive annoyance: out.raster is out of scope inside this conditional!
# Despite numerous attempted kludges, I cannot make plotting conditional!
# if (raster.plot) {

# 5.1 Fix layer names (which are plot labels in `layerplot`)

if (raster.brick) {
  # TODO: since using rasterVis::levelplot, get/set names with `layerNames`
  aqmeii.df <- as.data.frame(out.raster)
  names(aqmeii.df) <- formatC(
    as.numeric(sub(x=names(aqmeii.df), pattern="X", replacement="", fixed=TRUE)),
    format="g", digits=sigdigs)
  # start debugging
  cat(sprintf('layers.n==%s\n', layers.n))
  ## out.raster
  ## aqmeii.df
  ## cat('head and tail of aqmeii.df:\n')
  ## head(aqmeii.df)
  ## tail(aqmeii.df)
  cat('head and tail of names(aqmeii.df):\n')
  head(names(aqmeii.df))
  tail(names(aqmeii.df))
  #  end debugging
} # end if (raster.brick)

# 5.2 Create world map (which rasterVis will autocrop)

  cat('Loading map packages ...\n') # debugging
  library(maps)
  library(maptools) # needed for map2SpatialLines
  data(wrld_simpl)      # from maptools

  # TODO: save/restore after first creation!
  global.map <- map('world', plot=FALSE)
  global.map.shp <- map2SpatialLines(global.map, proj4string=global.proj)

# 5.3 Plot and display.

  cat('Loading plot packages ...\n') # debugging
  library(rasterVis)

  cat(sprintf(
    '%s: plotting %s may take awhile! (TODO: add progress control)\n',
    this.fn, out.pdf.fp))
#  if (raster.plot) {
  pdf(file=out.pdf.fp, width=out.pdf.width, height=out.pdf.height)

if (raster.brick) {
    rasterVis::levelplot(out.raster,
      margin=FALSE,
      names.attr=names(aqmeii.df),
      layout=c(1,layers.n), # all layers in one "atmospheric" column
    ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
  } else {
    rasterVis::levelplot(out.raster,
      margin=FALSE
    ) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
    # the plot crops the map!
  } # end if (raster.brick)
  # teardown
  dev.off()

# } # end if (raster.plot)

system(sprintf('ls -alh %s', out.pdf.fp)) # debugging
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.