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!

I guess this is a duplicate of: https://stackoverflow.com/questions/51292952/snap-a-point-to-the-closest-point-on-a-line-segment-using-sf , meaning that the solution provided there probably answers your question (although it's still somewhat cumbersome). — I_O, Jul 11 '23 at 21:37
Hmmm. I'm not able to make that solution work in this case. Not sure if it's the package or the shapefile I'm using but the first loop wont even run. — Alyssa C, Jul 11 '23 at 22:01
EDIT: I was able to run but it completely butchers the location of my points. I am trying to snap to coastline in Argentina and the linestring method puts it in Antarctica! — Alyssa C, Jul 11 '23 at 22:08
Well, both starting with "A", that's a first step ;-) See answer please. — I_O, Jul 11 '23 at 22:32

I_O · Answer 1 · 2023-07-11T22:49:29.900

one approach, using the names of your example data:

place points at the desired resolution along the coastline (using st_geod_segmentize as we're working with spherical coordinates):

library(lwgeom)
library(units)

points_along_coastline <- 
  st_geod_segmentize(coastline,
                     max_seg_length = set_units(10, km) ## adapt as fit
                     ) |>
  st_cast('POINT')

retrieve indices of nearest points for your offshore points:

indices_nearest <- st_nearest_feature(points_sf, points_along_coastline)

control:

library(ggplot2)

ggplot() +
  geom_sf(data = coastline) +
  geom_sf(data = points_along_coastline, alpha = .2, size = 1) +
  geom_sf(data = points_sf, col = 'blue') +
  geom_sf(data = points_along_coastline[indices_nearest], col = 'red') +
  coord_sf(xlim = c(-73, -70), ylim = c(-35, -33))

score 2 · Accepted Answer · answered Jul 12 '23 at 11:37

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

Neat. Still baffles me that/how `st_nearest_points` finds orthogonals to the coastline - or does it discretise the coastline under the hood? — I_O, Jul 12 '23 at 14:24
I have little idea... a quick check of the source code https://github.com/r-spatial/sf/blob/main/R/nearest.R#L46 says that it calls a S2 function - and that would get us into a rabbit hole of C++ code that quickly gets above my paygrade. — Jindra Lacko, Jul 12 '23 at 17:59
This solution works great!! super elegant. The only thing I wish it did was to retain the other columns from points_sf, which are metadata associated with the point -- essentially so I can plot the data as if they were the original ones. Any ideas? — Alyssa C, Aug 01 '23 at 20:15
I believe that the number and order of points is not changed, so you should be able to pass the metadata as a vector, perhaps after turning sfc to sf object... — Jindra Lacko, Aug 03 '23 at 07:40

Adjust points to closest point on a coastline in R

2 Answers2