如何按组计算地理质心？

Question

我需要计算按类别分组的不同多边形的地理质心。让我举个例子。假设我们有这个 sf 数据集

library("sf")
nc <- st_read(system.file("/shape/nc.shp", package="sf"))

假设我们将不同的多边形分组。

df <- data.frame(id=rep(1:10,rep=10))
nc <- cbind(nc,df)

如果我们绘制它，我们可以看到不同的多边形属于不同的类别。最重要的是，这些多边形并不总是连续的。

plot(nc["id"])

我想计算每个类别的地理质心。我考虑的选项之一是先分解不同的多边形，然后计算质心。但是，由于多边形是不连续的，所以这不起作用。

p <- nc %>% group_by(id) %>%
  summarise(geometry = st_union(geometry)) %>%
  ungroup()

此外，如果我尝试下面的代码，它会给出每个多边形的质心，而不是按组

p <- nc %>% group_by(id) %>%
    summarise(geometry = st_centroid(geometry),.groups = "keep")

知道如何按组获取质心吗？

Answer 1

如果您在摘要中 group_by、st_union，然后找到质心，这似乎可行。

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(tidyverse)

nc <- st_read(system.file("/shape/nc.shp", package="sf"))
#> Reading layer `nc' from data source 
#>   `/home/x/R/x86_64-pc-linux-gnu-library/4.0/sf/shape/nc.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27

df <- data.frame(id=rep(1:10,rep=10))
nc <- cbind(nc,df)


# new sf data.frame, get centroids by group
nc_grouped <- nc %>%
                group_by(id) %>%
                summarise(st_union(geometry)) %>%
                st_centroid() 
#> Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
#> geometries of x

# plot everything
ggplot() + 
  geom_sf(data = nc, aes(fill = id, alpha = .2)) +
  geom_sf(data = nc_grouped, aes(color = id), size = 3) + 
  scale_fill_viridis_c() +
  scale_color_viridis_c()

很难看到，但它们都在那里。

# plotting only id ==1, looks about right.
ggplot() + 
  geom_sf(data = nc[nc$id == 1,]) + 
  geom_sf(data = nc_grouped[nc_grouped$id == 1,])

只检查第1组，看起来是对的。



nc %>% group_by(id) %>%
  summarise(geometry = st_centroid(geometry),.groups = "keep")
#> Simple feature collection with 100 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -84.05986 ymin: 34.07671 xmax: -75.8095 ymax: 36.49111
#> Geodetic CRS:  NAD27
#> # A tibble: 100 × 2
#> # Groups:   id [10]
#>       id             geometry
#>    <int>          <POINT [°]>
#>  1     1  (-81.49823 36.4314)
#>  2     1 (-79.33478 36.39346)
#>  3     1 (-76.61642 36.14886)
#>  4     1 (-77.98691 35.96228)
#>  5     1 (-81.17403 35.91575)
#>  6     1 (-77.37784 35.58906)
#>  7     1 (-81.91783 35.39871)
#>  8     1 (-80.24928 35.31388)
#>  9     1 (-84.05986 35.13111)
#> 10     1 (-77.10388 35.13065)
#> # … with 90 more rows

^{由 reprex package (v2.0.1)}

于 2022-05-17 创建

如何按组计算地理质心？

How to calculate geographic centroids by group?

gis

r

spatial

rmaps

sf