Rotating Spatial Data

Outline

This report shows how to effectively rotate spatial data by an angle. This is motivated by the case where you have an analysis technique that only works on regular grids, and cannot be told to just ignore points outside the polygon by marking them as 'NA' or other missing data indicator.

The usual way to work would be to compute over the full grid and then crop to the polygon of interest, but this results in wasted computer time for the grid cells outside the polygon. Rotating a polygon to be more axis-aligned can significantly decrease the number of grid cells outside the polygon and hence shorten execution times.

The downside of this is that you end up with a grid that is not aligned to the original projection, but that may not be a problem. However using a projection for the rotation means reverse-projections can always be used to back-transform or project rasters to the original basis.

Note this is only good for small regions, as it depends on the cartesian nature of a cylindrical projection at close distances to the effective equator (not the real equator) of the projection. We can define this effective equator to go through the middle data.

Worked Example

First we'll get the 32nd region in the Scottish regions data supplied with the rgdal package.

require(sp)

## Loading required package: sp

require(rgdal)

## Loading required package: rgdal
## rgdal: version: 0.8-11, (SVN revision 479M)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.0, released 2013/04/24
## Path to GDAL shared files: /usr/share/gdal/1.10
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)

dsn <- system.file("vectors", package = "rgdal")[1]
scot <- readOGR(dsn = dsn, layer = "scot_BNG")

## OGR data source with driver: ESRI Shapefile 
## Source: "/nobackup/rowlings/RLibrary/R/x86_64-pc-linux-gnu-library/3.0/rgdal/vectors", layer: "scot_BNG"
## with 56 features and 13 fields
## Feature type: wkbPolygon with 2 dimensions

region = scot[32, ]
plot(region)
title(region$NAME)

This region is quite diagonal. A regular grid over its bounding box is going to be mostly outside the polygon.

Let's first create a function that returns a 1km grid of points over the bounding box:

make1kGrid = function(region) {
    b = bbox(region)
    xs = seq(b["x", "min"], b["x", "max"], by = 1000)
    ys = seq(b["y", "min"], b["y", "max"], by = 1000)
    pts = expand.grid(xs, ys)

    coordinates(pts) = pts
    proj4string(pts) = proj4string(region)
    return(as(pts, "SpatialPoints"))
}

Now plot the region, do a point-in-polygon test and plot the points coloured by their state:

plot(region)
pts = make1kGrid(region)
op = over(pts, as(region, "SpatialPolygons"))

points(pts, col = 1 + is.na(op))

In this case we have 3021 points in our grid and only 992 are in the polygon. We could be wasting computational time on 2029 points.

To improve this ratio we can rotate the polygon to better align with the grid axes. The elide function in the maptools package can rotate polygons but does it naively, and doesn't use map projections that could then be applied to other spatial data.

If we transform our data using an oblique mercator projection, we can effectively rotate the points using a coordinate reference system specified by a PROJ.4 string. The following function builds such a string using the centre of the bounding box of a spatial data object to define the location of the oblique projection's effective equator and the alpha parameter to define the tilt.

rotateProj = function(spobj, angle) {
    # get bounding box as spatial points object
    boxpts = SpatialPoints(t(bbox(spobj)), proj4string = CRS(proj4string(spobj)))
    # convert to lat-long
    boxLL = bbox(spTransform(boxpts, CRS("+init=epsg:4326")))
    # find the centre
    llc = apply(boxLL, 1, mean)
    # construct the proj4 string
    prj = paste0("+proj=omerc +lat_0=", llc[2], " +lonc=", llc[1], " +alpha=", 
        angle, " +gamma=0.0 +k=1.000000 +x_0=0.000 +y_0=0.000 +ellps=WGS84 +units=m ")
    # return as a CRS:
    CRS(prj)
}

Now we can transform our region. Eyeballing and trial and error gives us a rotation angle of -61 degrees.

# get the proj4 string
regionProj = rotateProj(region, -61)

# transform
regionR = spTransform(region, regionProj)


# get a 1km grid
ptsR = make1kGrid(regionR)

plot(regionR)
opR = over(ptsR, as(regionR, "SpatialPolygons"))
points(ptsR, col = 1 + is.na(opR))

Now we have only 2100 points in our grid (921 fewer) but still roughly the same number of points (1004) in the polygon. In this case we are now only wasting time on 1096 points.

You can now use that projection string to transform points or warp rasters.