如何过滤 OSM 形状内的所有点？

Question

我正在尝试从 CitiBikes 数据集中查找曼哈顿内的所有车站。我可以通过以下方式从 OpenStreetMap 获得曼哈顿的形状

q <- opq(bbox="Manhattan")
q <- add_osm_feature(q, key="name", value="Manhattan")
q <- add_osm_feature(q, key="boundary", value="administrative")
administrative <- osmdata_sf (q)
manhattan = administrative$osm_multipolygons[1]
plot(manhattan, main="Manhattan borough")

我可以通过

将数据转换为简单几何点

cbike <- read.csv(file="201903-citibike-tripdata.csv", stringsAsFactors=F, 
                  sep=",", na.strings=c("NA","NaN", "NULL"))
cbike$start.sf <- cbike %>%
  select(lat=start.station.latitude, lon=start.station.longitude) %>%
  st_as_sf(coords=c("lon", "lat"), crs=st_crs(manhattan))

但我不知道哪些点位于manhattan 内。我假设我想使用 st_intersection 或 st_intersects 来获得一个布尔向量，但我对它的工作原理有点迷茫：任何谓词都显示为空。

> inside = cbike$start.sf[1,1]
> inside
Simple feature collection with 1 feature and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -74.00945 ymin: 40.71107 xmax: -74.00945 ymax: 40.71107
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
                    geometry
1 POINT (-74.00945 40.71107)

这是富尔顿街上的一个点，显然在曼哈顿内部。但是，st_intersects(inside, manhattan) 和 st_within(inside, manhattan) 都是空的。 （st_intersects(inside, inside) 和st_intersects(manhattan, manhattan) 都是TRUE，所以我认为不是点或多面体一般不能相交）

（这是osmdata、sf 和dplyr 包。）

Answer 1

如果您只想选择曼哈顿起点站，可以只从 bbox 的坐标中完成：

xmin <- as.numeric(sub("(.*),.*,.*,.*", "\\1", administrative[[1]]))
xmax <- as.numeric(sub(".*,.*,(.*),.*", "\\1", administrative[[1]]))
ymin <- as.numeric(sub(".*,.*,.*,(.*)", "\\1", administrative[[1]]))
ymax <- as.numeric(sub(".*,(.*),.*,.*", "\\1", administrative[[1]]))

cbike$start_sf <- cbike %>%
  select(lat=`Start Station Latitude`, lon=`Start Station Longitude`) %>% 
subset(., (lat>=xmin & lat<=xmax & lon<=ymin & lon>=ymax)) %>%
  st_as_sf(coords=c("lon", "lat"), crs=st_crs(manhattan))

Answer 2

你检查过曼哈顿几何的有效性吗？查看https://www.r-spatial.org/r/2017/03/19/invalid.html#corrup-or-invalid-geometries！您似乎没有使用 lwgeom 库。安装它（可能来自源代码）、加载，然后在 manhattan 上运行 st_make_valid 会得到一个有效的几何图形，该几何图形与曼哈顿的点有非空相交：

manhattan <- st_make_valid(manhattan)

cbike$start.sf <- cbike %>%
  select(lat=start.station.latitude, lon=start.station.longitude) %>%
  st_as_sf(coords=c("lon", "lat"), crs=st_crs(manhattan))

cbike$starts.inside.manhattan <- st_intersects(cbike$start.sf, manhattan, sparse=FALSE)[,1]

（st_intersects的返回值是n×1的两两交集矩阵，只需要第一列）

> cbike$starts.inside.manhattan
   [1]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
  [21]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE