删除polyton之外的空间点

Question

我正在尝试使用下面的代码从预定义多边形外部删除数据点。

但是，当我绘制删除内容与保留内容时，有一些点看起来像是某些点被错误删除（即蓝色多边形内的红点）和一些未删除的点（即冰岛附近多边形外部的灰色点） ）。如何更改此代码以正确删除或保留这些区域？

library(tidyverse)
library(openxlsx)
library(sp)
library(sf)
library(rnaturalearth)

# download data
zooWide <- read.xlsx("https://dassh.ac.uk/downloads/mba_rdm/NWAtl_CPR%20data.xlsx", 
                     cols = c(1:6, 8:20)) %>% rename(YEAR = Year)

# create list of spatial points for samples
sp <- SpatialPoints(dplyr::select(zooWide, Longitude, Latitude),
                    proj4string = CRS("+proj=longlat +datum=WGS84")) %>% st_as_sf()

# set boundaries for retaining samples
poly <- st_polygon(list(as.matrix(data.frame(lon = c(-6, -6, -50, -50, -64, -64, -54, -54, -6),
                                              lat = c(45, 66, 66, 70, 70, 51, 51, 45, 45))))) %>%
  st_sfc(crs = st_crs("+proj=longlat +datum=WGS84"))

# determine if sample locations fall within polygon
samples = apply(st_intersects(sp, poly, sparse = TRUE), 1, any)

# filter all samples to those within polygon
zooRetained <- filter(zooWide, samples == TRUE)
zooRemoved <- filter(zooWide, samples == FALSE)

# basemap
world <- ne_countries(scale = "medium", returnclass = "sf")

# map
ggplot() +
  # samples
  geom_point(data = zooRemoved, aes(x = Longitude, y = Latitude), alpha = 0.1, col = "red") +
  geom_point(data = zooRetained, aes(x = Longitude, y = Latitude), alpha = 0.1) +
  # region
  geom_sf(data = poly, aes(geometry = geometry), col = "navy", fill = "transparent", linewidth = 2) + 
  # basemap
  geom_sf(data = world, colour = NA, fill = "grey75") +
  scale_y_continuous(labels = ~.x, limits = c(43, 72), breaks = seq(45, 70, 5)) +
  scale_x_continuous(labels = ~.x, limits = c(-68, -2), breaks = seq(-70, 0, 10)) +
  labs(x = NULL, y = NULL)

Answer 1

由于

poly

仅由角定义，因此球体上水平边缘的最短路径将不会遵循平行线，至少在这些纬度上不会。虽然通过向多边形边缘添加更多点很容易解决这个问题，但我们可以使用

st_segmentize()

来实现这一点。

library(dplyr, warn.conflicts = FALSE)
library(openxlsx)
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(rnaturalearth)
#> Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`
library(ggplot2)

# by default sf uses spherical geometires for shapes with geographic coordinates,
# if you just define corners, edges will follow the spehere in a way you probably 
# did not expect; to keep vertiacl edges aligned with parallels, we could add more points
# to poly with st_segmentize()
poly <- 
  st_polygon(list(cbind(lon = c(-6, -6, -50, -50, -64, -64, -54, -54, -6),
                        lat = c(45, 66, 66, 70, 70, 51, 51, 45, 45)))) %>%
  st_segmentize(dfMaxLength = .1) %>% 
  st_sfc(crs = "WGS84")

zooWide <- read.xlsx("https://dassh.ac.uk/downloads/mba_rdm/NWAtl_CPR%20data.xlsx", 
                     cols = c(1:6, 8:20)) %>% rename(YEAR = Year)

# frame to sf object, add a feature indicating point intersactions
zooWide_sf <- zooWide %>% 
  st_as_sf(coords = c("Longitude", "Latitude"), crs = "WGS84") %>% 
  mutate(retained = st_intersects(geometry, poly, sparse = FALSE)[,1])

# some Natural Earth polygons are bit broken, one workaround is disabling s2;
# though you might consider some other source for shapes instead, 
# i.e. giscoR::gisco_get_countries()
sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off
world <- ne_countries(scale = "medium", returnclass = "sf")$geometry %>% 
  st_crop(c(xmin=-68, ymin=43, xmax=-2, ymax=72))
#> although coordinates are longitude/latitude, st_intersection assumes that they
#> are planar

ggplot() +
  geom_sf(data = world, colour = NA, fill = "grey75") +
  geom_sf(data = zooWide_sf, aes(colour = retained), alpha = .1, show.legend = FALSE) +
  geom_sf(data = poly, fill = NA, colour = "navy", linewidth = 2) +
  scale_color_manual(values = c("FALSE" = "red", "TRUE" = "black")) +
  coord_sf(expand = FALSE)

^{创建于 2023-11-18，使用 reprex v2.0.2}