时间戳

首页 > 解决方案 > 对于给定的点，找到最近的负值网格单元，R

r - 对于给定的点，找到最近的负值网格单元，R

问题描述

我在 R 中有一个栅格类型对象，它是一个纬度/经度网格，每个单元格都有一个值，即该单元格的深度或高度（从 marmap 下载getNOAA.bathy()，bathymap在玩具示例中）。

然后我有一个带有经纬度坐标（points）的点列表。对于每个点，我想找到最接近具有负（深度）值的点的网格单元。如果有多个单元格，那么我想随机选择其中一个。一旦我有/一个最近的负网格单元格，然后我想找到并保存一个随机的 lat,long 落在这个单元格内。

这里的最终目标是将每个点“移动”到最近的单元格中具有负值的随机位置。

请注意，当前未投影这些点和深海图。它们只是十进制度的纬度和经度。

这是一个可以使用的玩具设置：

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

标签： rgridgeospatialrasternearest-neighbor

解决方案

这是一个解决方案。基本上，对于每个原始地理点，我做了以下工作：

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

推荐阅读