如何高效合并Pandas DataFrame实现三角网格H5数据集栅格化？

假设你的数据结构如下：

• df_CellsAll：每行对应一个三角形，包含cell_id（用于关联数值）、node_id_1/node_id_2/node_id_3三个顶点ID列
• df_NodesAll：包含node_id、X、Y顶点坐标
• df_Var：包含cell_id、value，每行对应一个三角形的数值

放弃循环，用Pandas矢量化操作实现，效率提升显著：

1. 批量合并顶点坐标
   依次将三个顶点ID列与df_NodesAll合并，提取对应坐标：

   # 合并第一个顶点坐标并重命名列
   df_cells = df_CellsAll.merge(
       df_NodesAll.rename(columns={"X": "X1", "Y": "Y1"}),
       left_on="node_id_1", right_on="node_id", how="left"
   )
   # 合并第二个顶点坐标
   df_cells = df_cells.merge(
       df_NodesAll.rename(columns={"X": "X2", "Y": "Y2"}),
       left_on="node_id_2", right_on="node_id", how="left"
   )
   # 合并第三个顶点坐标
   df_cells = df_cells.merge(
       df_NodesAll.rename(columns={"X": "X3", "Y": "Y3"}),
       left_on="node_id_3", right_on="node_id", how="left"
   )
   

   若没有cell_id，可直接用行索引作为关联键，确保与df_Var顺序一致。

2. 矢量化计算质心
   直接通过列运算生成质心坐标，无需遍历每行：

   df_cells["centroid_X"] = (df_cells["X1"] + df_cells["X2"] + df_cells["X3"]) / 3
   df_cells["centroid_Y"] = (df_cells["Y1"] + df_cells["Y2"] + df_cells["Y3"]) / 3
   

3. 关联数值字段
   合并df_Var得到完整的三角形信息表：

   df_Triangles = df_cells.merge(df_Var, on="cell_id", how="left")
   

基于质心的IDW/Kriging插值并非最优选择，利用已有三角网的拓扑关系直接栅格化，效率和精度都更高：

方案1：GeoPandas + Rasterio 直接栅格化 #

1. 转换为GeoDataFrame
   从df_Triangles生成三角形多边形几何：

   import geopandas as gpd
   from shapely.geometry import Polygon
   
   # 批量生成三角形几何（列表推导比apply更高效）
   geometries = [
       Polygon([(x1, y1), (x2, y2), (x3, y3)])
       for x1, y1, x2, y2, x3, y3 in zip(
           df_Triangles["X1"], df_Triangles["Y1"],
           df_Triangles["X2"], df_Triangles["Y2"],
           df_Triangles["X3"], df_Triangles["Y3"]
       )
   ]
   gdf_triangles = gpd.GeoDataFrame(df_Triangles, geometry=geometries, crs="EPSG:XXX") # 替换为你的坐标系
   

2. 生成目标栅格
   定义栅格参数后，直接通过rasterize函数栅格化：

   import rasterio
   import numpy as np
   from rasterio.transform import from_bounds
   
   # 计算数据范围
   xmin, ymin, xmax, ymax = gdf_triangles.total_bounds
   resolution = 100  # 替换为目标分辨率
   width = int((xmax - xmin) / resolution)
   height = int((ymax - ymin) / resolution)
   transform = from_bounds(xmin, ymin, xmax, ymax, width, height)
   
   # 写入栅格文件
   with rasterio.open(
       "output_raster.tif",
       "w",
       driver="GTiff",
       height=height,
       width=width,
       count=1,
       dtype="float32",
       crs=gdf_triangles.crs,
       transform=transform,
   ) as dst:
       rasterio.features.rasterize(
           ((geom, val) for geom, val in zip(gdf_triangles.geometry, gdf_triangles["value"])),
           out=dst.read(1),
           transform=transform,
           fill=np.nan,
           dtype="float32"
       )
   

   此方法直接利用三角网几何赋值，无额外插值损耗，速度快。

方案2：GDAL 原生TIN栅格化 #

处理超大规模数据时，GDAL的C++底层实现效率更高：

from osgeo import gdal
import geopandas as gpd

# 先将GeoDataFrame保存为临时矢量文件（或使用内存数据集）
gdf_triangles.to_file("temp_triangles.shp")

# 调用GDAL Grid工具生成栅格
xmin, ymin, xmax, ymax = gdf_triangles.total_bounds
resolution = 100
width = int((xmax - xmin) / resolution)
height = int((ymax - ymin) / resolution)

gdal.Grid(
    "output_raster.tif",
    "temp_triangles.shp",
    format="GTiff",
    algorithm="linear",  # TIN线性插值
    width=width,
    height=height,
    outputBounds=[xmin, ymin, xmax, ymax],
    zfield="value"  # 数值字段名
)

方案对比 #

• 质心插值会丢失三角网拓扑信息，结果精度低于直接TIN栅格化
• TIN直接栅格化的时间复杂度仅与三角网数量、栅格分辨率相关，整体效率远高于“质心生成+插值”流程

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