How to make a raster gradient in any direction in R

Question

How can I create a raster, where the values increase in a certain direction?

I know how to do the "main" directions with full 90 degree steps, as illustrated in the image. But I would need some help to do lets say degree 27, or 154?

Is there a function for that?

My idea was to always create the "North" gradient raster, rotate the input raster to "face Norh", crop/mask the gradient and assign values to the input raster and rotate back. But rotating a raster also doesn't seem to be too trivial.

Anyone some suggesions or magical functions I am not aware of?

Below is some code to produce the plot and play around.

library(raster)
library(sp)
# library(maptools)
Polygon1 <- Polygon(rbind(c(4498482, 2668272), c(4498482, 2669343),
                          c(4499991, 2669343), c(4499991, 2668272)))
Polygon1 <- Polygons(list(Polygon1), 1)
Polygon1 <- SpatialPolygons(list(Polygon1))

## make raster
r <- raster(extent(Polygon1), ncol = 180, nrow = 180)
polyraster <- rasterize(Polygon1, r)
values(polyraster) <- 1

par(mfrow = c(2, 2))

## Make Raster gradient for direction 0
polymat <- as.matrix(polyraster)
polymat <- apply(polymat, 1, function(x) seq(1, ncol(polymat), 1)) 
polymat <- raster(polymat,
                  xmn = extent(polyraster)[1],
                  xmx = extent(polyraster)[2],
                  ymn = extent(polyraster)[3],
                  ymx = extent(polyraster)[4])
plot(polymat, main="North 0°")

## Make Raster gradient for direction 180
polymat <- as.matrix(polyraster)
polymat <- apply(polymat, 1, function(x) rev(seq(1, ncol(polymat), 1))) 
polymat <- raster(polymat,
                  xmn = extent(polyraster)[1],
                  xmx = extent(polyraster)[2],
                  ymn = extent(polyraster)[3],
                  ymx = extent(polyraster)[4])
plot(polymat, main="South 180°")

## Make Raster gradient for direction 270
polymat <- as.matrix(polyraster)
polymat <- apply(polymat, 1, function(x) seq(1, ncol(polymat), 1)) 
polymat <- raster(t(polymat),
                  xmn = extent(polyraster)[1],
                  xmx = extent(polyraster)[2],
                  ymn = extent(polyraster)[3],
                  ymx = extent(polyraster)[4])
plot(polymat, main="East 270°")

## Make Raster gradient for direction 90
polymat <- as.matrix(polyraster)
polymat <- apply(polymat, 1, function(x) rev(seq(1, ncol(polymat), 1))) 
polymat <- raster(t(polymat),
                  xmn = extent(polyraster)[1],
                  xmx = extent(polyraster)[2],
                  ymn = extent(polyraster)[3],
                  ymx = extent(polyraster)[4])
plot(polymat, main="West 90°")



## Rotate Polygon with maptools::elide
poly_rotate <- maptools::elide(Polygon1, rotate = 27)
r <- raster(extent(poly_rotate), ncol = 180, nrow = 180)
polyraster <- rasterize(Polygon1, r)
values(polyraster) <- 1
polymat <- as.matrix(polyraster)
polymat <- apply(polymat, 1, function(x) seq(1, ncol(polymat), 1)) 
polymat <- raster(polymat,
                  xmn = extent(polyraster)[1],
                  xmx = extent(polyraster)[2],
                  ymn = extent(polyraster)[3],
                  ymx = extent(polyraster)[4])
plot(polymat, main="North 0°")
plot(poly_rotate, add=T)
poly_rot_clip <- crop(polymat, extent(poly_rotate))
poly_rot_clip <- mask(poly_rot_clip, poly_rotate)
plot(poly_rot_clip)
# How to rotate back?

Answer 1

A linear combination of X and Y coordinates with scale factors can create a slope in any direction. For example, using a grid from -1 to +1 of length 100 in both directions:

> x = seq(-1,1,len=100)
> y = seq(-1,1,len=100)
> a=1;b=2
> r12 = raster(outer(x,y,function(x,y){a*x+b*y}))

produces a slope dependent on the ratio of a/b, something like:

> 180*atan2(b,a)/pi
[1] 63.43495

in degrees from the horizontal.

Answer 2

Thats the resulting functions I came up with, if anyone needs something similar. All credit goes to @Spacedman!

## Data ##############
library(raster)
library(sp)
Polygon1 <- Polygon(rbind(c(4498482, 2668272), c(4498482, 2669343),
                          c(4499991, 2669343), c(4499991, 2668272)))
Polygon1 <- Polygons(list(Polygon1), 1)
Polygon1 <- SpatialPolygons(list(Polygon1))

## Funcs ##############
inv_atan2 <- function(angle) {
  ## deg to rad
  radian = angle * pi / 180 
  ## invert atan2
  x1 = cos(radian)
  y1 = sin(radian)
  return(cbind(x1, y1))
}
poly_gradient <- function(shape, angle) {
  a = inv_atan2(angle)[1]
  b = inv_atan2(angle)[2]
  x = seq(-1, 1, len = 100)
  y = seq(-1, 1, len = 100)
  ras = raster(outer(x, y, function(x,y) {a * x + b * y}))
  extent(ras) <- extent(shape)
  projection(ras) <- proj4string(shape)
  invisible(ras)
}

## Tests ##############
r12 <- poly_gradient(Polygon1, angle = 189)
plot(Polygon1)
plot(r12, add = T)