使用 sf 创建表示子组边界框的多边形

Question

我正在尝试使用 sf 包和 piped/tidyverse 工作流程根据另一列中定义的组生成边界框。我认为它应该像下面这样工作，但似乎 st_bbox 不尊重群体。

我希望收到三个代表来自 a、b 和 c 的点周围的边界框的多边形记录，但我却收到三个代表所有点的边界框的多边形记录。

library(dplyr)
library(sf)

a <- data.frame(group=rep('a',100), lon=rnorm(100,11,.2), lat=rnorm(100,53,.2))
b <- data.frame(group=rep('b',100), lon=rnorm(100,11.5,.2), lat=rnorm(100,53.5,.2))
c <- data.frame(group=rep('c',100), lon=rnorm(100,12,.2), lat=rnorm(100,54,.2))
dat <- rbind(a,b,c)

pts <- dat %>% st_as_sf(coords=c('lon','lat'),crs=4326) 

pts %>%
  group_by(group) %>%
  summarize(geometry = st_as_sfc(st_bbox(geometry)))

这个returns:

Simple feature collection with 3 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 10.34313 ymin: 52.43993 xmax: 12.54254 ymax: 54.54012
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
# A tibble: 3 x 2
  group                                                                                     geometry
  <fct>                                                                                <POLYGON [°]>
1 a     ((10.34313 52.43993, 12.54254 52.43993, 12.54254 54.54012, 10.34313 54.54012, 10.34313 52...
2 b     ((10.34313 52.43993, 12.54254 52.43993, 12.54254 54.54012, 10.34313 54.54012, 10.34313 52...
3 c     ((10.34313 52.43993, 12.54254 52.43993, 12.54254 54.54012, 10.34313 54.54012, 10.34313 52...

Answer 1

一种选择是通过 tidyr::nest 和 purrr::map 使用嵌套数据框。我还使用了包装函数来简化 map 调用

library(tidyverse)
box_sf <- pts %>% 
  group_by(group) %>%
  nest()

bbox_wrap <- function(x) st_as_sfc(st_bbox(x))
box_sf <- box_sf %>% 
  mutate(bbox = map(data, bbox_wrap))

这将为您提供一个边界框列表作为数据框的一列。如果你想转换回 sf 对象，你可以这样做：

box_sf %>% 
  mutate(geometry = st_as_sfc(do.call(rbind, bbox))) %>% 
  select(-data, -bbox) %>% 
  st_as_sf()

似乎有点迂回，我希望看到使用 group_by 的解决方案，如您最初预期的那样

Answer 2

st_bbox() 函数似乎无法与 group_by() 一起使用，因为它从 sf_points 中提取 bbox 属性，而 sf_points 并未为每个单独的组定义。

一种方法是使用以下方法手动创建边界框：

library(dplyr)
library(sf)
library(ggplot2)
library(tidyr)

# function calculates angle with respect to polygon centroid.
# we need this to order the polygon correctly
calc_angle <- function(lon,lat) {
  cent_lon <- mean(lon)
  cent_lat <- mean(lat)
  ang <- atan2(lat - cent_lat, lon - cent_lon)

  return(ang)
}

bbox <-dat %>%
  group_by(group) %>%
  summarise(xmin = min(lon),ymin = min(lat), xmax=max(lon),  ymax = max(lat)) %>%
  gather(x,lon,c('xmin','xmax')) %>%
  gather(y,lat,c('ymin','ymax')) %>%
  st_as_sf(coords=c('lon','lat'),crs=4326,remove=F) %>%
  group_by(group) %>%
  mutate(angle = calc_angle(lon,lat)) %>%
  arrange(angle) %>%
  summarise(do_union=FALSE) %>%
  st_cast('POLYGON')

本质上，我们通过获取每个组的 xmin、xmax、ymin、ymax 来计算我们自己的 bbox。然后我们收集以仅保留 x 和 y 值，并在转换为 polygon.

之前按顺时针顺序对点进行排序

看起来有点乱，但这是使用 group_by() 解决此问题的一种方法。

pts <- dat %>% 
  st_as_sf(coords=c('lon','lat'), crs=4326, remove=F)

ggplot(pts) + geom_sf(aes(col=group)) +
  geom_sf(data=bbox, aes(col=group), fill=NA)

使用 sf 创建表示子组边界框的多边形

Create polygons representing bounding boxes for subgroups using sf

r

dplyr

sf