Create conical gradient around a point using R

I have a point on a map that I am simulating bird movements from. I'd like to create a raster grid around this point with a conical gradient, similar to the picture below.

The idea is that the birds generally move in one compass direction. Say, for the sake of this example, that the birds generally only move in a northwesterly direction. So - if we just match the picture above as a guide, I'd like to create a raster where:

• Any grid cells at 300° from the central point = white = a raster value of 1;
• Any grid cells at 120° from the central point = black = a raster value of 0;
• Any grid cells at angles in between 120 and 300 are interpolated to be between 0 or 1.

How can I do this in R? The layer of dust and cobwebs in the 'basic trigonometry' part of my brain is thick and I am struggling!

I am imagining something along the lines of:

``````library(sf)
library(terra)
library(dplyr)

theta <- 120 # my angle of interest

pt <- st_point(c(-126, 53)) %>%
st_sfc(crs = 4326) %>%
st_transform(crs = 3005)

# Create 30km buffer

# Rasterize it to 500m resolution
cone <- vect(p30)
cone <- rast(cone, resolution = 500)

# Roughly based off answers here
# https://stackoverflow.com/questions/5300938/calculating-the-position-of-points-in-a-circle
# Generate coordinates for a circle of points around my central point `pt`
# Given a radius `r`, angle (theta) `t`, and circle's center `(h, k)`:
circ_coords <- function(r, t, h, k){
t <- t * pi/180
# Calculate coordinates
x <- r * cos(t) + h
y <- r * sin(t) + k
z <- c(x, y)
names(z) <- c('X', 'Y')
return(z)
}

circle <- lapply(seq(theta, theta+360),
circ_coords,
r = radius, # 30km buffer
h = st_coordinates(pt)[1],
k = st_coordinates(pt)[2])
circle <- data.frame(do.call(rbind, circle))

# Center of 'cone'
center <- cellFromXY(cone, st_coordinates(pt))
# Edges of 'cone'
edges <- cellFromXY(cone, circle)

# Just to visualise...
values(cone) <- 0
cone[center] <- 1
cone[edges] <- 2
plot(cone)

# Interpolate the edge values
# The first record in `circle` represents our angle
# of interest where the value should be 1; the 181st
# record of `circle` is 180° from there, where the
# value should be 0. The last record of `circle` is
# back to our starting point and has a value of 1.
circle\$value <- NA
circle[["value"]][1] <- 1
circle[["value"]][181] <- 0
circle[["value"]][361] <- 1

# Now interpolate the in-betweenies
circle\$value <- zoo::na.approx(circle\$value)

# Reassign the raster values for the edge cells
cone[edges] <- circle\$value

# Visualize
plot(cone)

# Now do angular interpolation??
# Perhaps something similar as to
# https://stackoverflow.com/questions/2708476/rotation-interpolation
# and
# https://stackoverflow.com/questions/44523568/r-angular-distance-weighting-interpolation-function
# ?

# And ideally output a result that looks similar to the picture above.
``````

This is as far as I have gotten - the edges look great, but how can I interpolate the middle??

The angle is independent of radius inside the circle. This is probably easiest done by getting the cell xy coordinates and doing trigonometry over those coordinates, then sticking the values into the raster:

``````makecone <- function(pt, radius, theta, res){
theta = pi*theta/180

## Rasterize it to 500m resolution
cone <- vect(p30)
cone <- rast(cone, resolution = res)
xy = crds(cone)
xy[,1] = xy[,1] - st_coordinates(pt)[1]
xy[,2] = xy[,2] - st_coordinates(pt)[2]

v = (atan2(xy[,1], xy[,2]) + theta)

## scale 0-1
cone[] = (1+sin(v))/2
cone
}
``````

For testing:

``````c1 = makecone(pt, radius, 120, 500)
plot(c1)
plot(c1>0.99)
``````

The second plot is just to make sure the "1" is in the right direction.

• This was excellent, thank you. I have modified the radians conversion slightly to ensure the correct angle in degrees is used: # Convert degrees to radians # 1) Multiply theta by -1 to ensure radians go clockwise (radians are measured counter-clockwise, by default) # 2) Multiply by pi/180 to convert to radians # 3) Subtract 270° to rotate the angle clockwise by 270° (radians by default are measured from the x-axis, i.e. "east", rather than "north" on a compass-rose) `theta <- (-1 * theta) * (pi/180) - (3 * pi/2)`
– srha
Commented Aug 31, 2023 at 21:14