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

[Michael Sumner]

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 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)


## 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 <mdsumner at gmail.com> wrote:
> 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



-- 
Michael Sumner
Hobart, Australia
e-mail: mdsumner at gmail.com