StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
text
Body
copied!<h1>summary</h1> <p>My question </p> <pre><code>How to limit display of an R `raster::plot` to the bounds of a `Raster*` object? </code></pre> <p>was based on my incorrect diagnosis of the problem. The solution was to set the bounds (or <code>extent</code>) of the underlying <code>Raster*</code> <em>data</em> object (i.e., the source of the data for the plot) correctly, and particularly</p> <ol> <li>to get the bounds of the template object (used to create the data object) by querying its underlying file (make no assumptions!)</li> <li>to get the CRS of the template object by querying its underlying file (make no assumptions!)</li> <li>to set the bounds and CRS for the template object directly</li> </ol> <p>But</p> <ul> <li>don't try to just set the resolution on the <code>Extent</code>; for me at least, that hangs</li> <li>there is still one minor bounding problem with the <code>raster::plot</code> map; at least, relative to the map I'm using with <code>fields::image.plot</code></li> </ul> <h1>details</h1> <p>The answer to this question is actually the answer to a different question: how to correctly set the <code>extent</code> of a <code>Raster*</code> object? because the <code>raster::plot</code> 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 <code>raster</code> developer, was trying to convey in <a href="https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016501.html" rel="nofollow noreferrer">this post</a>, but if so, I did not perceive it at the time.</p> <p>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 <a href="http://cran.r-project.org/web/packages/M3/M3.pdf" rel="nofollow noreferrer">R package <code>M3</code></a>, which is useful for dealing with data input to and output by CMAQ and WRF. The suggestions were</p> <ol> <li>Use <code>project.lonlat.to.M3</code> for the regridding task.</li> <li>The plot problem may be related to the assumption by the generating model (CMAQ) of a spherical projection.</li> </ol> <p>The second suggestion seemed plausible, since I knew that I was getting and setting a CRS with <code>+ellps=WGS84</code>. (See profusely-commented R in <a href="http://github.com/TomRoche/GEIA_to_netCDF" rel="nofollow noreferrer">this repository</a>, specifically <a href="http://github.com/TomRoche/GEIA_to_netCDF/blob/master/regrid.global.to.AQMEII.r" rel="nofollow noreferrer"><code>regrid.global.to.AQMEII.r</code></a>.) So I resolved to look into that.</p> <p>I knew the first suggestion would only work with difficulty (since it only projects points), but looked at the <a href="http://github.com/TomRoche/GEIA_to_netCDF" rel="nofollow noreferrer">M3 doc</a> just to be sure. The following functions in its ToC immediately caught my eye:</p> <ol> <li><p><code>get.proj.info.M3</code>, 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:</p> <pre><code># 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 </code></pre></li> <li><p><code>get.grid.info.M3</code>, which can, with a bit of effort, return the bounds/extent of the template file:</p> <pre><code>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) </code></pre></li> </ol> <p>One can then set those bounds on the template <code>Raster*</code></p> <pre><code>template.in.raster <- raster(template.in.fp, ...) template.raster <- projectExtent(template.in.raster, crs=out.crs) template.raster@extent <- template.extents </code></pre> <p>and use the template to regrid the input <code>Raster*</code></p> <pre><code>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) </code></pre> <p>(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 <code>raster::plot</code> </p> <pre><code>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) </code></pre> <p>was nearly as expected, with one exception: <img src="https://i.stack.imgur.com/lEzSZ.png" alt="excess north-south plot extents with raster::plot"> 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. <a href="http://github.com/TomRoche/cornbeltN2O/wiki/images/Cmaq_aqmeii_domain.png" rel="nofollow noreferrer">this</a>. Interestingly, when I plot with <code>fields::image.plot</code> like</p> <pre><code># 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 ) </code></pre> <p>I don't get that problem: <img src="https://i.stack.imgur.com/OPXO6.png" alt="fields::image.plot after problem fixed">. So I'll probably use <code>fields::image.plot</code> for display, at least for the moment.</p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload