How to set the display bounds of a raster::plot?

Question

summary

How to limit display of an R raster::plot to the bounds of a Raster* object? Why I ask:

I'm an R beginner who must

1. convert unprojected (lon-lat) global spatial data from ASCII to netCDF (solved)
2. "regrid" it: i.e., restricted the data to a region, project it, and downscale (decrease the size of the gridcells) it (mostly solved)

Unfortunately, in scientific work, one mostly QAs data transformations such as this by visual inspection. Step 1 above looks good. But though step 2 now looks much better, I really want to plot only the bounds (or extents) of the regridded RasterLayer. I have R code driven by bash scripts in a public git repository (here) that does the conversion step, and plots the converted output, apparently correctly. However my attempts to plot the output of the regridding with raster::plot are not quite right (see first 3 pages of output here).

How to fix?

details

background

I need to take some data (a global marine emissions inventory (EI)), combine it with other data (mostly other EIs), and input that to a model. The model wants to consume netCDF, and it wants that netCDF over a spatial domain slightly bigger than the contiguous US (CONUS). The borders of this image are the boundaries of the domain. (Note that part, but only part, of the domain is oceanic.) The model also wants that netCDF projected a certain way (LCC at 12-km resolution). I'm treating this as 2 independent problems:

1. convert the global EI from its native ASCII format to netCDF
2. "regrid" the netCDF from global/unprojected to the downscaled (finer-resolution), projected subdomain.

I'm attempting to solve these problems with code I have archived at this public git repository. If you clone that git repo, and then configure/run the bash driver GEIA_to_netCDF.sh as described for the "first example" in the README (and presuming your R is appropriately configured, notably with packages=raster, ncdf4) it will output netCDF and plot (using fields::image.plot) that output to a PDF, which should look like this: . The distribution of the output data appears correct, so problem 1 seems solved.

problem

Solving problem 2 (aka the "second example" in the README) seems to require R package=raster. My bash driver regrid_GEIA_netCDF.sh calls regrid.global.to.AQMEII.r (both in repository), which runs without error, and displays a PDF of its output. GEIA_N2O_oceanic_regrid.pdf comprises 4 pages, corresponding to 4 blocks of code at the end of regrid.global.to.AQMEII.r. Only the first 3 pages are relevant here (the 4th page tries to use fields::image.plot and has bigger problems).

Page=1/4 results from a simple raster::plot of the regridded output, plus a projected version of a CONUS map from wrld_simpl:

plot(out.raster)
plot(map.us.proj, add=TRUE)

Unfortunately the output appears global, or nearly so: But the desired domain, to which the output was regridded, is much smaller: . So I have 3 questions (1 main question, 2 followups):

1. Does the image displayed by raster::plot by default display only the extents of the RasterLayer given in its arguments?
2. If so, what is wrong with my regrid? (more below)
3. If not (i.e., if raster::plot by default displays more than the extents of its RasterLayer), how to make raster::plot display only the extents of that RasterLayer?

The code in regrid.global.to.AQMEII.r (here) that does the regridding seems to get the correct output:

out.crs <- '+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556000 +y_0=-1728000'

Note that CRS appears to match the definition of the output domain given here. (Once at the linked page, scroll past the map.)

out.raster <-
  projectRaster(
from=in.raster, to=template.raster, method='bilinear', crs=out.crs,
overwrite=TRUE, progress='window', format='CDF',
# args from writeRaster
NAflag=-999.0,  # match emi_n2o:missing_value,_FillValue (TODO: copy)
varname=data.var.name, 
varunit='ton N2O-N/yr',
longname='N2O emissions',
xname='COL',
yname='ROW',
filename=out.fp)
out.raster
# made with CRS
#> class       : RasterLayer 
#> dimensions  : 299, 459, 137241  (nrow, ncol, ncell)
#> resolution  : 53369.55, 56883.69  (x, y)
#> extent      : -14802449, 9694173, -6258782, 10749443  (xmin, xmax, ymin, ymax)
# ??? why still proj=longlat ???
#> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#> data source : /home/rtd/code/R/GEIA_to_netCDF/GEIA_N2O_oceanic_regrid.nc
#> names       : N2O.emissions 
#> zvar        : emi_n2o 

But, as noted, the regrid output (out.raster) shows itself to be lon-lat in its CRS: I'm not sure why that is, or if it implies out.raster is global in extents.

I also tried to restrict the plot itself in two ways:

Firstly, I tried adding an extents to the plot, i.e.,

plot(out.raster, ext=template.raster)
plot(map.us.proj, add=TRUE)

which generates page=2/4 of the PDF. Unfortunately, that does not change the output at all: AFAICS, it is completely identical to page=1/4 (above).

Secondly, I tried using raster::crop to bound the regridded netCDF itself, before plotting that, using the same RasterLayer object (template.raster) that I used to bound the regrid:

template.raster.extent <- extent(template.raster)
out.raster.crop <-
  # overwrite the uncropped file
  crop(out.raster, template.raster.extent, filename=out.fp, overwrite=TRUE)
...
plot(out.raster.crop)
plot(map.us.proj, add=TRUE)

which generates page=3/4 of the PDF. Unfortunately, that it is also apparently completely identical to page=1/4 (above).

(In case you're wondering, page 4 of the PDF was generated using fields::image.plot, which has a different problem, described here, unless StackOverflow whacks that link.)

Your assistance to this R newbie is appreciated!

Answer 1

summary

My question

How to limit display of an R `raster::plot` to the bounds of a `Raster*` object?

was based on my incorrect diagnosis of the problem. The solution was to set the bounds (or extent) of the underlying Raster* data object (i.e., the source of the data for the plot) correctly, and particularly

1. to get the bounds of the template object (used to create the data object) by querying its underlying file (make no assumptions!)
2. to get the CRS of the template object by querying its underlying file (make no assumptions!)
3. to set the bounds and CRS for the template object directly

But

• don't try to just set the resolution on the Extent; for me at least, that hangs
• there is still one minor bounding problem with the raster::plot map; at least, relative to the map I'm using with fields::image.plot

details

The answer to this question is actually the answer to a different question: how to correctly set the extent of a Raster* object? because the raster::plot problem mostly went away once I properly set the extent of the object I wanted to plot. (With one minor exception--see below.) This may have been what Robert J. Hijmans, the main raster developer, was trying to convey in this post, but if so, I did not perceive it at the time.

I solved the problem after getting two suggestions which weren't themselves what I needed, but which set me on the proper path. In this case, that path led to the R package M3, which is useful for dealing with data input to and output by CMAQ and WRF. The suggestions were

1. Use project.lonlat.to.M3 for the regridding task.
2. The plot problem may be related to the assumption by the generating model (CMAQ) of a spherical projection.

The second suggestion seemed plausible, since I knew that I was getting and setting a CRS with +ellps=WGS84. (See profusely-commented R in this repository, specifically regrid.global.to.AQMEII.r.) So I resolved to look into that.

I knew the first suggestion would only work with difficulty (since it only projects points), but looked at the M3 doc just to be sure. The following functions in its ToC immediately caught my eye:

1. get.proj.info.M3, which not only returned a spherical PROJ.4 string from the template file, but one without the false eastings and northings which I believed I needed to supply:

   # 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
   

2. get.grid.info.M3, which can, with a bit of effort, return the bounds/extent of the template file:

   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)
   template.extents <-
     extent(extents.xmin, extents.xmax, extents.ymin, extents.ymax)
   

One can then set those bounds on the template Raster*

template.in.raster <- raster(template.in.fp, ...)
template.raster <- projectExtent(template.in.raster, crs=out.crs)
template.raster@extent <- template.extents

and use the template to regrid the input Raster*

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=-999.0,  # match emi_n2o:missing_value,_FillValue (TODO: copy)
    varname=data.var.name, 
    varunit='ton N2O-N/yr',
    longname='N2O emissions',
    xname='COL',
    yname='ROW',
    filename=out.fp)
# above fails to set CRS, so
out.raster@crs <- CRS(out.crs)

(As suggested above, regridding by setting the template took only 7 sec, but regridding by setting a grid resolution failed to complete after 2 hr, when I killed the job.) After that, the raster::plot

map.us.unproj <- wrld_simpl[wrld_simpl$ISO3 %in% c('CAN', 'MEX', 'USA'),]
map.us.proj <-
  spTransform(map.us.unproj, CRS(out.crs)) # projected
...
pdf(file=pdf.fp, width=5.5, height=4.25)
...
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)

was nearly as expected, with one exception: the map extends beyond the bounds of the data to the north and south (though not to the east and west), causing the image to insufficiently resemble published images of the domain, e.g. this. Interestingly, when I plot with fields::image.plot like

# see code in
# https://github.com/TomRoche/GEIA_to_netCDF/blob/master/plotLayersForTimestep.r
plot.raster(
  raster=out.raster,
  title=title,
  subtitle=subtitle,
  q.vec=probabilities.vec,
  colors,
  map.cmaq
)

I don't get that problem: . So I'll probably use fields::image.plot for display, at least for the moment.