如何获取matplotlib hexbin图中六边形的精确面积以计算质量面密度？

我想用matplotlib的hexbin绘制质量面密度图，输入是包含[经度(lon)、纬度(lat)、质量(mass)]的数千条数据的DataFrame。hexbin速度很快，我调用时设置了reduce_C_function=np.sum来获取每个六边形内的质量总和，现在需要将这个总和除以六边形面积得到单位面积质量，但没法获取生成的不规则六边形的精确面积，目前只能通过代码获取六边形中心位置，请问有什么方法可以获取精确面积？

代码示例 #

# 计算x方向六边形数量，使六边形宽度约为100米
nx = round(111000 * lenght_x / 100)

hb = df.plot.hexbin(
        x="lon",
        y="lat",
        C="mass",
        gridsize=nx,
        cmap="viridis",
        mincnt=3,
        ax=ax,
        reduce_C_function=np.sum,
    )

# 获取六边形中心
pollycollection = hb.get_children()[0]
centers = pollycollection.get_offsets()
x_c = [p[0] for p in centers]
y_c = [p[1] for p in centers]
plt.plot(x_c, y_c, "x-", color="red")

解决方法 #

方法1：通过六边形路径计算实际面积

PolyCollection对象的get_paths()可以获取每个六边形的顶点路径，结合地理坐标转换计算真实面积：

import matplotlib.path as mpath
import numpy as np

pollycollection = hb.get_children()[0]
paths = pollycollection.get_paths()
actual_areas = []

for idx, path in enumerate(paths):
    # 获取当前六边形的顶点经纬度
    vertices = path.vertices
    lon_list, lat_list = vertices[:, 0], vertices[:, 1]
    
    # 转换经纬度为米级坐标（基于中心纬度的近似转换）
    lat_center_rad = np.radians(y_c[idx])
    # 1度经度对应的米数：随纬度变化
    lon_per_m = 111320 * np.cos(lat_center_rad)
    # 1度纬度对应的米数：近似111000米
    lat_per_m = 111000
    
    # 将顶点转换成米坐标
    x_m = (lon_list - lon_list.mean()) * lon_per_m
    y_m = (lat_list - lat_list.mean()) * lat_per_m
    
    # 计算多边形面积
    area = mpath.Path(np.column_stack((x_m, y_m))).area
    actual_areas.append(abs(area))

方法2：利用gridsize计算理论面积（适合小范围数据）

hexbin的六边形是规则排列的，可通过数据范围和gridsize推导面积，再结合纬度转换：

import numpy as np

# 获取经纬度数据范围
x_min, x_max = df['lon'].min(), df['lon'].max()
y_min, y_max = df['lat'].min(), df['lat'].max()

# 计算单个六边形的经纬度宽度
hex_width_deg = (x_max - x_min) / nx
# 六边形的理论经纬度面积公式：正六边形面积 = (3√3/2) * (边长)²
hex_area_deg = (3 * np.sqrt(3) / 2) * (hex_width_deg / 2) ** 2

# 转换为实际地理面积（平方米）
lat_rad = np.radians(y_c)
actual_areas = hex_area_deg * 111000 * 111320 * np.cos(lat_rad)

方法3：先投影坐标再计算（最准确）

将经纬度转换为平面投影坐标系（如UTM），再用hexbin绘制，此时直接计算的面积就是真实地理面积：

from pyproj import Transformer
import matplotlib.path as mpath

# 根据数据所在区域选择对应的UTM投影（示例为UTM 33N）
transformer = Transformer.from_crs("EPSG:4326", "EPSG:32633")
df['x_utm'], df['y_utm'] = transformer.transform(df['lat'], df['lon'])

# 重新计算gridsize（基于100米宽度）
nx = round((df['x_utm'].max() - df['x_utm'].min()) / 100)

hb = df.plot.hexbin(
        x="x_utm",
        y="y_utm",
        C="mass",
        gridsize=nx,
        cmap="viridis",
        mincnt=3,
        ax=ax,
        reduce_C_function=np.sum,
    )

# 直接获取六边形面积（单位：平方米）
pollycollection = hb.get_children()[0]
paths = pollycollection.get_paths()
actual_areas = [abs(mpath.Path(path.vertices).area) for path in paths]

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