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, ...).
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)
terra::extractAlongfor coastlines after projection but without merging to multilinestring. – Chris Commented Jan 31 at 14:44