如何确定一个点是否位于跨越日期线的 sf 几何图形内？答案

【问题标题】：How to determine if a point lies within an sf geometry that spans the dateline?如何确定一个点是否位于跨越日期线的 sf 几何图形内？
【发布时间】：2021-03-25 05:13:11
【问题描述】：

使用 R 包 sf，我试图确定某些点是否出现在 shapefile 的范围内（在本例中为夏威夷、专属经济区）。有问题的 shapefile 可以在 here 找到。不幸的是，所讨论区域的边界跨越 +/-180 经度，我认为这让我很困惑。（我在 sf 网站上阅读了有关新版本中球面几何的一些业务，但我无法安装该版本。我认为我正在处理的多边形足够“平坦”，可以避免任何这些无论如何问题）。部分问题似乎是我的 shapefile 包含由日期线分解的多个几何图形，但我不确定如何组合它们。

如何使用 sf 判断某些点是否在 shapefile 中某个对象的边界内（恰好跨越日期线）？

我尝试了st_shift_longitude 的各种组合，但无济于事。我也尝试过转换为我认为的平面投影（2163），但没有奏效。

这是我目前正在尝试的方式：

library(sf)
library(maps)
library(ggplot2)
library(tidyverse)

# this is the shapefile from the link above
eez_unshifted <- read_sf("USMaritimeLimitsAndBoundariesSHP/USMaritimeLimitsNBoundaries.shp") %>%
  filter(OBJECTID == 1206) %>%
  st_transform(4326) 

eez_shifted <- read_sf("USMaritimeLimitsAndBoundariesSHP/USMaritimeLimitsNBoundaries.shp") %>%
  filter(OBJECTID == 1206) %>%
  st_transform(4326) %>%
  st_shift_longitude()

# four points, in and out of the geometry, on either side of the dateline
pnts <- tibble(x=c(-171.952474,176.251978,179.006220,-167.922929),y=c(25.561970,17.442716,28.463375,15.991429)) %>%
  st_as_sf(coords=c('x','y'),crs=st_crs(eez_unshifted))  


# these all return false for every point
st_within(pnts,eez_unshifted)
st_within(st_shift_longitude(pnts),eez_unshifted)
st_within(pnts,eez_shifted)
st_within(st_shift_longitude(pnts),eez_shifted)

# these also all return false for every point
st_intersects(pnts,eez_unshifted)
st_intersects(st_shift_longitude(pnts),eez_unshifted)
st_intersects(pnts,eez_shifted)
st_intersects(st_shift_longitude(pnts),eez_shifted)
  
# plot the data just to show that it looks right
wrld2 <- st_as_sf(maps::map('world2', plot=F, fill=T))
ggplot() +
  geom_sf(data=wrld2, fill='gray20',color="lightgrey",size=0.07) +
  geom_sf(data=eez_shifted) + 
  geom_sf(data=st_shift_longitude(pnts)) + 
  coord_sf(xlim=c(100,290), ylim=c(-60,60)) +
  xlab("Longitude") +
  ylab("Latitude")

【问题讨论】：

看起来如果我围绕我的专属经济区几何构造一个凸包，st_within 就像我希望的那样工作，表明 shapefile 不是“完整”几何。有什么方法可以“关闭”形状吗？

标签： r mapping tidyverse gis sf

【解决方案1】：

答案是确保您检查的几何图形是多边形：

> eez_poly <- st_polygonize(eez_shifted)
> st_within(pnts,eez_poly)
although coordinates are longitude/latitude, st_within assumes that they are planar
Sparse geometry binary predicate list of length 4, where the predicate was `within'
 1: 1
 2: (empty)
 3: 1
 4: (empty)

【讨论】：