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

[Tom Roche]

https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016501.html
> library(raster)
> # replacing your file with example data
> in.raster <- raster()
> in.raster[] <- 1:ncell(in.raster)
> 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")

I didn't know about the additional arguments (e.g., 'xmn'): they're not in the doc

http://cran.r-project.org/web/packages/raster/raster.pdf

> out.raster <- projectRaster(from=in.raster, to=template.raster, method='bilinear')

Actually

https://stat.ethz.ch/pipermail/r-sig-geo/2012-October/016500.html
>> 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
>>     varname=data.var.name, 
>>     varunit='ton N2O-N/yr',
>>     longname='N2O emissions',
>>     xname='COL',
>>     yname='ROW',
>>     filename=out.fp)

> library(maptools)
> data(wrld_simpl)
> map.us.unproj <- wrld_simpl[wrld_simpl$ISO3 %in% c('CAN', 'MEX', 'USA'),]
> map.us.proj <- spTransform(map.us.unproj, CRS=projection(out.raster, asText=F))

didn't know about that CRS parameterization, thanks

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

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

What? I want to plot the *same* extent as my data, which is the *same*
extent as was used to create the output, i.e.,

>> template.in.raster <- raster(template.in.fp, varname=data.var.name, band=template.band)
>> template.raster <- projectExtent(template.in.raster, crs=out.crs)
>> template.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)
>> # coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556000 +y_0=-1728000 +ellps=WGS84 

which not only I, but many others, know, from years of use of that
spatial domain by various air-quality models for an international
model-intercomparison project, AQMEII (about which see, e.g., the March
2012 _Atmospheric Environment_ or

http://aqmeii.jrc.ec.europa.eu/

) is a 12-km LCC grid with 299 rows and 459 columns (not

> temp2 <- raster(nrow=70, ncol=117,

(from where comes those row and col counts?)

> xmn=-5622887, xmx=621350, ymn=-3642132, ymx=339726.7, crs="+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=-2556000 +y_0=-1728000 +ellps=WGS84")

) and which plots like

https://github.com/TomRoche/cornbeltN2O/wiki/images/Cmaq_aqmeii_domain.png

Don't believe me? see, e.g., figure 1 of

http://dx.doi.org/10.1016/j.atmosenv.2011.12.041

So please explain: 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? Or am I missing something?

TIA, Tom Roche <Tom_Roche at pobox.com>