Calculate average distance to coastline in R

Question

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, ...).

Answer 1

Something like this could work.

Example data

library(terra)
coast <- rnaturalearth::ne_coastline(scale = 110, ret="sv") |> aggregate()
# the below is to make sure it also works in the Pacific Ocean
E <- crop(coast, ext(90, 180, -90, 90))
W <- crop(coast, ext(-180, 0, -90, 90))
coast <- rbind(shift(E, -360), coast, shift(W, 360)) |> aggregate()

# your points
pnts <- cbind(c(90,-30, 20), c(-20,20,10))
vpnts <- vect(pnts, crs=coast)

Solution

directs <- vect(c("LINESTRING (0 180, 0 -180)", # EW 
                  "LINESTRING (-180 -180, 180 180)", # SW-NE
                  "LINESTRING (-180 180, 180 -180)", # NW-SE
                  "LINESTRING (-360 0, 360 0)"), # NS
                       crs=coast)

mdist <- rep(NA, nrow(pnts))

i <- 1
#for (i in 1:nrow(pnts)) {
drs <- shift(directs, pnts[i,1], pnts[i,2])
edr <- erase(drs, coast) |> disagg()
j <- which(relate(edr, buffer(vpnts[i], 100), "intersects"))
d <- distance(as.points(edr[j]), vpnts[i], unit="km")
mdist[i] <- mean(d)
#}

d
[1,] 4671.684
[2,] 5251.762
[3,] 6366.009
[4,] 2238.266
[5,] 5722.547
[6,] 6037.589
[7,] 4306.038
[8,] 3041.791

mdist[i]
[1] 4704.461

Illustrated below. Note that the red lines are not the shortest paths (except in the NS direction).

countries <- rnaturalearth::ne_countries(scale = 110, ret="sv") |> aggregate()
plot(countries, col="light gray", background="azure", border="gray", buffer=FALSE)
lines(drs)
lines(edr[j], col="red", lwd=3)
points(vpnts[i], pch=1)
points(edr[j], col="blue")

# add shortest paths
p <- crds(as.points(edr[j]))
p <- terra::as.lines(cbind(pnts[i,1], p[,1], pnts[i,2], p[,2]), crs=crs(edr)) |> densify(100000)
lines(p, lty=2)

For the third point (on land):

A function to create directions could look like this

directions <- function(d) {
    d <- d * pi/180
    sd <- rep(cos(d), each=2) * c(-360, 360)
    cd <- rep(sin(d), each=2) * c(-180, 180)
    sdcd <- cbind(rep(1:length(d), each=2), sd, cd)
    vect(sdcd, "lines", crs="lonlat")
}

d <- directions(seq(0, 150, 30))
plot(countries, asp=2)
lines(d, col="red")

---

The following gets the last intersections (as requested in a comment).

library(terra)
conts <- rnaturalearth::ne_countries(scale = 110, ret="sv")
af <- conts[conts$continent == "Africa", ]
set.seed(1)
af <- sample(af, 25)
pnt <- cbind(20, 10)

# have separate segments for N and S etc.
directs <- vect(c("LINESTRING (0 180, 0 0)" , "LINESTRING(0 0, 0 -180)", # E and W 
                  "LINESTRING (-180 -180, 0 0)","LINESTRING(0 0, 180 180)", # SW and NE
                  "LINESTRING (-180 180, 0 0)", "LINESTRING(0 0, 180 -180)", # NW-SE
                  "LINESTRING (-360 0, 0 0)", "LINESTRING(0 0, 360 0)"), # N and S
                     crs=af)
values(directs) <- data.frame(ID=1:8)

z <- shift(directs, pnt[,1], pnt[,2]) |> crop(af)
g <- geom(z)
dist <- distance(g[, c("x", "y")], pnt, lonlat=T)
d <- data.frame(geom=g[,1], g[, c("x", "y")], dist=dist)
m <- merge(aggregate(dist ~ geom, max, data = d), d)
m
#  geom      dist          x          y
#1    1 1264026.4  20.000000  21.422728
#2    2 3522543.5  20.000000 -21.845463
#3    3  333993.1  17.858601   7.858601
#4    4 1845315.3  32.000000  22.000000
#5    5 3863310.6  -5.758871  35.758871
#6    6 3105606.6  39.841971  -9.841971
#7    7 3758767.3 -14.299010  10.000000
#8    8 3371388.6  50.761263  10.000000

plot(af, col="light gray", border="gray", lwd=2); lines(z, col="blue", lwd=2); points(z)
points(m[, c("x", "y")], col="red", pch=20, cex=2)