[R-sig-Geo] working with a rotated grid of SpatialPolygons

Waichler, Scott R Scott.Waichler at pnnl.gov
Thu Nov 14 05:00:18 CET 2013


I have an environmental process model whose domain is a regular grid that is not oriented in the cardinal directions.  I want to show an outline of the active part of the domain, as well as other subgroupings of cells within the domain, on a map figure.  These will be "blocky" polygons oriented at an angle to E-W, N-S.  I assume I should do the rotation on (x,y,z) with rotate3d() of rgl.  If there is a way to rotate sp objects directly, do tell.

library(rgl)
library(sp)
x <- y <- seq(0.5, 9.5, by=1)  # cell centers of a toy grid
xy <- expand.grid(x, y)
xyz <- matrix(data = 0, nrow=nrow(xy), ncol=3)
xyz[,1] <- xy[,1]
xyz[,2] <- xy[,2]
# rotate the grid 45 degrees counterclockwise
xyz.r <- rotate3d(obj=xyz, angle=-pi/4, x=0, y=0, z=1)
# plot the cell centers (not vertices) of the original and rotated grids
plot(xyz[,1], xyz[,2], col="blue", type="p", asp=1, xlim=c(-6, 10), ylim=c(0, 15), main="grid cell centers")
points(xyz.r[,1], xyz.r[,2], col="red", type="p")

# I have learned this for defining polygons based on the original grid:
grd <- GridTopology(cellcentre.offset=c(0.5,0.5), cellsize=c(1,1), cells.dim=c(10,10))
grdp <- as.SpatialPolygons.GridTopology(grd)

How can I efficiently extract the vertices of all the grdp polygons (each one a square) so that I can rotate them?  I think str() and slot() are needed but I don't understand them.  Once I have the rotated vertex coordinates, how do I use them to define the new SpatialPolygons?  Once I have SpatialPolygons, I want to use unionSpatialPolygons(), so I also need to understand how to carry through an attribute for classifying grid cells.

I could easily rotate vertices instead of cell centers, but I wouldn't know how to make grid cell polygons from the matrix of coordinates.

Scott Waichler
Pacific Northwest National Laboratory, Washington, USA
scott.waichler at pnnl.gov



More information about the R-sig-Geo mailing list