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, ...).
Solution
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)
- How to use purrr:map() and rlang to emulate a pipe chain
- magick annotate picture with arrows
- Import png files and convert to animation(.mp4) in R
- Convert a column with text files into separate images in R
- ggplot2 x-axis with many hours for each of many days. Is there a way to span the dates in the x-axis over the hours for that day?
- R odbc::odbcListDrivers() does not list dirver in /opt/homebrew/etc/odbcinst.ini
- Change size and aspect ratio without distortion
- How can I get an R environment via Sys.getenv() with GitHub Actions using secrets?
- ggplot geom_point color based on both x and y axis values
- Perform a random binomial draw for each row in R without rowwise()
- Scrape the university name (in QS World University Rankings website) with R
- ggplot for linear-log regression model?
- How to join (merge) data frames (inner, outer, left, right)
- Avoid rescaling while binning using scale_*_steps
- How can I mock a function globally using testthat?
- R - could not find function "cld"
- Create multiple lagged variables with different offsets
- Expanding dataframe to include non existing values
- Split string to columns based on paragraph ending from ocr'd image
- from magick-image to rasterBrick
- How to remove repeated elements in a vector, similar to 'set' in Python
- Rename multiple variables at once using dplyr
- Calculate average distance to coastline in R
- Reading large multi-part table from file and combing its parts into one tibble
- Processing multiple images with Magick (in R) with transformations
- R: Convert/Read 3D Matrix into a 'magick' object and vice versa
- Error using magick R to import PDF
- Method in R to crop whitespace on svg file
- Read table from a website into R Studio and create a dataframe with the info
- Perspective transformation using R and magick