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

[Tom Roche]

Tom Roche Wed, 19 Dec 2012 13:52:38 -0500
>> How to best/easiest range-shift lon-lat data?

Matt Landis Wed, 19 Dec 2012 14:09:57 -0500
> [raster::shift] will work, but [raster::rotate] is even easier

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>