1

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.

enter image description here

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?
3

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.

enter image description here

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)
1
  • 1
    Nice - I didn't have to time to flesh out a full solution like yours but glad you did! – Spacedman Feb 11 '19 at 15:33

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.