R SF包中的多边形内部重心问题

Question

我需要为多边形添加标签，通常使用质心，但是质心不在多边形内部。我发现了这个问题Calculate Centroid WITHIN / INSIDE a SpatialPolygon，但我正在使用sf包。

以下是一个玩具数据

rm(list = ls(all = TRUE)) #start with empty workspace

library(sf)
library(tidyverse)
library(ggrepel)

pol <- st_polygon(list(rbind(c(144, 655),c(115, 666)
                         ,c(97, 660),c(86, 640)
                         ,c(83, 610),c(97, 583)
                         ,c(154, 578),c(140, 560)
                         ,c(72, 566),c(59, 600)
                         ,c(65, 634),c(86, 678)
                         ,c(145, 678),c(144, 655)))) %>%
  st_sfc()

a = data.frame(NAME = "A")
st_geometry(a) = pol

a <- a  %>% 
  mutate(lon = map_dbl(geometry, ~st_centroid(.x)[[1]]),
     lat = map_dbl(geometry, ~st_centroid(.x)[[2]]))

ggplot() +
  geom_sf(data = a, fill = "orange") +
  geom_label_repel(data = a, aes(x = lon, y = lat, label = NAME)) 

由此产生以下结果:

- user3357059

---

你可以用~st_point_on_surface替换~st_centroid。也就是说，如果你不关心多边形的真实质心。 - Mitch

这个问题涉及到了PostGIS中st_PointOnSurface函数的更多信息，请参考https://gis.stackexchange.com/questions/76498/how-is-st-pointonsurface-calculated。 - Mitch

Answer 1

简单的答案是将st_centroid替换为st_point_on_surface。在质心位于多边形内部的情况下，这不会返回真正的质心。

a2 <- a  %>% 
  mutate(lon = map_dbl(geometry, ~st_point_on_surface(.x)[[1]]),
         lat = map_dbl(geometry, ~st_point_on_surface(.x)[[2]]))

ggplot() +
  ggplot2::geom_sf(data = a2, fill = "orange") +
  geom_label_repel(data = a2, aes(x = lon, y = lat, label = NAME))

或者

如果多边形的质心在多边形内部，则使用该质心，否则在多边形内寻找一个点。

st_centroid_within_poly <- function (poly) {

  # check if centroid is in polygon
  centroid <- poly %>% st_centroid() 
  in_poly <- st_within(centroid, poly, sparse = F)[[1]] 

  # if it is, return that centroid
  if (in_poly) return(centroid) 

  # if not, calculate a point on the surface and return that
  centroid_in_poly <- st_point_on_surface(poly) 
  return(centroid_in_poly)
}

a3 <- a  %>% 
  mutate(lon = map_dbl(geometry, ~st_centroid_within_poly(.x)[[1]]),
         lat = map_dbl(geometry, ~st_centroid_within_poly(.x)[[2]]))

ggplot() +
  ggplot2::geom_sf(data = a3, fill = "orange") +
  geom_label_repel(data = a3, aes(x = lon, y = lat, label = NAME)) 

上述函数 st_centroid_within_polygon 是从你所引用的 问题 中针对 sf 包进行改编的。有关 st_point_on_surface 工作方式的更全面的说明可以在此处找到。

- Mitch

Answer 2

在Mitch的回答基础上进行拓展，因为上面提供的st_centroid_within_poly函数只适用于单个多边形。

要在多个多边形上使用，请使用：

ST_Centroid(ST_Union(geom))

这将合并所有多边形并返回它们的中心点。

最初的回答：

st_centroid_within_poly <- function (poly) {

  # check if centroid is in polygon
  ctrd <- st_centroid(poly, of_largest_polygon = TRUE)
  in_poly <- diag(st_within(ctrd, poly, sparse = F))

  # replace geometries that are not within polygon with st_point_on_surface()
  st_geometry(ctrd[!in_poly,]) <- st_geometry(st_point_on_surface(poly[!in_poly,]))

  ctrd
}

- M.M. Middeldorp

---

在st_geometry(ctrd[!in_poly,]) <- st_geometry(st_point_on_surface(poly[!in_poly,]))处失败，出现poly[!in_poly, ] : incorrect number of dimensions错误。 - CoderGuy123