问题

我有德国邮政编码的多边形形状数据。对于每个邮政编码，我想计算从其质心到其边界的最大距离，并在地图上对其中一些进行说明。我找到了一个通过sf包和st_cast()以及st_distance()来计算这个最大值的帖子。我的数据以sf数据框的形式存在。

library(sf)
library(tidyverse)

# 获取德国邮政编码形状多边形
URL <- "https://downloads.suche-postleitzahl.org/v2/public/plz-5stellig.shp.zip"

# 使用GDAL虚拟文件系统从远程URL加载压缩的shapefile
GER_postcode <- paste0("/vsizip//vsicurl/", URL) %>% read_sf()

# 将plz列转换为数字
GER_postcode$plz <- as.numeric(GER_postcode$plz)

# 过滤特定的邮政编码
test <- GER_postcode %>% filter(plz == 15232)

# 计算距离
distances <- test %>% 
  st_cast("POINT") %>% 
  st_distance(st_centroid(test))

# 最大距离：
max_dist <- max(distances)
max_dist

ggplot() +
  geom_sf(data = test, size = 1, shape = 16, alpha = 0.1) + # 形状
  geom_sf(data = st_centroid(test)) + # 质心
  theme_bw()

找到的最大值在哪里（1297.496 [米]），如何在地图上显示连接？

<details>
<summary>英文:</summary>

I have polygon shape data for German postcodes. For each postcode I like to calculate the maximum distance from centroid to its border and illustrate this on a map for some of them. I found a post which calculates this maximum via `sf` package and `st_cast()` and `st_distance()`. My data comes as sf dataframe. 

https://stackoverflow.com/questions/48447680/how-to-compute-the-greatest-distance-between-a-centroid-and-the-edge-of-the-poly

    library(sf)
    library(tidyverse)

    # Get German postcode shape polygons
    URL &lt;- &quot;https://downloads.suche-postleitzahl.org/v2/public/plz-5stellig.shp.zip&quot;
    
    # use GDAL virtual file systems to load zipped shapefile from remote url
    GER_postcode &lt;- paste0(&quot;/vsizip//vsicurl/&quot;, URL) %&gt;%  read_sf()
    
    # convert numeric
    GER_postcode$plz &lt;- as.numeric(GER_postcode$plz)

    # filter a specific postcode
    test &lt;- GER_postcode %&gt;% filter(plz == 15232)
    
    # compute distances 
    distances &lt;- test %&gt;% 
      st_cast(&quot;POINT&quot;) %&gt;% 
      st_distance(st_centroid(test))

    # maximum dist:
    max_dist &lt;- max(distances)
    max_dist

    ggplot() +
      geom_sf(data = test, size = 1, shape = 16, alpha = 0.1) + # shape
      geom_sf(data = st_centroid(test)) + # centroid 
      theme_bw()

Where exactly is the found maximum (1297.496 [m]) and how can I show the connection on the map?

[![enter image description here][1]][1]


  [1]: https://i.stack.imgur.com/INXdl.png

</details>


# 答案1
**得分**: 5

你的代码通过将边界的MULTIPOLYGON转换为组成该多边形的点集来计算最大距离，然后计算到每个点的距离。

因此，我们可以找到其中距离最远的点，创建一个包含该点和质心的sf数据框，并使用st_cast()将它们汇总成一个LINESTRING。

```r
# 创建包含边界点的sf对象
border_points <- test %>%
    st_cast("POINT")

# 计算距离
distances <- border_points |>
    st_distance(st_centroid(test))

max_dist_linestring <- border_points |>
    filter(as.logical(distances == max(distances))) |>
    bind_rows(st_centroid(test)) |>
    summarise() |>
    st_cast("LINESTRING")

ggplot() +
    geom_sf(data = test, size = 1, shape = 16, alpha = 0.1) + # 形状
    geom_sf(data = st_centroid(test)) + # 质心
    geom_sf(data = max_dist_linestring) +
    theme_bw()

关于计算质心和投影的注意事项

你的数据是以经度/纬度格式表示的。st_crs(GER_postcode)返回4326，即WGS84，一个经纬度系统。然而，st_centroid()对经纬度数据不会给出准确的结果。

你应该将数据转换为投影坐标系，即平面坐标系。由于你的数据是德国的，你可以使用DE_ETRS89。你可以这样做：

GER_postcode <- st_transform(GER_postcode, crs = 25831)

如果你选择不同的坐标系，请确保st_is_longlat(GER_postcode)为FALSE。这将获得更准确的最大距离。在你发布的示例中，这会导致大约10米的差距。但是，根据位置的不同，你可能会得到完全不正确的结果（即不是实际的最远距离）。请参阅此处的伦敦投影与地理缓冲绘图以获取更多信息。

# create sf object of border points
border_points <- test %>%
    st_cast("POINT")

# compute distances
distances <- border_points |>
    st_distance(st_centroid(test))

max_dist_linestring <- border_points |>
    filter(as.logical(distances == max(distances))) |>
    bind_rows(st_centroid(test)) |>
    summarise() |>
    st_cast("LINESTRING")

ggplot() +
    geom_sf(data = test, size = 1, shape = 16, alpha = 0.1) + # shape
    geom_sf(data = st_centroid(test)) + # centroid
    geom_sf(data = max_dist_linestring) +
    theme_bw()

GER_postcode <- st_transform(GER_postcode, crs = 25831)