[R] [newbie] failure to plot a RasterLayer with raster::plot or fields::image.plot

[Tom Roche]

summary: spatial data to be input to a regional-scale environmental
model must (1) be converted to netCDF and then (2) "regridded" (cropped,
projected, increased resolution). In a public git repository

https://github.com/TomRoche/GEIA_to_NetCDF

I have R code (with bash drivers) that does the conversion step, and
plots the converted output, apparently correctly. However attempts to
plot the output of the regridding step (twice with raster::plot, once
with fields::image.plot) are very wrong. How to fix?

details:

I need to take some data (a global marine emissions inventory (EI)),
combine it with other data (mostly other EIs), and input that to a
model. The model wants to consume netCDF, and it wants that netCDF
over a certain domain (a bit bigger than the contiguous US (CONUS)),
and it wants that netCDF projected a certain way (LCC at 12-km
resolution). I'm doing this in 2 steps

1. convert the global EI from its native ASCII format to netCDF

2. "regrid" the netCDF from global/unprojected to a finer-resolution,
   projected subdomain.

and doing this @

https://github.com/TomRoche/GEIA_to_NetCDF

If you clone that git repo, and then configure/run GEIA_to_netCDF.sh
as described for the "first example" @

https://github.com/TomRoche/GEIA_to_netCDF/blob/master/README.md

(and presuming your R is appropriately configured, notably with
packages=fields, ncdf4) it will output netCDF and plot the
output to a PDF, which is also available @

https://github.com/downloads/TomRoche/GEIA_to_netCDF/GEIA_N2O_oceanic.pdf

). The distribution of the output data appears very similar to the input
data, and the output plot appears correct, so problem 1 seems solved.

Problem 2 (aka the "second example") is somewhat more difficult. Driver
regrid_GEIA_netCDF.sh calls regrid.global.to.AQMEII.r, which runs
without error. It also *seems* to regrid properly, in the sense that all
the API rules seem to have been followed. But something is very wrong
with my attempt to plot, which is currently available @

https://github.com/downloads/TomRoche/GEIA_to_netCDF/GEIA_N2O_oceanic_regrid.pdf

GEIA_N2O_oceanic_regrid.pdf comprises 3 pages, corresponding to 3 blocks of code at the end of regrid.global.to.AQMEII.r:

* page 1: results from a simple raster::plot of the regridded raster,
  plus a projected version of a CONUS map from wrld_simpl:

map.us.unproj <- wrld_simpl[wrld_simpl$ISO3 == 'USA', ]  # unprojected
map.us.proj <- spTransform(map.us.unproj, CRS(out.crs))  # projected

  Unfortunately the output is ... odd. The EI data is represented by a
  thin vertical smear down the middle of the plot space. The US map is
  oddly distorted, almost like seeing the Atlantic coast edge-on. The
  good news is, if one enlarges the page enough (e.g., to 800%), one
  sees

+ the space inside the US-map borders is white, which is good, since
  the data is marine

+ there is a roughly Canadian-shaped white blob oriented correctly
  relative to the CONUS blob, and similarly a roughly Mexico-shaped
  blob

+ there are green blotches corresponding to the elevated emissions off
  the Pacific coasts of Canada and Mexico, and they are situated
  similar to their position in the unprojected world map

https://github.com/downloads/TomRoche/GEIA_to_netCDF/GEIA_N2O_oceanic.pdf

  So I'm thinking that the data might have regridded correctly
  (excepting the extents problem--more below), but that I'm plotting it
  very wrongly :-( How to fix?

* page 2: results from a simple raster::plot of the regridded raster,
  with plot extents set (plus the CONUS map):

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

  Unfortunately AFAICS that made no difference whatsoever; page 2 is
  idential to page 1 :-(

* page 3: I have previously successfully used fields::image.plot
  (e.g., in the first problem, above), so I reverted to that. But the
  plot is, in some ways, even worse than the previous pages. The
  western hemisphere, esp CONUS, is much more recognizable. But the
  data is much worse: it's diagonal streaks, like old-fashioned TV
  noise. Something is very wrong :-(

The worst part is that, in none of the above 3 pages does the output
appear to be appropriately bounded. The other outputs to the model
are restricted to the AQMEII-NA domain (see map of extents @

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

) so I need to get this data similarly restricted. That the data is not
restricted to the appropriate spatial subdomain also seems to outweigh
the positive indicators above regarding the regridding (e.g., map
position, position of data overall, position of emission peaks).

Hence I'd appreciate very much if someone with more experience with
spatial plotting would take a look at the output and the code, and
diagnose the problem. Note that both output and code will remain public,
hopefully to be useful to the next person attempting this sort of task.
And of course I will be sure to annotate the code and project docs to
give credit to whoever helps fix this!

thanks in advance, Tom Roche <Tom_Roche at pobox.com>