1

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) 

Before snapping points

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) 

enter image description here

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) 

enter image description here

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.

3
  • Might be easier to start by fixing your projection errors, because those look like big errors that would undermine the validity of further analysis. – Mox Apr 9 '20 at 22:54
  • I appreciate the suggestion but a) the errors stem from the underlying data and I cannot correct for it, and b) at the resolution/extent we are modelling over the errors are negligible, except when they prevent me from acquiring a point estimate of the environmental conditions. – André.B Apr 9 '20 at 23:35
  • Also, to clarify, I am planning on dropping the points that are no where near land. It is just the ones close to the coast I care about. My example included a high tolerance just to demonstrate the issue. – André.B Apr 10 '20 at 0:37
1
+50

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

enter image description here

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.

1
  • I am not sure what I was going wrong, but I had tried the code from @TimSalabim post and I couldn't get it to work! Thank you for the help and putting it all together! – André.B May 6 '20 at 21:51
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.

0
2

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.