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

[Tim Appelhans]

Ben
here's some code I used for a paper.

Is this producing what you want?



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

-- 
#####################################
Tim Appelhans
Department of Geography
Environmental Informatics
Philipps Universität Marburg
Deutschhausstraße 12
35032 Marburg (Paketpost: 35037 Marburg)
Germany

Tel +49 (0) 6421 28-25957

http://environmentalinformatics-marburg.de/