如何计算等面积不规则全球网格的每个单元面积？

问题背景 #

我正在使用一套每个单元面积相等的全球网格（经纬度并非等距分布），仅拥有每个网格单元的中心点经纬度数据，需要计算网格单元的面积，但rasterize()仅适用于规则网格，因此需要适配非规则等面积网格的计算方法。

示例数据 #

data<- structure(list(longitude = c(-179.75, -179.75, -179.75, -179.75, 
-179.75, -179.75, -179.75, -179.75, -179.75, -179.75, -179.75, 
-179.75, -179.75, -179.75, -179.75, -179.75, -179.75, -179.75, 
-179.75, -179.75, -179.75, -179.75, -179.75, -179.75, -179.75, 
-179.75, -179.75, -179.75, -179.75, -179.75, -179.75, -179.75, 
-179.75, -179.75, -179.75, -179.75, -179.75, -179.75, -179.75, 
-179.75, -179.75, -179.75, -179.75, -179.75, -179.75, -179.75, 
-179.75, -179.75, -179.75, -179.75), latitude = c(-78.07014, 
-75.83873, -75.83873, -75.83873, -75.83873, -75.83873, -75.83873, 
-75.83873, -75.83873, -75.83873, -75.83873, -75.83873, -75.83873, 
-75.83873, -73.91449, -73.91449, -73.91449, -73.91449, -73.91449, 
-73.91449, -73.91449, -73.91449, -73.91449, -73.91449, -73.91449, 
-72.19456, -72.19456, -72.19456, -72.19456, -72.19456, -72.19456, 
-72.19456, -72.19456, -72.19456, -72.19456, -70.62343, -70.62343, 
-70.62343, -70.62343, -70.62343, -70.62343, -70.62343, -70.62343, 
-70.62343, -69.16692, -69.16692, -69.16692, -69.16692, -69.16692, 
-69.16692), m1 = c(2.23745321326184, 2.24038445612542, 2.3247780533395, 
1.66920796399225, 1.97447588253218, 1.80588800929381, 1.30442582075186, 
1.5896484492615, 1.44699114881535, 0.928007327323556, 0.895906387553201, 
0.333973772256837, 0.219532053834574, 0.183666645885356, 0.436727878024647, 
0.401417546671545, 0.425731067333535, 0.335518595275659, 0.259337721461039, 
0.1902264040069, 0.191442551684625, 0.133596076817284, 0.123893680481917, 
0.133860384950249, 0.120127435493534, 0.126049456057823, 0.106771777795723, 
0.147704945619801, 0.147752100865098, 0.117718960396799, 0.0991150581382824, 
0.0978849741864964, 0.114689160004564, 0.0933802168779779, 0.10092018847029, 
0.105032230450153, 0.0722648215256436, 0.107400133838986, 0.108854539405822, 
0.11324350965921, 0.0977313130700045, 0.103936014165899, 0.0947261373131327, 
0.118145300616623, 0.0876853839421345, 0.0801285545866838, 0.100545927078964, 
0.116151395720277, 0.14843963006203, 0.12462962629185)), row.names = c(NA, 
50L), class = "data.frame")

最优计算方法 #

方法1：通过相邻中心点推导网格边界，计算多边形面积 #

等面积网格的单元边界通常是相邻中心点的垂直平分线。操作步骤：

• 为每个中心点找到其所有相邻的网格中心点（需匹配网格拓扑结构，比如同一经度带的上下点、同一纬度带的左右点）
• 计算每对相邻中心点的球面垂直平分线
• 这些平分线的交点构成当前网格单元的多边形顶点
• 使用geosphere包的areaPolygon()函数计算该球面多边形的面积

示例代码片段：

library(geosphere)

# 假设已获取某网格单元的顶点经纬度矩阵（按顺时针/逆时针顺序排列）
cell_vertices <- matrix(c(
  -180, -79.15,
  -179.5, -79.15,
  -179.5, -77,
  -180, -77
), ncol=2, byrow=TRUE)

# 计算面积，单位为平方米
cell_area <- areaPolygon(cell_vertices)

方法2：利用已知网格类型的内置工具 #

如果你的网格属于标准等面积网格类型（比如HEALPix、ISEA格网），可直接用对应R包的内置函数：

• HEALPix：使用healpy包，通过分辨率参数直接获取单元面积
• ISEA网格：使用dggridR包，dg_area()函数可直接返回每个网格单元的面积

示例（dggridR）：

library(dggridR)

# 创建对应分辨率的ISEA网格
dggs <- dgconstruct(res = 6) # 替换为你的网格实际分辨率

# 获取单个网格单元的面积（单位：平方米）
cell_area <- dg_area(dggs)

方法3：用sf包构造近似网格多边形并计算面积 #

通过中心点生成球面Voronoi多边形，近似网格单元后计算面积：

library(sf)

# 将数据转为sf点对象（WGS84坐标系）
points_sf <- st_as_sf(data, coords = c("longitude", "latitude"), crs = 4326)

# 生成球面Voronoi多边形
voronoi <- st_voronoi(st_union(points_sf), crs = st_crs(4326))
voronoi_sf <- st_sf(st_cast(voronoi, "POLYGON"))

# 转换为等面积投影（如EPSG:6933）后计算面积
voronoi_sf <- st_transform(voronoi_sf, crs = 6933)
voronoi_sf$area_m2 <- st_area(voronoi_sf)

注意：Voronoi多边形是近似方法，若网格有严格拓扑结构，优先使用方法1或方法2。

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