[R-sig-Geo] Pole centric map with Raster?
Ben Tupper
btupper at bigelow.org
Sat Mar 28 00:41:07 CET 2015
Hi,
On Mar 27, 2015, at 5:55 PM, Edzer Pebesma <edzer.pebesma at uni-muenster.de> wrote:
>
>
> 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.
>
Yes, in the ggplot2 example the volcano matrix is transformed into a data frame, 'mat', of [lon,lat,z]. I followed this nice post:
http://www.numbertheory.nl/2011/11/08/drawing-polar-centered-spatial-maps-using-ggplot2/
I think the following replicates what I did before (less the coastline). I don't know how to control the boundaries in ggplot2. I'll try what you suggest with SpatialPolygonsDataFrame - I would prefer to use just one graphics paradigm, such as lattice/spplot, if possible.
Thanks!
Ben
### START
library(ggplot2)
# convert a matrix to data frame
as.tile <- function(x = volcano,
bb = c(ll.lat = 45, ll.lon = -90, ur.lat = 75, ur.lon = -40) ){
if (!is.matrix(x)) stop("input must be a matrix")
d <- dim(x)
xx <- rep(seq(from = bb[['ll.lon']], to = bb[['ur.lon']], length = d[2]), each = d[1])
yy <- rep(seq(from = bb[['ll.lat']], to = bb[['ur.lat']], length = d[1]), d[2])
data.frame(lon = xx, lat = yy, z = as.vector(x))
}
# custom theme without axes and annotations
theme_polar <- function(){
list(
theme_bw(base_size=18),
theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
legend.position = "none"),
labs(x='',y=''))
}
ylim <- c(30,90)
xlim <- c(-180,180)
projection <- "azequidist"
orientation <- c(90, -30, 0)
mat <- as.tile()
ggplot() +
geom_tile(data = mat, aes(x = lon, y = lat, fill = z)) +
coord_map(projection = projection, orientation = orientation, ylim = ylim, xlim = xlim) +
scale_x_continuous(breaks=(-6:6) * 30, expand = c(0,0)) +
scale_y_continuous(breaks = seq(from = ylim[1], to = ylim[2], by = 15), limits = ylim, expand = c(0,0)) +
theme_polar()
#### END
>>
>>
>> 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
>
> _______________________________________________
> 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