如何使用 sf 库查找多面体的质心

Question

我有一个多面体的 shapefile，我想找到每个多面体的质心，但我得到的是整个 shapefile 的质心。

我使用以下代码行将我的 csv 文件转换为 shapefile（请参阅此问题 Converting CSV file to shapefile - but want polygon not points）：

df <- as.data.frame(read_csv("/Users/alan/Desktop/shapes.csv"))
df1 <- lapply(split(df, df$shape), function(x) { coords <- as.matrix(cbind(x$longitude, 
x$latitude)); rbind(coords, coords[1,])}) 
Coord_Ref <- st_crs(3035)
plot_locations_df <-  st_multipolygon(x=list(df1))
plot_locations_df <-  st_sfc(plot_locations_df, crs = Coord_Ref)
plot(st_geometry(plot_locations_df))
plot(st_geometry(plot_locations_df, col = sf.colors(12, categorical=TRUE), border='grey', 
axes=TRUE))
plot(st_geometry(st_centroid(plot_locations_df)), pch = 3, col = 'red', add=TRUE)

我的数据框如下所示：

structure(list(shape = c(1.1, 1.1, 1.1, 1.1, 2.1, 2.1, 2.1, 2.1, 
3.1, 3.1, 3.1, 3.1, 4.1, 4.1, 4.1, 4.1, 5.1, 5.1, 5.1, 5.1, 6.1, 
6.1, 6.1, 6.1, 7.1, 7.1, 7.1, 7.1, 8.1, 8.1, 8.1, 8.1, 9.1, 9.1, 
9.1, 9.1), longitude = c(43, 43, 40, 40, 23, 23, 20, 20, 25, 
25, 38, 38, 25, 25, 38, 38, 45, 50, 50, 45, 65, 60, 60, 65, 60, 
60, 80, 80, 60, 60, 80, 80, 20, 20, 80, 80), latitude = c(10, 
13, 13, 10, 10, 13, 13, 10, 10, 10.3, 10.3, 10, 12.7, 13, 13, 
12.7, 13, 13, 10, 10, 13, 13, 10, 10, 9.8, 9.5, 9.5, 9.8, 6, 
5.7, 5.7, 6, 5, 4.5, 4.5, 5)), row.names = c(NA, 36L), class = "data.frame")

Answer 1

您可以使用sfheaders直接从您的 .csv 创建多边形

library(sfheaders)
library(sf)

sf <- sfheaders::sf_polygon(
  obj = df
  , polygon_id = "shape"
  , x = "longitude"
  , y = "latitude"
)

## Then you can get the centroids of each polygon
sf_centers <- sf::st_centroid( sf )

Answer 2

这里有两个不同的问题。首先，您用于创建sf 对象的方法会导致几何无效（st_is_valid(plot_locations_df) 返回 FALSE）。要获得有效的 MULTIPOLYGON，您可以使用：

df1 <- lapply(split(df, df$shape), function(x) { coords <- as.matrix(cbind(x$longitude, 
x$latitude)); list(rbind(coords, coords[1,]))}) 
names(df1)<- NULL
Coord_Ref <- st_crs(3035)
plot_locations_df <-  st_sfc(st_multipolygon(x=df1), crs = Coord_Ref)

st_is_valid(plot_locations_df)
[1] TRUE

但是，这仍然无济于事，因为您的几何图形仍然是一个 MULTIPOLYGON（即，由多个 POLYGONS 组成的单个要素），并且 MULTIPOLYGON 的质心是一个单点，计算时考虑了它的所有多边形。

plot_locations_df
> Geometry set for 1 feature 
> geometry type:  MULTIPOLYGON
> dimension:      XY
> bbox:           xmin: 20 ymin: 4.5 xmax: 80 ymax: 13
> projected CRS:  ETRS89-extended / LAEA Europe
> MULTIPOLYGON (((43 10, 43 13, 40 13, 40 10, 43 ...

st_centroid(plot_locations_df)
> Geometry set for 1 feature 
> geometry type:  POINT
> dimension:      XY
> bbox:           xmin: 49.10736 ymin: 8.969325 xmax: 49.10736 ymax: 8.969325
> projected CRS:  ETRS89-extended / LAEA Europe
> POINT (49.10736 8.969325)

为了得到你想要的，你必须通过重铸到 POLYGON 来分割多边形：

plot_locations_df <- st_cast(plot_locations_df, "POLYGON")
plot_locations_df

> Geometry set for 9 features 
> geometry type:  POLYGON
> dimension:      XY
> bbox:           xmin: 20 ymin: 4.5 xmax: 80 ymax: 13
> projected CRS:  ETRS89-extended / LAEA Europe
> First 5 geometries:
> POLYGON ((43 10, 43 13, 40 13, 40 10, 43 10))
> POLYGON ((23 10, 23 13, 20 13, 20 10, 23 10))
> POLYGON ((25 10, 25 10.3, 38 10.3, 38 10, 25 10))
> POLYGON ((25 12.7, 25 13, 38 13, 38 12.7, 25 12...
> POLYGON ((45 13, 50 13, 50 10, 45 10, 45 13))

st_centroid(plot_locations_df)
> Geometry set for 9 features 
> geometry type:  POINT
> dimension:      XY
> bbox:           xmin: 21.5 ymin: 4.75 xmax: 70 ymax: 12.85
> projected CRS:  ETRS89-extended / LAEA Europe
> First 5 geometries:
> POINT (41.5 11.5)
> POINT (21.5 11.5)
> POINT (31.5 10.15)
> POINT (31.5 12.85)
> POINT (47.5 11.5)

plot(st_geometry(plot_locations_df))
plot(st_geometry(plot_locations_df, col = sf.colors(12, categorical=TRUE), border='grey', 
axes=TRUE))
plot(st_geometry(st_centroid(plot_locations_df)), pch = 3, col = 'red', add=TRUE)

HTH！