for a given point, find closest grid cell with negative value, R

Question

I have a raster type object in R that is a lat/long grid and each cell has a value that is the depth or elevation of that cell (downloaded from marmap's getNOAA.bathy(), bathymap in toy example).

Then I have a list of points with lat,long coordinates (points). For each point, I want to find the closest grid cell to the point that has a negative (depth) value. If there are multiple cells, then I want to pick one of them randomly. Once I have the/one nearest negative grid cell, then I want to find and save a random lat,long that falls within this cell.

The end goal here is to "move" each point to a random location within the nearest cell that has a negative value.

Note that these points nor the bathymap are not currently projected. They are just lats and lons in decimal degrees.

Here is a toy set-up to work with:

library(marmap)
library(ggplot2)

#get a map
bathymap <- getNOAA.bathy(lat1 = -30, lon1 = 65, lat2 = 55, lon2 = 155, resolution = 10, keep = FALSE, antimeridian = FALSE, path = NULL)

#sample points
points = data.frame("x" = c(138.605,139), "y" = c(35.0833,35.6))

#visualize the problem/set-up
ggplot() + 
  
  #build map
  geom_raster(data = bathymap, aes(x=x, y=y, fill=z)) +
  geom_contour(data = bathymap, aes(x=x, y=y, z=z),
               breaks=0, #contour for land
               colour="black", size=0.3) +
  scale_fill_gradient2(low="skyblue", mid="white", high="navy", midpoint = 0) +
  
  #add points
  geom_point(data = points, aes(x=x, y=y), color = "blue", alpha = 1, size = 1) +
 
  lims(x = c(138,139.5),
       y = c(34,36)) +
  coord_fixed() +
  geom_text(data = bathymap, aes(x=x, y=y, fill=z, label=z), size = 2)```

Answer 1

Here is a solution. Basically, for each original geographical point, I did the following:

1. I used marmap::dist2isobath() to locate the closest point with a negative value along a given isobath.
2. Since the location of isobaths can be rather imprecise with low resolution bathymetries, a negative isobath sometimes crosses cells with positive values. So, I checked the depth of each point identified above with marmap::get.depth(). If the depth returned is negative, this is a success, if not, I go back to step 1 and look for the closest point along a deeper isobath (hence the depth_increment argument). This is thus an incremental process that stops only when a point is finally located in a negative cell.
3. I packaged all this in the find_closest_negative() function that you can apply to every point in a dataset with purrr::map2_df().

The function returns a table with:

1. the original coordinates of each point,
2. the coordinates of the closest point in the closest negative cell, and
3. the depth at the new location.

Points that were originally located in negative cells are not affected by the function. Be careful: it can be rather time-consuming with large bathymetries and many points...

library(marmap)
library(tidyverse)

# Custom function to get location of closest point along a given isobath
# bathy: bathymetric object from marmap::getNOAA.bathy
# x, y: longitude and latitude of a single point
# isobath: negative integer. Minimum depth at which the point should be located
# depth_increment: positive integer to look for deeper cells. Big values allow
# the function to run faster but might lead to more imprecise results
find_closest_negative <- function(bathy, x, y, isobath = -1, depth_increment = 50) {
  # Duplicate point to avoid weird dist2isobath() error message
  point <- data.frame(x, y)
  a <- point[c(1,1), ]
  depth_a <- get.depth(bathy, a, locator = FALSE)
  depth <- unique(depth_a$depth)
  if (depth  < 0) {
    return(unique(depth_a))
  } else {
    while(depth >= 0 ) {
      b <- dist2isobath(bathy, a, isobath = isobath)
      depth_b <- get.depth(bathy, b[,4:5], locator = FALSE)
      depth <- unique(depth_b$depth)
      isobath <- isobath - depth_increment
    }
    return(unique(depth_b))
  }
}

# Get (a small) bathymetric data
bathymap <- getNOAA.bathy(lat1 = 34, lon1 = 138, lat2 = 36, lon2 = 140, 
                          resolution = 10)

# Sample points
points <- data.frame("x" = c(138.605, 139, 138.5, 138.85, 138.95, 139.5), 
                     "y" = c(35.0833, 35.6, 34.2, 35.0, 34.8, 34.4))

# For each point, find the closest location within a negative cell
points2 <- points %>% 
  mutate(map2_df(x, y, .f = ~find_closest_negative(bathymap, .x, .y,
                                                   isobath = -1,
                                                   depth_increment = 50)))

points2
#>         x       y      lon      lat depth
#> 1 138.605 35.0833 138.5455 34.99248 -1111
#> 2 139.000 35.6000 139.2727 35.36346  -311
#> 3 138.500 34.2000 138.5000 34.20000 -2717
#> 4 138.850 35.0000 138.7882 35.00619  -101
#> 5 138.950 34.8000 139.0053 34.79250  -501
#> 6 139.500 34.4000 139.5000 34.40000  -181

# Some data wrangling to plot lines between pairs of points below
origin <- points2 %>% 
  mutate(id = row_number()) %>% 
  select(lon = x, lat = y, id)

destination <- points2 %>% 
  mutate(id = row_number()) %>% 
  select(lon, lat, id)

global <- bind_rows(origin, destination)

# Plot the original points and their new location
bathymap %>% 
  ggplot(aes(x, y)) + 
  geom_tile(aes(fill = z)) +
  geom_contour(aes(z = z), breaks = 0, colour = 1, size = 0.3) +
  geom_text(aes(label = z), size = 2) +
  geom_point(data = points2, aes(x = lon, y = lat), color = "red", size = 2) +
  geom_point(data = points2, aes(x = x, y = y), color = "blue") +
  geom_line(data = global, aes(x = lon, y = lat, group = id), size = 0.4) +
  scale_fill_gradient2(low = "skyblue", mid = "white", high = "navy",
                       midpoint = 0) +
  coord_fixed()