Summary
How to limit the display of R raster::plot boundaries of the Raster * object? Why am I asking:
I'm a newbie R who should
- convert non-projected (lon-lat) global spatial data from ASCII to netCDF (allowed)
- "regrid" it: i.e. limit the data to the area, project it and reduce the scale (reduce the size of the gridcells), it (basically solved)
Unfortunately, in scientific work, mainly QA data transformations, such as visual inspection. Step 1 above looks good. But, although step 2 now looks much better, I really want to draw only the borders (or extents ) of the registered RasterLayer . I have R code driven by bash scripts in a public git repository ( here ) that does the conversion step and displays the converted output, apparently correctly . However, my attempts to construct regression output using raster::plot not entirely correct (see the first 3 pages of output here ).
How to fix?
more details
background
I need to take some data (global marine emissions inventory (EI)), combine it with other data (mostly other EI) and enter them into the model. The model wants to consume netCDF, and it wants the netCDF above the spatial domain to be slightly larger than the adjacent US (CONUS). Borders of this image 
are the boundaries of the area. (Note that the part, but only part of the domain, is oceanic.) The model also wants the netCDF to project in a certain way ( LCC at 12 -km). I see this as two independent problems:
- convert global EI from its own ASCII format to netCDF
- "regrid" netCDF from global / non-projected to a smaller (thinner) projected subdomain.
I am trying to solve these problems with the code that I archived into this public git repository . If you clone this git repository and then configure / run the bash driver GEIA_to_netCDF.sh as described for the "first example" in README (and assuming your R is configured accordingly, especially with packages = raster , ncdf4 ), it will output netCDF and a graph (using fields::image.plot ) that are output to PDF , which should look like this: 
. The output distribution looks right, so problem 1 seems to be resolved.
Problem
The solution to problem 2 (like the "second example" in README ) seems to require R package = raster . My bash driver regrid_GEIA_netCDF.sh calls regrid.global.to.AQMEII.r (both in the repository ), which works without errors and displays a PDF file of its output. GEIA_N2O_oceanic_regrid.pdf contains 4 pages corresponding to the 4 code blocks at the end of regrid.global.to.AQMEII.r . Only the first 3 pages are important here (the fourth page tries to use fields::image.plot and has big problems).
Page = 1/4 is obtained from the simple raster::plot of the wrld_simpl output, as well as the projected version of the CONUS map from wrld_simpl :
plot(out.raster) plot(map.us.proj, add=TRUE)
Unfortunately, the output looks global or almost like this: 
But the desired domain to which the output was processed is much smaller: . So, I have 3 questions (1 main question, 2 observations):
- Does the image displayed by
raster::plot show by default only the RasterLayer extents specified in its arguments? - If so, what happened to my reggae? (below)
- If not (i.e. if
raster::plot defaults to display more than the extents of its RasterLayer ), how can raster::plot display only extents of this RasterLayer ?
The code is in regrid.global.to.AQMEII.r ( here ), which apparently returns the correct result:
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 seems to match the output domain definition specified here here . (Once on the linked page, scroll 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 already noted, the regression output ( out.raster ) shows that it is lon-lat in its CRS: I'm not sure why this is, or if it implies that out.raster is global in extents.
I also tried to limit the plot itself in two ways:
Firstly, I tried to add extents to the plot, i.e.
plot(out.raster, ext=template.raster) plot(map.us.proj, add=TRUE)
which creates the page <2/4 PDF . Unfortunately, this does not change the output at all: AFAICS, it is completely identical to the page = 1/4 (see above).
Secondly, I tried using raster::crop to bind the most corrected netCDF before building it using the same RasterLayer object ( template.raster ) that I used to bind regression:
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 a page = 3/4 pdf . Unfortunately, this also seems to be completely identical to page = 1/4 (see above).
(If you're interested, page 4 of the PDF was generated using fields::image.plot , which has another problem described here if StackOverflow doesn't hit the link.)
Your help to this newbie R is appreciated!