Calculate average distance to coastline in R

I'm using R and the terra package for spatial data operations.

For a number of points on Earth, I want to calculate the average distance to a coastline. So far, I have managed to calculate the minimum distance to a coastline:

# libraries
library(tidyverse)
library(sf)
library(terra)
library(rnaturalearth)


# generate a pixel raster at 200 km resolution
rasterPx_200km     <- terra::rast(resolution = 200000,              # 200km x 200 km resolution
                                  xmin = -17367530, xmax=17432470,  # custom extent (without Antarctica)
                                  ymin=-6045743, ymax=7154257,
                                  crs = "EPSG:6933")                # metric WGS84 


# convert into polygons
rasterPoly_200km      <- terra::as.polygons(rasterPx_200km, crs = "epsg:6933")  # raster to poly
rasterPoly_200km$rid  <- 1:nrow(rasterPoly_200km)                               # add raster id



# get center points
centroidPoints          <- terra::centroids(rasterPoly_200km)   # SpatVector      
centroidPoints_df       <- centroidPoints %>% sf::st_as_sf()    # data frame

# get coastline from Natural Earth
coastline         <- rnaturalearth::ne_coastline(scale = 110)                      
coastline_union   <- sf::st_union(coastline)                        # union into MultiLineString
coastline_proj    <- sf::st_transform(coastline_union, crs = 6933)  # same projection as centroidPoints


# calculate shortest distance from coast
minDistanceCoast    <- st_distance(coastline_proj, # coastline as one multiline element
                                   centroidPoints_df, # points
                                   unit="m", 
                                   pairwise = TRUE)   # vector (columns contain distances)

# data frame with distance to coast
minDistanceCoast_df <- minDistanceCoast  %>%
  t()                                    %>%   # transpose to set up columns vs. rows for data frame
  as.data.frame()                        %>%   # make data frame
  rename("minDistanceCoast" = 1)               # rename first column


units(minDistanceCoast_df$minDistanceCoast) <- "km"  # set to kilometres (automatically makes a conversion)


# full data frame with raster IDs
minDistanceCoast_full <- centroidPoints_df %>%
  cbind(minDistanceCoast_df)

The resulting data frame minDistanceCoast_full contains the id of each point (called rid) and the minimum distance from the coastline (without differentiating for land and water pixels). How can I adapt my code to calculate the average distance for each point in centroidPoints_df to the coast line object coastline_proj? The result should be a data frame containing the columns rid and avgDistanceCoast.

My idea would be to generate rays from the points to a number of different directions and measure the distance whenever they touch a coastline; then averaging the distances. It would be good if the number of rays to calculate the average could be set in a function. I'd like to start with 8 distances in the cardinal and intercardinal directions (N, NE, E, ...).

This is a cleaned-up answer that shows how to get the average distance to the coast from a point on land, when traveling along paths of constant bearing.

First a function to get the paths

get_dirs <- function(n) {
    dirs <- vect(cbind(0,0), crs="local") |> 
            buffer(360, quads=round(n/4)) |> crds() |> unique()
    dirs <- cbind(0, dirs[,1], 0, dirs[,2]) |> as.lines()
    values(dirs) <- data.frame(ID=1:nrow(dirs))
    crs(dirs) <- "+proj=longlat"
    dirs
}

And a function to get the average distance

get_dist <- function(land, directions, lonlat) {
    v <- vect(lonlat, crs="lonlat")
    i <- which(relate(v, land, "intersects"))
    if (length(i) == 0) return(NA)
    drs <- shift(directions, lonlat[1], lonlat[2])
    edr <- crop(drs, land[i]) |> disagg()
#   i <- which(relate(edr, buffer(v, 10), "intersects"))
    i <- which(round(distance(edr, v), -1) == 0)
    e <- densify(edr[i], 0.1, flat=TRUE)
    d <- perim(e) / 1000
    mean(d)
}

Example data, here using Africa. Remove the country boundaries, and disaggregate to separate islands/continents

library(terra)
land <- rnaturalearth::ne_countries(scale = 110, ret="sv")
africa <- land[land$continent == "Africa", ] |> aggregate() |> disagg()
# example points
pnts <- cbind(c(-10, 20, 20), c(10, -20, 20))

And a solution using 360 directions

dirs <- get_dirs(360)
sapply(1:nrow(pnts), \(i) get_dist(africa, dirs, pnts[i,,drop=F]))
# [1] 1804.092 2176.476 2766.308

Illustrated for the first point:

v <- vect(pnts[1,,drop=FALSE], crs="lonlat")
drs <- shift(dirs, pnts[1,1], pnts[1,2])
edr <- crop(drs, africa) |> disagg()
i <- which(round(distance(edr, v), -1) == 0)
plot(africa, col="light gray", border=NA, background="azure")
lines(edr, col="blue"); lines(edr[i], col="red")
points(edr[i], cex=.25); points(v)