Блог

У меня есть объект растрового типа в R , который представляет собой сетку lat/long, и каждая ячейка имеет значение, равное глубине или высоте этой ячейки (загружено с marmap getNOAA.bathy() , bathymap в примере игрушки).

Затем у меня есть список точек с широкими,длинными координатами ( points ). Для каждой точки я хочу найти ближайшую ячейку сетки к точке, которая имеет отрицательное значение (глубина). Если есть несколько ячеек, то я хочу выбрать одну из них случайным образом. Как только у меня будет/одна ближайшая отрицательная ячейка сетки, я хочу найти и сохранить случайный lat,длинный, который попадает в эту ячейку.

Конечная цель здесь состоит в том, чтобы «переместить» каждую точку в случайное место в ближайшей ячейке, имеющей отрицательное значение.

Обратите внимание, что ни эти точки, ни батимап в настоящее время не проецируются. Это просто латы и лоны в десятичных градусах.

Вот набор игрушек для работы с:

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

Вот решение. В принципе, для каждой исходной географической точки я сделал следующее:

1. Я использовал marmap::dist2isobath() , чтобы найти ближайшую точку с отрицательным значением вдоль заданной изобаты.
2. Поскольку расположение изобат может быть довольно неточным при батиметрии с низким разрешением, отрицательная изобата иногда пересекает ячейки с положительными значениями. Итак, я проверил глубину каждой точки, указанной выше marmap::get.depth() . Если возвращенная глубина отрицательна, это успех, если нет, я возвращаюсь к шагу 1 и ищу ближайшую точку вдоль более глубокой изобаты (отсюда depth_increment и аргумент). Таким образом, это постепенный процесс, который прекращается только тогда, когда точка, наконец, находится в отрицательной ячейке.
3. Я упаковал все это в find_closest_negative() функцию, которую вы можете применить к каждой точке набора purrr::map2_df() данных .

Функция возвращает таблицу с:

1. исходные координаты каждой точки,
2. координаты ближайшей точки в ближайшей отрицательной ячейке и
3. глубина на новом месте.

Точки, которые изначально были расположены в отрицательных ячейках, не подвержены влиянию функции. Будьте осторожны: это может занять довольно много времени при больших батиметриях и множестве точек…

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