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_ro...@pobox.com> ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
