如何用R/Python在大区域创建等面积空间网格？

我在创建覆盖非洲大陆的12km空间网格时遇到精度难题：多数工具基于投影坐标系（CRS），在大区域场景下精度不足，我需要能直接基于地理CRS（如WGS84）生成网格的工具。以下是用R语言演示的问题：

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.2.3, PROJ 7.2.1; sf_use_s2() is TRUE
library(magrittr)

# 非洲的边界框
africa_bbox <- rbind(c(-26, 55), c(-36, 38))
dimnames(africa_bbox) <- list(c("lon", "lat"), c("min", "max"))
africa_bbox %<>% t()
print(africa_bbox)
#>     lon lat
#> min -26 -36
#> max  55  38

# 创建几何对象
africa_sfc <- africa_bbox %>% 
  as.data.frame() %>% 
  st_as_sf(coords = c("lon", "lat"), crs = "EPSG:4326") %>% 
  st_bbox() %>% st_as_sfc() 
print(africa_sfc)
#> Geometry set for 1 feature 
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -26 ymin: -36 xmax: 55 ymax: 38
#> Geodetic CRS:  WGS 84
#> POLYGON ((-26 -36, 55 -36, 55 38, -26 38, -26 -...

st_area(africa_sfc)   # 区域面积
#> 7.706798e+13 [m^2]
# 地理CRS下该方法无法正常工作
st_make_grid(africa_sfc, cellsize = c(12000, 12000)) 
#> Geometry set for 1 feature 
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -26 ymin: -36 xmax: 11974 ymax: 11964
#> Geodetic CRS:  WGS 84
#> POLYGON ((-26 -36, 11974 -36, 11974 11964, -26 ...

# 要正常生成网格，需要转换到米制投影坐标系
# 使用位于非洲中部的UTM 34N坐标系
africa_sfc_metric <- africa_sfc %>% st_transform("EPSG:32634") 
print(africa_sfc_metric)
#> Geometry set for 1 feature 
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -3842510 ymin: -5189967 xmax: 3613422 ymax: 5419593
#> Projected CRS: WGS 84 / UTM zone 34N
#> POLYGON ((-3842510 -5189967, 3613422 -4567059, ...

# 生成12km网格
africa_12km <- st_make_grid(africa_sfc_metric, cellsize = c(12000, 12000)) 
head(africa_12km, 3)
#> Geometry set for 3 features 
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -3842510 ymin: -5189967 xmax: -3806510 ymax: -5177967
#> Projected CRS: WGS 84 / UTM zone 34N
#> POLYGON ((-3842510 -5189967, -3830510 -5189967,...
#> POLYGON ((-3830510 -5189967, -3818510 -5189967,...
#> POLYGON ((-3818510 -5189967, -3806510 -5189967,...

length(africa_12km) # 网格数量
#> [1] 550470
areas = st_area(africa_12km)
all(unclass(signif(areas, 4)) == 12000^2) # 检查单个网格面积
#> [1] TRUE
sum(areas) / st_area(africa_sfc) # 网格总面积比实际区域大2.85%
#> 1.028542 [1]

# 对比计算：非洲大陆外围12km缓冲区的面积占比
perimeter_12km_area <- africa_sfc %>% st_cast("MULTILINESTRING") %>% st_length() %>% multiply_by(12000)
# 仅占实际面积的0.5%，可见2.85%的误差已十分显著
perimeter_12km_area / st_area(africa_sfc)
#> 0.004725914 [1/m]

^{创建于2022-12-02，使用reprex包(v2.0.1)}

我可以基于Haversine公式自行开发工具，但想先了解是否已有现成的R或Python开源解决方案。

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