如何利用OpenStreetMap实现GPS点与陆地多边形的相交检测？

问题 #

手上有数千个GPS点位，部分在陆地、部分在海上，想用OpenStreetMap（OSM）的高分辨率陆地数据区分两类点位，但OSM数据存在几何方向倒置/拓扑孔洞问题，导致所有点位都被误判为海上。

复现代码 #

library(data.table)
library(leaflet)
library(sf)

pts  <- st_as_sf(data,
                 coords = c("longitude", "latitude"),
                 crs    = 4326)

sf_use_s2(FALSE)
# 下载land-polygons-split-4326.zip文件
land <- st_read("land_polygons.shp") 

land_poly <- land[st_geometry_type(land) %in% c("POLYGON", "MULTIPOLYGON"), ]
land_poly <- st_cast(land_poly, "MULTIPOLYGON") 
landm <- st_make_valid(land_poly)
st_crs(landm) 
bbox <- st_as_sfc(st_bbox(c(xmin = -2.5, xmax = 0, 
                            ymin = 36, ymax = 38)), crs = 4326)
st_crs(bbox) <- st_crs(landm)
land_bbox <- st_crop(landm, bbox)

on_land <- st_intersects(pts, land_bbox, sparse = FALSE)[, 1]
pts_os <- pts[!on_land, ]

测试数据 #

data <- structure(list(latitude = c(37.74492, 37.72914, 37.72934, 37.72892, 
37.72896, 37.72916, 37.72876, 37.73037, 37.73095, 37.73687, 37.74537, 
37.75453, 37.7285, 37.72851, 37.72719, 37.72886, 37.72716, 37.72805, 
37.7273, 37.72729, 37.72831, 37.72916, 37.7691, 37.72895, 37.72849, 
37.72842, 37.74227, 37.72788, 37.74085, 37.72928, 37.72702, 37.75486, 
37.72932, 37.72941, 37.72931, 37.7349, 37.7293, 37.72928, 37.72933, 
37.73162, 37.72934, 37.73019, 37.72742, 37.73006, 37.72956, 37.729, 
37.7293, 37.72905, 37.7293, 37.72611, 37.72924, 37.72927, 37.72613, 
37.72934, 37.72621, 37.7263, 37.7292, 37.72923, 37.72924, 37.72929, 
37.72931, 37.72929, 37.78725, 37.73048, 37.85666, 37.72901, 37.73468, 
37.77225, 37.73039, 37.72892, 37.72933, 37.72932, 37.72909, 37.72937, 
37.72919, 37.72943, 37.72953, 37.76244, 37.74427, 37.72896, 37.72899, 
37.72889, 37.72923, 37.7377, 37.72998, 37.72184, 37.72933, 37.74585, 
37.75051, 37.71287, 37.72912, 37.72924, 37.72942, 37.72931, 37.72922, 
37.72158, 37.72931, 37.76327, 37.77759, 37.72931), longitude = c(-0.71772, 
-0.7059, -0.70633, -0.70623, -0.70872, -0.70676, -0.70895, -0.70347, 
-0.70258, -0.70034, -0.69904, -0.70021, -0.70486, -0.70495, -0.70606, 
-0.70636, -0.70615, -0.70418, -0.70588, -0.70322, -0.70501, -0.70643, 
-0.71751, -0.70493, -0.7049, -0.70471, -0.71993, -0.70301, -0.70871, 
-0.70643, -0.70493, -0.72041, -0.7063, -0.70635, -0.70619, -0.70494, 
-0.70627, -0.70627, -0.70626, -0.70992, -0.70627, -0.70703, -0.7049, 
-0.70925, -0.70639, -0.70672, -0.70638, -0.70725, -0.70666, -0.70652, 
-0.70681, -0.70675, -0.70661, -0.70673, -0.70669, -0.70677, -0.70648, 
-0.70656, -0.70648, -0.70655, -0.70674, -0.7066, -0.71919, -0.70615, 
-0.73799, -0.70658, -0.71353, -0.71072, -0.70673, -0.70643, -0.70628, 
-0.70622, -0.70639, -0.70627, -0.70635, -0.70643, -0.70646, -0.72753, 
-0.71354, -0.7062, -0.7063, -0.70623, -0.70643, -0.71149, -0.70928, 
-0.72652, -0.70652, -0.71616, -0.71856, -0.73371, -0.70643, -0.70613, 
-0.70634, -0.70617, -0.70624, -0.71741, -0.70626, -0.73031, -0.73075, 
-0.70641)), row.names = c(NA, -100L), class = c("data.table", 
"data.frame"))

问题根源 #

OSM陆地多边形的拓扑错误本质是多边形环的方向不符合GIS规范：默认外多边形应为顺时针方向（平面坐标系），孔洞为逆时针方向。如果方向倒置，会导致GIS将整个地球识别为“孔洞”，陆地反而成为外围的负空间，最终所有点都被判定在“海上”。

修复方案 #

方案1：修正多边形环方向 #

转换到平面坐标系（如Web墨卡托3857），检查并反转不符合方向的多边形环，再转回WGS84：

library(data.table)
library(sf)

# 加载点位数据
pts  <- st_as_sf(data, coords = c("longitude", "latitude"), crs = 4326)

# 加载并处理OSM陆地数据
sf_use_s2(TRUE) # 启用S2球面几何，提升WGS84下的拓扑处理能力
land <- st_read("land_polygons.shp") 

# 过滤有效多边形并修复基础拓扑
land_poly <- land[st_geometry_type(land) %in% c("POLYGON", "MULTIPOLYGON"), ]
land_poly <- st_cast(land_poly, "MULTIPOLYGON") 
landm <- st_make_valid(land_poly)

# 裁剪到目标区域
bbox <- st_bbox(c(xmin = -2.5, xmax = 0, ymin = 36, ymax = 38), crs = 4326)
land_bbox <- st_crop(landm, bbox)

# 修正多边形环方向
land_bbox <- st_transform(land_bbox, 3857) # 转平面坐标系
land_bbox_fixed <- st_geometry(land_bbox) %>%
  lapply(function(multipoly) {
    lapply(multipoly, function(poly) {
      # 检查外环方向：顺时针为正，逆时针则反转
      if (st_is_closed(poly) && st_area(poly) > 0) {
        if (st_orientation(poly) == 1) {
          poly <- st_reverse(poly)
        }
      }
      poly
    }) %>%
      st_multipolygon()
  }) %>%
  st_sfc(crs = 3857) %>%
  st_transform(4326) %>%
  st_sf()

# 用st_within判断点是否在陆地内部（比st_intersects更严格）
on_land <- st_within(pts, land_bbox_fixed, sparse = FALSE)[, 1]
pts_os <- pts[!on_land, ]

# 验证结果
table(on_land) # 应同时存在TRUE（陆地）和FALSE（海上）的点位

方案2：使用osmdata包获取预处理的OSM数据 #

osmdata包会自动处理部分OSM拓扑问题，无需手动修复方向：

library(osmdata)
library(sf)

# 获取目标区域的陆地数据
land <- opq(bbox = c(-2.5, 36, 0, 38)) %>%
  add_osm_feature(key = "natural", value = "land") %>%
  osmdata_sf() %>%
  .$osm_multipolygons

# 后续点位判断逻辑同方案1
on_land <- st_within(pts, land, sparse = FALSE)[, 1]

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