[R-sig-Geo] Pole centric map with Raster?

[Ben Tupper]

Hello,

I'm learning how to use Raster* objects and ultimately hope to draw rasters from a (north) polar perspective.  I can do something like this using tiling in ggplot2:

https://dl.dropboxusercontent.com/u/8433654/polar-map-with-matrix.png


I have hopes of doing similar using raster, sp and lattice.  But I'm drowning in information as neophytes are prone to do.  Below is a self-contained example...  I'm not getting too far.  It creates two plots:

raster without using projectRaster: https://dl.dropboxusercontent.com/u/8433654/raster_no_project.png

raster using projectRaster: https://dl.dropboxusercontent.com/u/8433654/raster_with_project.png

Neither is what I would like.  How does one get the map centered on the pole and draw a (warped) raster using spplot?

Thanks!
Ben


##### START
library(raster)
library(maps)
library(maptools)
library(rgdal)

# create a Raster from 'volcano'
my_proj <- "+proj=stere +lat_0=90 +lon_0=-30 +x_0=0 +y_0=0" 
r <- raster(volcano,
   xmn = -90, xmx = 0, ymn = 50, ymx = 85, 
   crs = my_proj)
# class       : RasterLayer 
# dimensions  : 87, 61, 5307  (nrow, ncol, ncell)
# resolution  : 1.47541, 0.4022989  (x, y)
# extent      : -90, 0, 50, 85  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=stere +lat_0=90 +lon_0=-30 +x_0=0 +y_0=0 +ellps=WGS84 
# data source : in memory
# names       : layer 
# values      : 94, 195  (min, max)   

pr <- projectRaster(r, crs=my_proj)
# class       : RasterLayer 
# dimensions  : 123, 71, 8733  (nrow, ncol, ncell)
# resolution  : 1.48, 0.402  (x, y)
# extent      : -97.4, 7.68, 42.79, 92.236  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=stere +lat_0=90 +lon_0=-30 +x_0=0 +y_0=0 +ellps=WGS84 
# data source : in memory
# names       : layer 
# values      : 93.88372, 194.8761  (min, max)

# create a coastline and transform to the projection
mp <- map2SpatialLines(map(plot = FALSE, interior = FALSE), proj4string = CRS(my_proj))
mp <- spTransform(mp, CRSobj = CRS(my_proj))
# class       : SpatialLines 
# features    : 2904 
# extent      : -179.9572, 190.2908, -85.44308, 83.57391  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=stere +lat_0=90 +lon_0=-30 +x_0=0 +y_0=0 +ellps=WGS84 

# show the map with the 'unprojected' raster
plot_r <- spplot(r,
   sp.layout = list( sp.lines, mp, first = FALSE)
)

# show the map with the projected raster
plot_pr <- spplot(pr,
   sp.layout = list( sp.lines, mp, first = FALSE)   
)

print(plot_r)
print(plot_pr)
##### END



> sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] rgdal_0.9-2     maptools_0.8-34 maps_2.3-9      raster_2.3-33   sp_1.0-17      

loaded via a namespace (and not attached):
[1] foreign_0.8-63  grid_3.1.2      lattice_0.20-30 tools_3.1.2 


Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org