求助：基于多边形面积比例的sf点采样实现方案

我需要在包含12个不同面积的MULTIPOLYGON要素的sf对象中采样约10000个点，该对象包含"ID"和"area_sqkm"字段。要求按多边形面积比例分配采样数量，即面积大的多边形采样更多点，面积小的采样更少。

我尝试了以下代码：

pol_sample = st_read("./Data/Polygons/Polygons_merged.shp")
pol_sample = st_transform(pol_sample, 4326)
pol_sample = st_sample(pol_sample, 10000)

但该代码仅进行随机采样，未考虑多边形面积因素。

附sf对象输出信息：

Simple feature collection with 12 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -121.1789 ymin: -38.13566 xmax: 153.2769 ymax: 42.17854
Geodetic CRS:  GCS_unknown
First 10 features:
ID  area_sqkm                       geometry
1 2380240.67 MULTIPOLYGON (((77.46847 11...
2 1609475.32 MULTIPOLYGON (((148.3452 -2...
3  609946.15 MULTIPOLYGON (((118.0327 27...
4  895408.10 MULTIPOLYGON (((36.70649 -1...
5  999426.57 MULTIPOLYGON (((35.4854 -1....
6 4961657.01 MULTIPOLYGON (((-64.37801 -...
7 4930984.79 MULTIPOLYGON (((-93.46362 1...
8 1010392.03 MULTIPOLYGON (((-68.31586 1...
9   79048.80 MULTIPOLYGON (((-80.51816 2...
10   47379.37 MULTIPOLYGON (((-69.62635 1...

核心思路是先根据每个多边形的面积占比计算各自需要采样的点数，再对每个多边形单独采样，最后合并结果。

1. 计算面积占比与采样点数 #

先将数据转换为平面投影坐标系（地理坐标系如4326计算面积会有偏差，推荐使用等面积投影如EPSG:6933，或根据数据范围选择对应UTM分区），再基于area_sqkm字段计算每个多边形的采样数量，最后调整总点数确保接近10000：

# 读取数据
pol_sample = st_read("./Data/Polygons/Polygons_merged.shp")
# 转换为等面积投影（EPSG:6933是全球等面积圆柱投影）
pol_sample = st_transform(pol_sample, 6933)

# 计算每个多边形的采样点数
total_points = 10000
pol_sample$sample_n = round(pol_sample$area_sqkm / sum(pol_sample$area_sqkm) * total_points)

# 调整总点数，避免四舍五入导致总和偏差
diff = total_points - sum(pol_sample$sample_n)
if (diff != 0) {
  # 给面积最大的多边形补/减差值
  max_idx = which.max(pol_sample$area_sqkm)
  pol_sample$sample_n[max_idx] = pol_sample$sample_n[max_idx] + diff
}

2. 按分配点数逐个采样并合并 #

使用循环或purrr工具对每个多边形单独采样，再合并所有采样点：

library(purrr)

# 对每个多边形执行采样
sampled_list = map2(pol_sample$geometry, pol_sample$sample_n, function(geom, n) {
  if (n > 0) {
    st_sample(geom, size = n, type = "random")
  } else {
    st_sfc() # 采样数为0时返回空几何集合
  }
})

# 合并所有采样点为单个几何集合
sampled_geometry = do.call(c, sampled_list)
# 转换为sf对象
sampled_points_sf = st_sf(geometry = sampled_geometry)

# 可选：关联原多边形的ID和面积字段
sampled_points_sf = st_join(sampled_points_sf, pol_sample[, c("ID", "area_sqkm")], join = st_within)

3. 可选：转换回地理坐标系 #

如果需要最终结果使用WGS84（EPSG:4326）：

sampled_points_sf = st_transform(sampled_points_sf, 4326)

内容的提问来源于stack exchange，提问作者sdorji