Hello,
Perhaps the problem could be related with the CRS setting you are using.
Instead of:
global.map.proj <- CRS('+proj=latlon +ellps=WGS84')
you should try:
global.map.proj <- CRS('+proj=longlat +ellps=WGS84')
Besides, following your code I think you are not setting the CRS of your
RasterBrick.
I have tried successfully the next toy example:
global.map.proj <- CRS('+proj=longlat +ellps=WGS84')
global.map <- map('world', plot=FALSE)
global.map.shp <- map2SpatialLines(global.map, proj4string=global.map.proj)
r <- raster()
r <- init(r, fun=function(x)1)
projection(r) <- global.map.proj
levelplot(r) + layer(sp.lines(global.map.shp))
Besides, you can modify the theme to fill the sea (NA values) with color:
myTheme <- modifyList(RdBuTheme(),
list(panel.background=list(col='skyblue3')))
## you won't see any differences here because r does not contains NA cells.
levelplot(r, par.settings=myTheme) + layer(sp.lines(global.map.shp))
Finally, you may find useful these sources for the countries boundaries:
http://www.diva-gis.org/Data
http://www.naturalearthdata.com/downloads/110m-cultural-vectors/
Best,
Oscar.
2012/12/19 Michael Sumner <mdsum...@gmail.com>
> Here's a crack at a full example, sorry the clip/merge polygon stuff
> is a mess still.
>
> library(maptools)
> library(raster)
>
> ## first, do the jiggery pokery of wrld_simpl to
> ## convert from Atlantic to Pacific view
> xrange <- c(0, 360)
> yrange <- c(-90, 90)
>
> require(maptools)
> require(raster)
> require(rgeos)
>
> ext <- extent(xrange[1], xrange[2], yrange[1], yrange[2])
> data(wrld_simpl, package = "maptools")
>
> ## this is a bit more general than needed for this example, since I
> ## wanted arbitrary crops originally
>
> 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)
>
>
> ## so, "world" is longitudes [0, 360], which is pretty much essential
> ## given the input data shown
>
> r <- raster(nrows = 180, ncols = 360, xmn = 0, xmx = 360, ymn = -90,
> ymx = 90, crs = "+proj=longlat +ellps=WGS84 +over")
>
> ## fake a raster so we have something equivalent to your data (but
> simple/r)
> rast <- rasterize(world, r)
>
> ## fill the ocean parts with something . . .
> rast[is.na(rast)] <- 1:sum(is.na(getValues(rast)))
>
> library(rasterVis)
> levelplot(rast, margin = FALSE) + layer(sp.polygons(world, lwd=0.8,
> col='darkgray'))
>
>
>
>
>
> On Tue, Dec 18, 2012 at 8:55 PM, Michael Sumner <mdsum...@gmail.com>
> wrote:
> > I meant "wrld_simpl2 version like the maps package has".
> >
> > On Tue, Dec 18, 2012 at 8:54 PM, Michael Sumner <mdsum...@gmail.com>
> wrote:
> >> 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
> >
> >
> >
> > --
> > Michael Sumner
> > Hobart, Australia
> > e-mail: mdsum...@gmail.com
>
>
>
> --
> 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
>
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
