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

[Barry Rowlingson]

On Thu, Nov 14, 2013 at 4:00 AM, Waichler, Scott R
<Scott.Waichler at pnnl.gov> wrote:
> 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.

Ooh thats a lot of questions. If you are working on a regular grid
rotated at a known angle to your coordinate system, it might be
possible to define a coordinate system where your points are aligned
with the coordinate axes.

I do something like this using an oblique mercator projection, as
explained here:

http://rpubs.com/geospacedman/rotatespatial

The usual mercator projection imagines a cylinder wrapped around the
earth touching at the equator, a transverse mercator has the cylinder
wrapped at right-angles to that, along a line of latitude (and its
180-degree opposite), and an oblique mercator imagines the cylinder
wrapped at any great circle.

Is your data in a real geographic coordinate system, like UTM or is
this simply a model with no base on the ground:?

Barry