3

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)

enter image description here

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?

  • 1
    Possible workaround: use oblique Mercator instead, center point and alpha but set gamma to 0. omerc also supports definition by two points, but may reorient so north is "up" or use gamma to set that. There's also a no_rot parameter that might work. – mkennedy Mar 3 '16 at 0:30
  • Nice suggestion. I've tried with e.g. +proj=omerc +lat_0=0 +lonc=0 +alpha=45, with various alphas and +no_rot argument but also no luck. – geotheory Mar 3 '16 at 0:54
  • 1
    Maybe try switching to +gamma? – mkennedy Mar 3 '16 at 1:14
  • Yes that does start to work. Unclear how the params alpha and gamma interrelate though: I seem to have similar issues to this question. I'm working through it.. – geotheory Mar 3 '16 at 11:56
  • Alpha is applied near the beginning within the mathematics of the projection. Gamma is applied at the end to re-orient the rotated cylinder--often used to have geodetic/true north "up" at the center of the projection. – mkennedy Mar 3 '16 at 15:47
4

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).

Figure

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 
  • Awesome thanks @whuber. What would you recommend as the best reference for Proj4? Preferably something outlining all valid arguments for a particular projection.. I've been trawling through the official manual and remotesensing but not found anything particularly helpful. – geotheory Mar 4 '16 at 17:25
  • Great answer! You can also use rasterImage with its angle and interpolation arguments, but you need to correct for the rotation centre (bottom left corner). – mdsumner Apr 7 '16 at 0:17
  • @whuber, thank you for your answer; however, I have question, so I am trying to rotate a raster file with a defined proj4string(r) <- "+proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0", do I need to change something in your rotate function? I am interested in keeping the same proj4string. – Perro Sep 26 '18 at 16:27
  • @Perro You shouldn't need to change anything. The projected raster will not have an appropriate projection string associated with it, though, nor should it--after all, after this rotation it won't be in the same projection as the input raster. It looks like the comments after the figure fully explain this. – whuber Sep 26 '18 at 16:30
  • @whuber thank your your response. I know understand what you are saying about the input and out projection system. Based on that, I was wondering if it is possible to rotate this image that I have (lets say 90deg) without affecting its projection system? – Perro Sep 26 '18 at 16:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.