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

Edzer Pebesma edzer.pebesma at uni-muenster.de
Fri Mar 27 22:55:58 CET 2015



On 03/27/2015 10:14 PM, Tim Appelhans wrote:
> Ben
> here's some code I used for a paper.
> 
> Is this producing what you want?
> 
Cool!

The ggplot2 example of Ben looks as if it plots cells that form a
rectangular grid in long-lat onto a polar stereographic projection. One
could do this by converting the raster to a SpatialPolygonsDataFrame,
rgdal::spTransform this to the new projection, and then plot. It might
be improved by not plotting the cell (polygon) boundaries.

> 
> 
> library("remote")
> 
> ###### EXAMPLE I ##########################################################
> ###########################################################################
> library("rworldmap")
> library("rgdal")
> library("rgeos")
> library("gridExtra")
> library("RColorBrewer")
> 
> ## load data
> data("vdendool")
> data("coastsCoarse")
> 
> ## calculate 4 leading modes
> nh_modes <- eot(x = vdendool, y = NULL, n = 4, reduce.both = FALSE,
>                 standardised = FALSE, verbose = TRUE)
> 
> ## create coastal outlines
> ster <- CRS("+proj=stere +lat_0=90 +lon_0=-45")
> xmin <- -180
> xmax <- 180
> ymin <- 20
> ymax <- 90     # Coordinates for bounding box
> bb <- cbind(x = c(xmin, xmin, xmax, xmax, xmin),
>             y = c(ymin, ymax, ymax, ymin, ymin))    #Create bounding box
> SP <- SpatialPolygons(list(Polygons(list(Polygon(bb)), "1")),
>                       proj4string = CRS(proj4string(coastsCoarse)))
> 
> gI <- gIntersects(coastsCoarse, SP, byid = TRUE)
> out <- vector(mode = "list", length = length(which(gI)))
> ii <- 1
> 
> for (i in seq(along = gI)) if (gI[i]) {
>   out[[ii]] <- gIntersection(coastsCoarse[i, ], SP)
>   row.names(out[[ii]]) <- row.names(coastsCoarse)[i]
>   ii <- ii + 1
> }
> 
> nh_coasts <- do.call("rbind", out)
> nh_coasts_ster <- spTransform(nh_coasts, ster)
> 
> lout <- list("sp.lines", nh_coasts_ster,
>              col = "grey30", grid = TRUE)
> 
> ## define colours
> clrs <- colorRampPalette(rev(brewer.pal(9, "RdBu")))
> 
> ## project modes to polar stereographic CRS
> mode1 <- projectRaster(nh_modes[[1]]@r_predictor, crs = ster)
> 
> spplot(mode1, sp.layout = lout,
>        col.regions = clrs(1000), at = seq(-1, 1, 0.2),
>        par.settings = list(axis.line = list(col = 0)),
>        colorkey = list(height = 0.75, width = 1))
> 
> 
> Tim
> 
> 
> 
> On 27.03.2015 21:20, Ben Tupper wrote:
>> 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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi),  University of Münster,
Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
Journal of Statistical Software:   http://www.jstatsoft.org/
Computers & Geosciences:   http://elsevier.com/locate/cageo/
Spatial Statistics Society http://www.spatialstatistics.info

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20150327/90bd9802/attachment.bin>


More information about the R-sig-Geo mailing list