problems drawing a world map on a levelplot
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 at polygons)))
dfw <- data.frame(x = seq_len(length(worldw at 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_Roche at 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_Roche at pobox.com>
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Michael Sumner Hobart, Australia e-mail: mdsumner at gmail.com