3

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.

cone gradient

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
radius <- 30000 # my radius length of interest (in meters)

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

# Create 30km buffer
p30 <- st_buffer(pt, radius)

# 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){
  # Convert degrees to radians
  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 edges of the circle are interpolated, but not the middle!

1 Answer 1

3

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
    p30 <- st_buffer(pt, radius)
    
    ## 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 = mask(cone, vect(p30))
    cone
}

For testing:

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

enter image description here

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

1
  • 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
    Aug 31, 2023 at 21:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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