[R-sig-Geo] problems drawing a world map on a levelplot
Michael Sumner
mdsumner at gmail.com
Wed Dec 19 05:55:10 CET 2012
I meant "wrld_simpl2 version like the maps package has".
On Tue, Dec 18, 2012 at 8:54 PM, Michael Sumner <mdsumner at 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 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
--
Michael Sumner
Hobart, Australia
e-mail: mdsumner at gmail.com
More information about the R-sig-Geo
mailing list