[R-sig-Geo] problems drawing a world map on a levelplot

Tom Roche Tom_Roche at pobox.com
Wed Dec 19 05:15:48 CET 2012


summary: I've tried several ways to draw wrld_simpl over global data
plotted with rasterVis::levelplot, but keep encountering the same
problem:

- only the eastern hemisphere of the map plots

- the eastern hemisphere of the map overlays the western hemisphere of
  the data

How to fix? or otherwise provide geographical context?

details:

(I'll attempt to produce a "complete, self-contained example" after I
make the code publicly viewable, but for now, I'm hoping the problem is
sufficiently simple to diagnose.)

I have atmospheric data in a netCDF file specifying gas concentrations
over a 3D space with dimensions longitude, latitude, and (vertical)
level, such that

* the horizontal extents are global 

* the vertical coordinates are pressures (specifically hybrid
  sigma-pressure)

package=raster provides very usable API for loading the data into a
RasterBrick, and package=rasterVis makes it relatively easy to levelplot
the data in a single "atmospheric column" with the highest pressures at
the bottom. However I also want to provide some geographic context
(e.g., to overlay the data with national boundaries), and have failed to
date.

I have done this with maptools::wrld_simpl successfully in other
contexts, i.e.,

* loading the data with package=ncdf4 (which is considerably more work)
* converting the data to a data.frame
* lattice::levelplot-ing the data
* lattice::xyplot-ing the map

But when I do

data(wrld_simpl)
global.raster <- brick(in.fp, varname=data.var.name)
layers.n <- nbands(global.raster)
global.data.df <- as.data.frame(global.raster)
# some twiddling of names(global.data.df) omitted
global.map <- map('world', plot=FALSE)
global.map.df <- data.frame(lon=global.map$x, lat=global.map$y)
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
  margin=FALSE,
  names.attr=names(global.data.df),
  layout=c(1,layers.n), # all layers in one "atmospheric column"
) + xyplot(lat ~ lon, global.map.df, type='l', lty=1, lwd=1, col='black')
dev.off()

the data plots apparently correctly, but the map plots incorrectly, as
shown by this .png of the top level:

https://docs.google.com/open?id=0BzDAFHgIxRzKOVdlSDFMR29BdWM

As noted above, the map's eastern hemisphere overlays the data's western
hemisphere, and nothing overlays the data's eastern hemisphere.

I tried several variations on the above, but either got the same
results, or outright error. (Unfortunately the error messages are often
printed on the levelplot, but the messages are truncated (so they can't
be understood), and it's annoying to "fail slow" (after waiting for a
long plot to finish).) So after some research I sought to emulate the
means recommended @

http://rastervis.r-forge.r-project.org/#levelplot

data(wrld_simpl)
global.raster <- brick(in.fp, varname=data.var.name)
layers.n <- nbands(global.raster)
global.data.df <- as.data.frame(global.raster)
# some twiddling of names(global.data.df) omitted
global.map.proj <- CRS('+proj=latlon +ellps=WGS84')
global.map <- map('world', plot=FALSE)
global.map.shp <- map2SpatialLines(global.map, proj4string=global.map.proj)
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
  margin=FALSE,
  names.attr=names(global.data.df),
  layout=c(1,layers.n), # all layers in one "atmospheric column"
) + layer(sp.lines(global.map.shp, lwd=0.8, col='darkgray'))
dev.off()

This also fails, and very similarly (though the line colors differ). See
this .png:

https://docs.google.com/open?id=0BzDAFHgIxRzKYzd6cVB2TDE4RE0

I'm guessing I can either do something simple to wrld_simpl's longitude
values, or Something Completely Different that would Just Work Better,
but don't know. Your assistance is appreciated.

TIA, Tom Roche <Tom_Roche at pobox.com>



More information about the R-sig-Geo mailing list