Adjust points to closest point on a coastline in R

Question

I'm trying to find the closest point on a coastline to a given set of lat-long points in R. I have an sf object representing the coastline and a data frame with the lat-long points.

Here's an example of the data:

library(sf)
library(rnaturalearth)

# Load world coastline data
coastline <- ne_coastline(scale = "small", returnclass = "sf")

# Lat-long points
points <- data.frame(
  longitude = c(-71.5374, -72.1234, -70.9876),
  latitude = c(-33.865, -34.567, -33.456)
)
points_sf <- st_as_sf(points, coords = c("longitude", "latitude"), crs = 4326)

I want to find the closest point on the coastline for each point in the points_sf object, and create a new dataset of those coastal points that I can re-join to my old dataset and plot.

I have tried various methods, including using the distGeo function from the geosphere package and the st_nearest_feature function from the sf package. However, I'm encountering issues with incorrect dimensions or unexpected results.

Here's what I tried so far but keep encountering problems with the distances line.

library(sf)
library(geosphere)

# Calculate distances between each point and the coastline
distances <- distGeo(st_coordinates(points_sf), st_coordinates(coastline))

# Find the index of the nearest coastline point for each point
nearest_indices <- apply(distances, 1, which.min)

# Get the closest points on the coastline
closest_points <- st_coordinates(coastline[nearest_indices, ])

# Create an sf object with the closest points
closest_points_sf <- st_as_sf(data.frame(geometry = st_point(closest_points)),
                              crs = st_crs(coastline)) 

Any help you can offer would be amazing! I have spent a long time on this. Tahnk you!

Answer 1

The answer by I_O building on segmentized coastline is interesting and ingenious; for completeness sake I propose an alternative approach building on sf::st_nearest_points(). It will be somewhat more precise (given that it does not use sampling) and avoids dependency on {lwgeom}, making do with {sf} alone.

What this piece of code does is:

• unions all the world's coastlines to a single object (otherwise the calculation would be done by feature = many times)
• calculates the connection of origins and coast as lines
• casts the said lines to points (origin + destination)
• subsets the object to get odd numbered points (i.e. filter out the blue origins)

The outcome is otherwise hard to distinguish from the answer by I_O.

points_on_coastline <- st_union(coastline) %>% # one coastline to rule them all!
  st_nearest_points(points_sf) %>% # lines from points to coast
  st_cast("POINT") %>% # convert lines to points / origin, end
  .[seq(1, 2* nrow(points_sf), 2)] # just the coastal parts, thank you!
  
# plot borrowed from earlier answer by I_O
library(ggplot2)

ggplot() +
  geom_sf(data = coastline) +
  geom_sf(data = points_sf, col = 'blue') +
  geom_sf(data = points_on_coastline, col = "red") +
  coord_sf(xlim = c(-73, -70), ylim = c(-35, -33))