Skip to content

problems drawing a world map on a levelplot

11 messages · Michael Sumner, Tom Roche, Matthew Landis +1 more

#
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>
#
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:

  
    
#
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:

  
    
#
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:

  
    
#
https://stat.ethz.ch/pipermail/r-sig-geo/2012-December/017014.html
(rearranged)
True, but since it was also global data, brick(...) "did the right
thing." I have now set its CRS just to be sure.
Thanks again, but please note that I was following (see starred line)

http://rastervis.r-forge.r-project.org/#levelplot
** proj <- CRS('+proj=latlon +ellps=WGS84')
...
Perhaps you could fix that?

TIA, Tom Roche <Tom_Roche at pobox.com>
#
summary: How to best/easiest range-shift lon-lat data?

details:

https://stat.ethz.ch/pipermail/r-sig-geo/2012-December/017011.html
<self-administered dopeslap/> Of course! the problem is my input!
see starred line:

R session
* extent      : -1.25, 358.75, -90.94737, 90.94737  (xmin, xmax, ymin, ymax)
which is wrong:

- latitudes outside the range [-90, +90]?

- (as you note) 0-based longitudes are likely troublesome (I'm pretty
  sure the other data with which this must integrate is [-180, +180])

- ... and these longitudes aren't even 0-based!

I didn't even check :-( since the data's source is much more eminent
than a mere student like myself. But now I note

$ ncdump -h ./2008N2O_restart.nc | head
$ ncdump -v lat ./2008N2O_restart.nc | tail -n 28
I'm not sure why raster is reporting y range=[-90.94737, 90.94737],
except that 2 * 0.94737 = 1.89474 == the latitude resolution.

$ ncdump -v lon ./2008N2O_restart.nc | tail -n 15
Again, it seems raster is reporting range(x)=[-1.25, 358.75] only
because resolution(x)=2.5?. I'll defer to its superior wisdom.

Given my need to integrate this data with other more conventional,
longitudes=[-180, +180] datasets, it seems better to fix this data than
to fix wrld_simpl. So I'm wondering, how best to "shift" my data from
longitudes=[0, 360] to longitudes=[-180, +180]?

TIA, Tom Roche <Tom_Roche at pobox.com>
#
Tom Roche Wed, 19 Dec 2012 13:52:38 -0500
Matt Landis Wed, 19 Dec 2012 14:09:57 -0500
Thanks! Now this code-------------------------------------------------

global.proj <- CRS('+proj=longlat +ellps=WGS84')
global.raster <- brick(in.fp, varname=data.var.name, crs=global.proj)
global.raster
# class       : RasterBrick 
# dimensions  : 96, 144, 13824, 56  (nrow, ncol, ncell, nlayers)
# resolution  : 2.5, 1.894737  (x, y)
# extent      : -1.25, 358.75, -90.94737, 90.94737  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
layers.n <- nbands(global.raster)

# make longitudes NOT zero-based
global.raster <- 
  rotate(global.raster,
    overwrite=TRUE) # else levelplot does one layer per page?
global.raster
# class       : RasterBrick 
# dimensions  : 96, 144, 13824, 56  (nrow, ncol, ncell, nlayers)
# resolution  : 2.5, 1.894737  (x, y)
# extent      : -181.25, 178.75, -90.94737, 90.94737  (xmin, xmax, ymin, ymax)
# ...

# TODO: if using rasterVis::levelplot, get names with `layerNames`
global.df <- as.data.frame(global.raster)
names(global.df) <- formatC(as.numeric(sub(x=names(global.df), pattern="X", replacement="", fixed=TRUE)), format="g", digits=sigdigs)
# debugging
summary(unlist(global.df))
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
# 0.01459 0.22440 0.48350 0.36550 0.48530 0.50210
# matches known stats

global.map <- map('world', plot=FALSE)
global.map.shp <- map2SpatialLines(global.map, proj4string=global.proj)

# start plot driver
cat(sprintf(
  '%s: plotting %s (may take awhile! TODO: add progress control)\n',
  this.fn, pdf.fp))
pdf(file=pdf.fp, width=5, height=100)
# using rasterVis::levelplot, not lattice::levelplot
levelplot(global.raster,
  margin=FALSE,
  names.attr=names(global.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()

--------------------------------------------------------------produces

https://docs.google.com/open?id=0BzDAFHgIxRzKaGhMOE5BblpEeVE

with both map hemispheres overlaying the data.

thanks all! Tom Roche <Tom_Roche at pobox.com>