raster - Rotating ~90° using Two-Point Equidistant projection with Proj4

I'm trying to rotate a tif raster in R using Proj4's Two-Point Equidistant projection - I think the only way to rotate spatial geometries to custom angles. But I'm finding the operation gets very slow approaching 90°, and in fact fails for that angle (possibly 2 separate issues). This demonstrates:

require(raster)
require(rgdal)
# demo raster
r <- disaggregate(raster(volcano[1:60,1:60]), 4)

extent(r) <- c(-5,5,-5,5)
crs(r) <- "+proj=longlat +datum=WGS84"

# proj4 two-point equidistant rotate transformations
tpeqd_45 <- CRS("+proj=tpeqd +lat_1=2 +lon_1=-2 +lat_2=-2 +lon_2=2 +a=6371000 +units=m +no_defs")
tpeqd_80 <- CRS("+proj=tpeqd +lat_1=1.9 +lon_1=0.35 +lat_2=-1.9 +lon_2=-0.35 +a=6371000 +units=m +no_defs")
tpeqd_90 <- CRS("+proj=tpeqd +lat_1=2 +lon_1=0 +lat_2=-2 +lon_2=-0 +a=6371000 +units=m +no_defs")

system.time(r_45 <- projectRaster(r, crs=tpeqd_45))
   user  system elapsed 

  0.350   0.051   0.404 
par(mfcol = c(1,2)); image(r, asp=T); image(r_45, asp=T)

system.time(r_80 <- projectRaster(r, crs=tpeqd_80))
   user  system elapsed 
  3.857   0.879   4.877
system.time(r_90 <- projectRaster(r, crs=tpeqd_90))
Error in if (value[1] != nrow(x) | value[2] != ncol(x)) { : 
  missing value where TRUE/FALSE needed

In addition: Warning message:
In `dim<-`(`*tmp*`, value = c(nr, nc)) :
  NAs introduced by coercion to integer range

Note the increasing time to process for 80°. For a 10x bigger raster this crashes my Macbook. Does anyone know if this is a problem in the math for this projection and/or if there is any other way to rotate these in R?

Answer

Using an azimuthal projection makes it easy to rotate a raster (or any geographic layer, for that matter). Just pretend the data are in a polar azimuthal projection and change its central longitude by the negative of the angle.

Here, for example, is a bare-bones R function illustrating the procedure. Its arguments are a raster object x, an angle of rotation a (in degrees counterclockwise), and an optional output resolution.

• 
  
  

  The first line creates a copy of x and tells the software to interpret the coordinates in the copy as if they were in an azimuthal equidistant projection (at the north pole) with a reference longitude of zero.

  
  


• 
  

  The second line reprojects the copy to an equidistant azimuthal coordinate system with a different reference longitude.

  
  

The net effect is the desired rotation, as the maps in the figure show. (The resolution was also refined in the rotated images to verify the effect of the resolution argument.)

library(raster)
rotate <- function(x, angle=0, resolution=res(x)) {
  y <- x; crs(y) <- "+proj=aeqd +ellps=sphere +lat_0=90 +lon_0=0"

  projectRaster(y, res=resolution, 
                     crs=paste0("+proj=aeqd +ellps=sphere +lat_0=90 +lon_0=", -angle))
}

As an example, I created a raster and rotated it by various angles around its (0,0) coordinate (at its center).

Notice that it is irrelevant what coordinate system is actually used in the input raster. However, the output will be tagged with this artificial rotated azimuthal equidistant system. You can always re-specify it if you want, but in most applications this likely would not be meaningful.

The only limitation to this method is that coordinates in the input raster must be valid coordinates in the azimuthal system. This will rarely be a problem, but if it is, first shrink the coordinates, perform the rotation, and then expand the coordinates.

Here is the R code to create the maps.

x <- raster(matrix(1:(15*25), nrow=15), xmn=-1000, xmx=1000, ymn=-1000, ymx=1000)

plot(x, main="Original")
for (a in c(45, 90, 150))
  plot(rotate(x, a, resolution=10), main=paste("Rotated by", format(a, ndigits=0)))

---

This method, although convenient to code, is not blazingly fast: for some purposes (where resampling is not needed, for instance) you can do much better by applying a rotation matrix directly to the underlying data. But it's not too slow, either. On my system a rotation of a 600 by 1000 raster takes under 3 seconds:

> x <- raster(matrix(1:(600*1000), nrow=600), xmn=-1000, xmx=1000, ymn=-1000, ymx=1000)
> system.time(rotate(x, 33))

user  system elapsed 

2.73    0.17    2.90