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 - convert unprojected (lon-lat) global spatial data from ASCII to netCDF (solved) - "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…

Accepted Answer

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

- to get the bounds of the template object (used to create the data object) by querying its underlying file (make no assumptions!)

- to get the CRS of the template object by querying its underlying file (make no assumptions!)

- 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

- Use `project.lonlat.to.M3` for the regridding task.

- 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:

`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
', out.crs)) # debugging
# out.crs=+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000
```

`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…

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

Solution

Problem