I happened to write this yesterday, which switches the parts of wrld_simpl from < 0 longitude to the east (Pacific view). I haven't investigated what you need in detail, so just ignore if I'm off track.
It's not careful about data or anything, since all I needed was the raw geometry. It might be of use. "xrange" here can be anything from [-180, 360] though it's really assuming you don't traverse more than 360 in total. It would be nice to carefully clean up a copy of wrld_simpl like this and make a wrld_simpl version like the maps package has. Cheers, Mike. xrange <- c(140, 250) ##xrange <- c(0, 360) yrange <- c(-60, 20) require(maptools) require(raster) require(rgeos) ext <- extent(xrange[1], xrange[2], yrange[1], yrange[2]) data(wrld_simpl, package = "maptools") if (xrange[2] > 180.0) { eastrange <- c(xrange[1], 180.0) westrange <- c(-180.0, xrange[2] - 360) ew <- extent(westrange[1], westrange[2], yrange[1], yrange[2]) ee <- extent(eastrange[1], eastrange[2], yrange[1], yrange[2]) geome <- as(ee, "SpatialPolygons") geomw <- as(ew, "SpatialPolygons") ## why does this get dropped? proj4string(geome) <- CRS(proj4string(wrld_simpl)) proj4string(geomw) <- CRS(proj4string(wrld_simpl)) worlde <- gIntersection(wrld_simpl, geome) worldw <- gIntersection(wrld_simpl, geomw) worldw <- elide(worldw, shift = c(360.0, 0.0)) proj4string(worldw) <- CRS(proj4string(wrld_simpl)) dfe <- data.frame(x = seq_len(length(worlde@polygons))) dfw <- data.frame(x = seq_len(length(worldw@polygons))) rownames(dfw) <- as.character(as.numeric(rownames(dfe)) + nrow(dfe)) worlde <- SpatialPolygonsDataFrame(worlde, dfe, match.ID = FALSE) worldw <- SpatialPolygonsDataFrame(worldw, dfw, match.ID = FALSE) worldw <- spChFIDs(worldw, rownames(dfw)) world <- spRbind(worlde, worldw) } else { geom <- as(ext, "SpatialPolygons") proj4string(geom) <- CRS(proj4string(wrld_simpl)) ## TODO: make this spdf if needed world <- gIntersection(wrld_simpl, geom) } plot(world) On Tue, Dec 18, 2012 at 8:15 PM, Tom Roche <tom_ro...@pobox.com> wrote: > > 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_ro...@pobox.com> > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Michael Sumner Hobart, Australia e-mail: mdsum...@gmail.com _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo