[R-sig-Geo] how to plot Raster*? problems with raster::plot, fields::image.plot

[Tom Roche]

Evidence for a bug in raster::plot :

I tried adding the code you suggested. First, the additional args to
raster(...)

https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016501.html
>> template.raster <- raster(nrow=299, ncol=459, xmn=-14802449, xmx=9694173, ymn=-6258782, ymx=10749443, crs="+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556000 +y_0=-1728000 +ellps=WGS84")

but since I am creating my template.raster using projectExtent(...), for which those args are unused, I rolled back.

However the different call to spTransform(...) was quite suggestive ...

>> map.us.proj <- spTransform(map.us.unproj, CRS=projection(out.raster, asText=F))

... because it failed: in my plot output, the map did not draw.

This made me think about *why* that might have failed, which led me to a
question previously noted in a comment in my code (starred below):

https://github.com/TomRoche/GEIA_to_netCDF/blob/master/regrid.global.to.AQMEII.r
>>> 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
>>> #> 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

That would explain why your syntax

>> map.us.proj <- spTransform(map.us.unproj, CRS=projection(out.raster, asText=F))

failed: projectRaster(...), at least in this case, is failing to
properly record its projection in out.raster. (Perhaps some slot is not
properly assigned?) This also seems to suggest why raster::plot, when
called using out.raster

https://github.com/TomRoche/GEIA_to_netCDF/blob/master/regrid.global.to.AQMEII.r
>>> plot(out.raster)
or
>>> plot(out.raster, ext=template.raster)
or
>>> template.raster.extent <- extent(template.raster)
>>> out.raster.crop <-
>>>   crop(out.raster, template.raster.extent, filename=out.fp, overwrite=TRUE)
>>> plot(out.raster.crop)

is producing a plot with *global* extents

https://github.com/downloads/TomRoche/GEIA_to_netCDF/raster..plot_bounds_problem.png

instead of respecting the extents of its projection ... because
projectRaster(...) is failing to provide out.raster with the correct
projection!

Am I missing something? Because, frankly, your previous suggestion

>> I take it you want a smaller extent. Then use a template with a
>> better extent

makes no sense, as previously noted:

https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016502.html
> why do I need to choose *one* set of extents to get the regridded
> *data* that I need, and *another* set of extents to display that
> data with raster::plot?

FWIW, Tom Roche <Tom_Roche at pobox.com>