如何按 long/lat 值过滤 sfc 多边形?

How to filter sfc polygons by long/lat values?

假设t是:

t <- structure(list(structure(list(structure(c(-89.990791, -89.990772, 
-89.990901, -89.99092, -89.990791, 30.727025, 30.727083, 30.727114, 
30.727057, 30.727025), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", 
"sfg")), structure(list(structure(c(-89.991691, -89.991755, -89.991755, 
-89.991691, -89.991691, 30.716004, 30.716004, 30.715916, 30.715916, 
30.716004), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg"
))), class = c("sfc_POLYGON", "sfc"), precision = 0, bbox = structure(c(xmin = -89.991755, 
ymin = 30.715916, xmax = -89.990772, ymax = 30.727114), class = "bbox"), crs = structure(list(
    epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), n_empty = 0L)
> t
Geometry set for 2 features 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -89.99176 ymin: 30.71592 xmax: -89.99077 ymax: 30.72711
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
POLYGON ((-89.99079 30.72703, -89.99077 30.7270...
POLYGON ((-89.99169 30.716, -89.99175 30.716, -...

如何按 long/lat 边界过滤多边形?假设我想过滤掉 hax lat>30.72 的任何多边形(因此只保留第二个多边形)。有什么特定的函数可以用来过滤多边形吗?

不知道有没有现成的功能,就是长方形的“空间滤镜” 很容易建造。只需定义“角”,从中创建一个 bbox,转换 到多边形并找到包含/重叠的原始多边形 "filter area".

这是一个简单的例子:

library(sf)
polys_sf<-st_read(system.file("shape/nc.shp", package="sf") ) %>% 
  st_transform(crs="+init=epsg:4326")
plot(st_geometry(polys_sf))

定义“空间过滤器”

xmin <- -80
xmax <- -76
ymin <- 34
ymax <- 36

根据过滤器创建一个多边形。 (您可以对某些值使用 "NA",因此如果您只想过滤 "on the left",您可以将 xmax 设置为 NA)

filt_bbox <- sf::st_bbox(c(xmin = ifelse(is.na(xmin), -180, xmin), 
                           ymin = ifelse(is.na(ymin),  -90,  ymin), 
                           xmax = ifelse(is.na(xmax), +180, xmax), 
                           ymax = ifelse(is.na(ymax),  +90, ymax)), 
                         crs = st_crs(4326)) %>% 
  sf::st_as_sfc(.)

现在 "filter" 基于 bbox 多边形的原始数据集:如果您只想保留完全包含在定义区域中的多边形,请使用 st_within

find_data <- sf::st_within(polys_sf, filt_bbox)
#> although coordinates are longitude/latitude, st_within assumes that they are planar
filt_data <- polys_sf[which(lengths(find_data) != 0), ]

plot(filt_bbox)
plot(st_geometry(filt_data), add = TRUE, reset = FALSE)

如果要保留与定义区域相交的所有多边形,请使用 st_intersects

find_data <- sf::st_intersects(polys_sf, filt_bbox)
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
filt_data <- polys_sf[which(lengths(find_data) != 0), ]

plot(st_geometry(filt_data))
plot(filt_bbox, add = TRUE)

(很明显,如果你的多边形和“过滤范围”都是lat/long,这就有效,否则 你必须注意重新投影等。) reprex package (v0.2.1)

创建于 2018-11-15

您可以使用一个简单的函数更轻松地完成此操作:

filter_sf <- function(.data, xmin = NULL, xmax = NULL, ymin = NULL, ymax = NULL) {
  bb <- sf::st_bbox(.data)
  if (!is.null(xmin)) bb["xmin"] <- xmin
  if (!is.null(xmax)) bb["xmax"] <- xmax
  if (!is.null(ymin)) bb["ymin"] <- ymin
  if (!is.null(ymax)) bb["ymax"] <- ymax
  sf::st_filter(.data, sf::st_as_sfc(bb), .predicate = sf::st_within)
}

现在你可以做:

t %>%
  filter_sf(ymin = 30.72)

此方法适用于 POLYGONMULTIPOLYGON 几何图形。