R语言使用marmap水深数据查找给定点最近负值网格并生成随机坐标

完整实现代码 #

# 加载依赖包
library(marmap)
library(ggplot2)
library(geosphere)
library(dplyr)

# ------------- 测试数据（和提供的示例一致）-------------
bathymap <- getNOAA.bathy(lat1 = -30, lon1 = 65, lat2 = 55, lon2 = 155, resolution = 10, keep = FALSE, antimeridian = FALSE, path = NULL)
points = data.frame("x" = c(138.605,139), "y" = c(35.0833,35.6))

# ------------- 核心处理步骤-------------
# 1. 预处理水深栅格：提取所有水深（z<0）网格，计算栅格分辨率
bathy_df <- as.data.frame(bathymap) %>% filter(z < 0)
# 计算x（经度）、y（纬度）方向的栅格步长
res_x <- diff(sort(unique(bathy_df$x)))[1]
res_y <- diff(sort(unique(bathy_df$y)))[1]

# 2. 定义单点位处理函数
move_to_nearest_water <- function(point_x, point_y, bathy_df, res_x, res_y){
  # 计算当前点到所有水深网格中心的球面距离（单位：米）
  dists <- distHaversine(p1 = c(point_x, point_y), p2 = as.matrix(bathy_df[,c("x","y")]))
  # 筛选距离最小的网格集合
  min_dist <- min(dists)
  candidate_grids <- bathy_df[dists == min_dist,]
  # 多个候选网格时随机选一个
  selected_grid <- candidate_grids[sample(nrow(candidate_grids), 1),]
  # 在选中网格范围内生成随机坐标
  random_x <- runif(1, min = selected_grid$x - res_x/2, max = selected_grid$x + res_x/2)
  random_y <- runif(1, min = selected_grid$y - res_y/2, max = selected_grid$y + res_y/2)
  # 返回结果
  return(data.frame(
    original_x = point_x, original_y = point_y,
    new_x = random_x, new_y = random_y,
    selected_grid_z = selected_grid$z,
    distance_m = min_dist
  ))
}

# 3. 批量处理所有点位
result <- do.call(rbind, lapply(1:nrow(points), function(i){
  move_to_nearest_water(points$x[i], points$y[i], bathy_df, res_x, res_y)
}))

# ------------- 结果可视化验证-------------
ggplot() + 
  geom_raster(data = bathymap, aes(x=x, y=y, fill=z)) +
  geom_contour(data = bathymap, aes(x=x, y=y, z=z), breaks=0, colour="black", size=0.3) +
  scale_fill_gradient2(low="skyblue", mid="white", high="navy", midpoint = 0) +
  # 原始点位
  geom_point(data = points, aes(x=x, y=y), color = "blue", alpha = 1, size = 3) +
  # 移动后点位
  geom_point(data = result, aes(x=new_x, y=new_y), color = "red", alpha = 1, size = 3) +
  lims(x = c(138,139.5), y = c(34,36)) +
  coord_fixed()

关键逻辑说明 #

• 距离计算使用geosphere::distHaversine，适配十进制度经纬度的球面距离计算，比平面欧氏距离精度更高
• 筛选到多个距离相同的最近水深网格时，通过sample函数随机抽取，符合需求
• 随机坐标生成严格限定在选中网格的实际覆盖范围内，不会超出网格边界
• 输出结果包含原始坐标、新坐标、对应网格水深、移动距离四个核心字段，可直接导出使用

内容的提问来源于stack exchange，提问作者rt11