Snapping points outside a polygon to the polygons in R

Question

I have records of where certain species breed on mainland New Zealand and I want to model where these tend to be. However, partly due to some projection errors, as well as potential measurement errors, some of my sites are out in the ocean, rather than on land which should be impossible. As such, I need to snap these points to the closest segment of the closest polygon (there are multiple polygons representing different islands).

To do this I am working in R and I am using the sf package; I would like to stay in this environment, as the rest of my code is written to work with these objects.

Here is a bit of setup:

library(tidyverse)
library(sf)

#This loads a shapefile which I would like to use to snap my points to. I have hosted the files on github but I couldn't figure out how to call them directly, as the code I usually use for this didn't work with the read_sf() function. See the commented out version below and if you see what I did wrong I would happily edit this to work! Otherwise the shapefile is accessible here:  https://data.linz.govt.nz/layer/51153-nz-coastlines-and-islands-polygons-topo-150k/

#coastline <- read_sf(dsn = "https://raw.github.com/AndreMBellve/Coastline/tree/master/coastlineClean", layer = "nz-coastlines-and-islands-polygons-topo-150k")

#Reading in and transforming the CRS to be consistent with my coordinates for my sites.
coastline <- st_transform(read_sf("./data/Coastline/nz-coastlines-and-islands-polygons-topo-150k.shp"), crs = '+init=EPSG:2193')

#A sample of my dataframe
SeaCoor.df <- structure(list(Long = c(173.830876, 171.382936, 171.365341, 171.337832, 
171.340413, 171.36002, 171.366843, 171.345085, 171.325994, 173.07732, 
173.07732, 173.07732, 173.07732, 173.07732, 173.07732, 173.07732, 
172.68566, 172.77134, 173.13035, 172.695354, 173.05046, 172.86761, 
172.82686, 172.77778, 172.76675, 172.80786, 172.794857, 172.80966, 
172.85528, 172.87878, 172.87541, 172.90811, 172.90665, 172.9206, 
172.93519, 172.94691, 172.98296, 173.00158, 173.02016, 173.04214, 
173.08939, 173.10469, 173.10458, 173.10685, 173.10587, 173.10943, 
173.11913, 173.12599, 173.12874, 173.13035, 173.11956, 173.10934, 
173.11379, 173.10522, 173.0981, 173.10134, 173.08651, 173.07842, 
173.06943, 173.06205, 173.05389, 173.03132, 173.02989, 173.02282, 
173.01718, 172.99699, 172.9915, 172.98369, 172.96898, 172.96597, 
172.96412, 172.96034, 172.91313, 172.94635, 172.94729, 172.93802, 
172.92575, 172.89502, 172.88644, 172.87614, 172.8691, 172.86017, 
172.85726, 172.85082, 172.83752, 172.83202, 172.82902, 172.82434, 
172.82189, 172.81341, 172.82082, 172.80172, 172.79516, 172.78851, 
172.77855, 172.7828, 172.77065, 172.76014, 172.75593, 172.74937, 
172.74216, 172.74018, 172.73997, 172.74027, 172.73581, 172.72967, 
172.72121, 172.77666, 172.68965, 173.0599, 169.30167, 169.26356, 
169.2309, 168.88449, 168.88449, 168.88449, 168.88449, 168.88449, 
168.62211, 168.51365, 167.99494, 167.8171, 170.979492, 169.611599, 
167.8462, 167.99277, 168.56881, 168.53276, 168.56881, 168.56881, 
168.59681, 168.38814, 166.86198, 166.89506, 166.86198, 166.86198, 
166.89506, 166.89506, 166.89506, 166.89506, 166.89506, 167.63768, 
167.63768, 167.63768, 167.63768, 167.63768, 167.63768, 167.63768, 
167.59288, 167.68446, 167.63768, 168.206138, 168.230986, 168.211588, 
168.239311, 168.206138, 168.245191, 168.264889, 168.239311, 168.206138, 
168.211588, 168.264889, 168.239311, 167.7906, 167.8441, 167.88954, 
167.9883, 168.0017, 168.0195, 168.12068, 168.12977, 168.17831, 
168.17762, 168.2126, 168.21522, 168.2199, 168.0459, 167.66217, 
167.7041, 167.65586, 167.40461, 167.55494, 167.43718, 167.4162, 
167.43718, 167.43718, 167.40461, 167.3956, 167.43718, 167.40461, 
167.55494, 167.38833, 176.659641, 176.494675, 176.523342, 176.656036, 
176.414337, 176.402149, 176.347733, 176.326962, 176.704702, 176.240273, 
176.48798, 176.571836, 176.80315, 176.626511, 176.418972, 176.269455, 
176.347733, 176.626511, 176.566772, 176.704702, 176.240273, 176.48798, 
176.544199, 176.701698, 176.719551, 176.847954, 176.892586, 176.326962, 
176.18259, 176.24165, 176.17641, 176.22161, 176.18093, 176.18093, 
176.18093, 176.18093, 176.29623, 176.29623, 176.28079, 176.28079, 
176.31804, 176.29623, 176.31804, 176.28079, 176.33658, 176.28079, 
176.28092, 176.28092, 176.28092, 175.83443, 175.83443, 175.83443, 
176.01399, 176.01399, 176.01399, 176.01399, 176.80761, 176.81422, 
179.024443, 179.024443, 179.024443, 179.024443, 179.024443, 178.78107, 
178.78107, 178.78107, 178.78107, 178.78107, 178.78107, 166.61003, 
166.58097, 166.58097, 166.61964, 166.57329, 166.58097, 166.57329, 
166.61964, 166.60376, 166.58097, 166.61964, 166.58097, 166.58097, 
166.50304, 166.50149, 166.49343, 166.49343, 166.50149, 166.27482, 
166.27482, 166.27344, 166.27482, 166.30314, 166.2422, 166.27482, 
166.27482, 166.27482, 166.30314, 166.27482, 166.27344, 166.27482, 
166.27344, 166.26074, 166.14821, 166.14821, 165.89064, 166.12633, 
166.20049, 166.22109, 166.26074), Lat = c(-42.247994, -42.032578, 
-42.048609, -42.09203, -42.100929, -42.070531, -42.032992, -42.084387, 
-42.114524, -43.06163, -43.06163, -43.06163, -43.06163, -43.06163, 
-43.06163, -43.06163, -43.6349, -43.85663, -43.76366, -43.629541, 
-43.85687, -43.89607, -43.62676, -43.57436, -43.85641, -43.58683, 
-43.5937, -43.61281, -43.60349, -43.60775, -43.61915, -43.62452, 
-43.63117, -43.63896, -43.62694, -43.63263, -43.63862, -43.65639, 
-43.65446, -43.65135, -43.68795, -43.6959, -43.7004, -43.70626, 
-43.71426, -43.71885, -43.78783, -43.74152, -43.75343, -43.76366, 
-43.77574, -43.78454, -43.79511, -43.80406, -43.82165, -43.83277, 
-43.83558, -43.85366, -43.84623, -43.84864, -43.85526, -43.87302, 
-43.87529, -43.87432, -43.88061, -43.88305, -43.87587, -43.88651, 
-43.88791, -43.87281, -43.87234, -43.87036, -43.85421, -43.89171, 
-43.89721, -43.90012, -43.89894, -43.89251, -43.89647, -43.89573, 
-43.89183, -43.8893, -43.88162, -43.88565, -43.88991, -43.88181, 
-43.87977, -43.87896, -43.87866, -43.87934, -43.87324, -43.86854, 
-43.8665, -43.86619, -43.86415, -43.85653, -43.85947, -43.85325, 
-43.8509, -43.85458, -43.84691, -43.84217, -43.83666, -43.83437, 
-43.83409, -43.83258, -43.82908, -43.57181, -43.62878, -43.83855, 
-43.68369, -43.6961, -43.71268, -43.86064, -43.86064, -43.86064, 
-43.86064, -43.86064, -43.96016, -44.00356, -44.32841, -44.49668, 
-45.109398, -46.541206, -46.45796, -46.40591, -46.76869, -46.73782, 
-46.76869, -46.76869, -46.7501, -46.79583, -46.56985, -46.57322, 
-46.56985, -46.56985, -46.57322, -46.57322, -46.57322, -46.57322, 
-46.57322, -47.78846, -47.78846, -47.78846, -47.78846, -47.78846, 
-47.78846, -47.78846, -46.75172, -46.77682, -47.78846, -46.912993, 
-46.869118, -46.848165, -46.907891, -46.912993, -46.82374, -46.866067, 
-46.907891, -46.912993, -46.848165, -46.866067, -46.907891, -46.6915, 
-46.7028, -46.69362, -46.7695, -46.7888, -46.8083, -46.93475, 
-46.95386, -46.94205, -46.97392, -47.0416, -47.06059, -47.09674, 
-47.1198, -47.23868, -47.2613, -47.22606, -47.24073, -47.11769, 
-47.22051, -47.13956, -47.22051, -47.22051, -47.24073, -47.26532, 
-47.22051, -47.24073, -47.11769, -47.21402, -44.068655, -44.089589, 
-44.101147, -44.044167, -43.752993, -43.755411, -43.735879, -44.04935, 
-43.754233, -43.774317, -43.896408, -43.943395, -43.745305, -43.697417, 
-44.06292, -43.730918, -43.735879, -43.697417, -43.94253, -43.754233, 
-43.774317, -43.896408, -43.899067, -43.820409, -43.830255, -43.844432, 
-43.823877, -44.04935, -44.22724, -44.35632, -44.25215, -44.22698, 
-44.35017, -44.35017, -44.35017, -44.35017, -44.26622, -44.26622, 
-44.24003, -44.24003, -44.27802, -44.26622, -44.27802, -44.24003, 
-44.28663, -44.24003, -44.43272, -44.43272, -44.43272, -43.96254, 
-43.96254, -43.96254, -44.22317, -44.22317, -44.22317, -44.22317, 
-43.56422, -43.56379, -47.7512, -47.7512, -47.7512, -47.7512, 
-47.7512, -49.68453, -49.68453, -49.68453, -49.68453, -49.68453, 
-49.68453, -48.01981, -48.03035, -48.03035, -48.04438, -48.03566, 
-48.03035, -48.03566, -48.04438, -48.01074, -48.03035, -48.04438, 
-48.03035, -48.03035, -48.05387, -48.05528, -48.05866, -48.05866, 
-48.05528, -50.52571, -50.52571, -50.52609, -50.52571, -50.52958, 
-50.51534, -50.52571, -50.52571, -50.52571, -50.52958, -50.52571, 
-50.52609, -50.52571, -50.52609, -50.55347, -50.82306, -50.82306, 
-50.8312, -50.53034, -50.52696, -50.53133, -50.55347)), row.names = 1:301, class = "data.frame")

#Re-projecting for consistency with other rasters and making it a multipoint object for sf
SeaCoor.df %>%
  st_as_sf(coords = c("Long", "Lat")) %>%
  st_cast("MULTIPOINT") %>%
  st_set_crs(4326) %>% 
  st_transform('+init=EPSG:2193') -> SeaCoor.mp 

As far as I can see, everything up to this point works and is fine (although any advice on how it can be improved is welcome!). I plotted my data to check that it is doing what I think - here is the plot prior to snapping:

#Before snapping points
ggplot() + 
  geom_sf(data = coastline)  + 
  geom_sf(data = SeaCoor.mp) 

Now I snap my points:

#Snapping (I know the tolerance is huge)
newPoints <- st_snap(SeaCoor.mp, coastline, tolerance = 1000) 

ggplot() + 
  geom_sf(data = coastline)  + 
  geom_sf(data = newPoints) 

And that is where the problems seem to be occurring.

1. The points are snapping from inside my polygon, as well as outside - I know I didn't code it to exclude the inside, but that is because I don't know how! How do I only snap points that are outside any of the polygons boundaries (i.e. in the sea)?

2. The st_snap() appears to be snapping my points to specific points along the polygons boundary. I think this is because it is snapping them to 'beginnings' and/or 'ends' of polygon lines(?). How do I make it just snap to the nearest segment of the polygon, rather than beginning or ends?

___________________________________UPDATE____________________________________

I have tried the suggestions of @Tim Assal but it does not seem to have worked. Here is what I have tried and you can see that the same issues are still occurring.

# SnappingPoints ----------------------------------------------

#Creating a logical vector to subset the 'outside' points
outside <- sapply(st_intersects(SeaCoor.mp, coastline),function(x){length(x)==0})

#Adding more points to the coastline raster (points every 10 m)
coastline <- st_segmentize(st_cast(coastline, "MULTILINESTRING"), units::set_units(10, m))

#Snapping at sea points
SeaCoor.mp[outside,] <- st_snap(SeaCoor.mp[outside,], coastline, tolerance = 70.71) #70.71 is 50 by root 2 - the grid cells from my rasters are 50 × 50m

#This still appears to be snapping points from quite a distance. This does not appear to be an issue of units be used as the CRS that SeaCoor.mp is stored in is NZGD2000 which appears to be classed in meters.

#Checking post fix points
ggplot() + 
  geom_sf(data = coastline)  + 
  geom_sf(data = SeaCoor.mp) 

The points still appear to be snapping from a large distance away when they should only be moving 10 meters at the maximum! The ones near the top of the south island (the largest landmass on the map) move half way down the west coast.

Answer 1

This implements the function for snapping provided by @TimSalabim in the answer linked to by @dieghernan (https://stackoverflow.com/a/51300037) as well as the suggestion from @Tim Assal for excluding points on land.

# Code follows from yours immediately after preparation of data. Apologies 
# for bringing in tmap for plotting, but it's what I've been using lately 
# so it was way quicker. 
library(tmap)

# Identify points outside the polygons
outside1 <- sapply(st_intersects(SeaCoor.mp, coastline), function(x){length(x)==0})
seapts1 <- SeaCoor.mp[outside1,]

# Create basemap
bm <- tm_shape(coastline) +
  tm_borders()

# Create map of all points before snap
m1 <- bm +
  tm_shape(SeaCoor.mp) +
  tm_bubbles(size = 0.15, col = 'red') +
  tm_layout(legend.show = FALSE, main.title = 'All points pre snap')

# Create map of points in sea before snap
m2 <- bm +
  tm_shape(seapts1) +
  tm_bubbles(size = 0.15, col = 'red') +
  tm_layout(legend.show = FALSE, main.title = 'Points in sea pre snap')

# Function from @TimSalabim copy+pasted directly from the link provided 
# by @dieghernan: https://stackoverflow.com/a/51300037
st_snap_points = function(x, y, max_dist = 1000) {
  
  if (inherits(x, "sf")) n = nrow(x)
  if (inherits(x, "sfc")) n = length(x)
  
  out = do.call(c,
                lapply(seq(n), function(i) {
#Consider adding an st_buffer() around y in the next line with a super small negative value in your relevant units if you are having floating point math issues--see comments!
                  nrst = st_nearest_points(st_geometry(x)[i], y)
                  nrst_len = st_length(nrst)
                  nrst_mn = which.min(nrst_len)
                  if (as.vector(nrst_len[nrst_mn]) > max_dist) return(st_geometry(x)[i])
                  return(st_cast(nrst[nrst_mn], "POINT")[2])
                })
  )
  return(out)
}

# Perform snap setting maximum distance to 70 m
SeaCoor.mp[outside1,] <- st_snap_points(SeaCoor.mp[outside1,], coastline, 70) %>% 
  st_cast('MULTIPOINT') %>% 
  st_as_sf()

# Identify points that are still in the sea after snapping
outside2 <- sapply(st_intersects(SeaCoor.mp, coastline), function(x){length(x)==0})
seapts2 <- SeaCoor.mp[outside2,]

# Map of all points after snap
m3 <- bm +
  tm_shape(SeaCoor.mp) +
  tm_bubbles(size = 0.15, col = 'red') +
  tm_layout(legend.show = FALSE, main.title = 'All points after snap')

# Map of all points in sea after snap
m4 <-bm +
  tm_shape(seapts2) +
  tm_bubbles(size = 0.15, col = 'red') +
  tm_layout(legend.show = FALSE, main.title = 'Points in sea after snap')

# Call to plot
tmap_arrange(m1, m2, m3, m4)

As you can see the two maps with all points are virtually indistinguishable, while the two maps with points in the sea show that a few of them have nonetheless been snapped to the polygons, which is to be expected with the relatively short maximum snapping distance of 70 m.

Answer 2

Regarding question 1, you could select points that are located in the sea:

outside <- sapply(st_intersects(SeaCoor.mp, coastline),function(x){length(x)==0})

This would give you a logical (TRUE/FALSE) vector that you can subset with:

sea_points <- SeaCoor.mp[outside, ]
sea_points

Then you could snap sea_points to the coastline, then join back to the land points.

Answer 3

For Q2, see issue https://github.com/r-spatial/sf/issues/792 and potential workaround on https://stackoverflow.com/a/51300037.

I would suggest to st_segmentize(st_cast(YOURSHAPE,"LINESTRING")) in order to add more points to the segments and try to snap to that segmentized shape.