[R-sig-Geo] Pole centric map with Raster?
Ben Tupper
btupper at bigelow.org
Wed Apr 1 14:48:32 CEST 2015
Hello again,
Thanks again for your help - it really gave me the push I needed. Below is a stripped down version of your script that projects the datasets::volcano matrix over the north pole. In case it is of some use to others, I have added some notes one what I think is happening. If you run the script you'll see that I am still fiddling with the graticules (reaching 90N and the one missing at 180).
Thanks also for making me aware of the remote package - that looks really useful.
Cheers and thanks,
Ben
### START
library(raster)
library(rgdal)
library(rgeos)
library(rworldmap)
library(RColorBrewer)
data("coastsCoarse", package = 'rworldmap')
# the crs of coastsCoarse has a leading space, this trims it
crs(coastsCoarse) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
# the desrired projection
ster <- "+proj=stere +lat_0=90 +lon_0=-30"
# map view limits
xmn <- -180
xmx <- 180
ymn <- 30
ymx <- 90
# create a raster from datasets::volcano
r <- raster::raster(volcano,
xmn = xmn, xmx = xmx, ymn = ymn, ymx = ymx,
crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
# if the raster extends into opposite hemisphere it helps to crop it
#r <- crop(b, extent(xmn, xmx, ymn, ymx))
# create a polygon used for trimming the map segments
# bounding box from lower left clockwise
bb <- cbind(x = c(xmn, xmn, xmx, xmx, xmn),
y = c(ymn, ymx, ymx, ymn, ymn))
# convert to Spatial
SP <- sp::SpatialPolygons(
list(sp::Polygons(
list(sp::Polygon(bb)), "1")),
proj4string = sp::CRS(sp::proj4string(coastsCoarse))
)
# compute where map bounding box and map segements overlap
# I think this amounts to cropping the map segments
gI <- rgeos::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]] <- rgeos::gIntersection(coastsCoarse[i, ], SP)
row.names(out[[ii]]) <- row.names(coastsCoarse)[i]
ii <- ii + 1
}
}
nh_coasts <- do.call("rbind", out)
# project the coastline map segments
nh_coasts_ster <- sp::spTransform(nh_coasts, CRS(ster))
# make graticules and project
gr <- sp::gridlines(nh_coasts,
easts = seq(xmn, xmx, by = 20),
norths = seq(ymn, ymx, by = 10),
ndiscr = length(xmn:xmx))
pgr <- sp::spTransform(gr, CRS(ster))
# project the raster volcano data
projr <- projectRaster(r, crs = ster)
# sp.layout elements
lout <- list(
list("sp.lines", nh_coasts_ster, col = "grey30", grid = TRUE),
list("sp.lines", pgr, col = 'grey50', lty = 'dashed'))
# pretty colors function from RColorBrewer
clrs <- colorRampPalette(brewer.pal(9, "YlOrRd"))
# compute steps for coloring
data_range <- range(values(projr), na.rm = TRUE)
steps <- 10
polar_map <- sp::spplot(projr, sp.layout = lout,
col.regions = clrs(steps),
at = seq(data_range[1], data_range[2], length = steps),
par.settings = list(axis.line = list(col = 1)),
colorkey = list(height = 0.75, width = 1))
print(polar_map)
### END
On Mar 27, 2015, at 5:14 PM, Tim Appelhans <tim.appelhans at gmail.com> wrote:
> 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/
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org
More information about the R-sig-Geo
mailing list