Python 3.6中基于GeoDataFrame点数据创建10m分辨率栅格（禁用make_geocube）

问题描述 #

现有包含geometry（点要素）和Value列的geopandas.geodataframe.GeoDataFrame对象，需生成10m分辨率的TIFF栅格文件，栅格对应点位的取值为Value列的值。因库版本限制无法使用make_geocube，寻求替代方案。

示例数据集：

解决方案：基于Rasterio + GeoPandas实现 #

可以通过rasterio手动构建栅格元数据、赋值点位值后写入TIFF文件，步骤如下：

1. 导入依赖库 #

import geopandas as gpd
import rasterio
from rasterio.transform import from_origin
import numpy as np

2. 准备GeoDataFrame数据 #

假设已加载好GeoDataFrame（若未加载，可通过gpd.read_file()读取或手动创建），需确保数据的CRS正确（示例数据为UTM坐标系，需对应设置）：

# 示例：手动创建GeoDataFrame（实际可替换为你的数据加载逻辑）
data = {
    'geometry': [
        'POINT (520595.000 5720335.000)', 'POINT (520605.000 5720335.000)',
        'POINT (520615.000 5720335.000)', 'POINT (520625.000 5720335.000)',
        'POINT (520635.000 5720335.000)', 'POINT (520645.000 5720335.000)',
        'POINT (520655.000 5720335.000)', 'POINT (520665.000 5720335.000)',
        'POINT (520675.000 5720335.000)', 'POINT (520685.000 5720335.000)',
        'POINT (520695.000 5720335.000)', 'POINT (520705.000 5720335.000)',
        'POINT (520715.000 5720335.000)', 'POINT (520725.000 5720335.000)',
        'POINT (520735.000 5720335.000)'
    ],
    'Value': [
        536.678345, 637.052185, 1230.553955, 944.970642, 1094.613281,
        1123.185181, 849.37634, 1333.459839, 492.866608, 960.957214,
        539.401978, 573.015625, 970.386536, 390.315094, 642.036865
    ]
}
gdf = gpd.GeoDataFrame(data, geometry=gpd.GeoSeries.from_wkt(data['geometry']))
# 设置CRS（示例为UTM 32N，需根据你的数据实际坐标系调整）
gdf.crs = 'EPSG:32632'

3. 定义栅格参数 #

• 分辨率：10m
• 计算栅格的边界范围（基于点数据的extent）
• 计算栅格的行列数

# 分辨率（单位：米，与CRS单位一致）
resolution = 10

# 获取点数据的边界
xmin, ymin, xmax, ymax = gdf.total_bounds

# 调整边界，使得每个点对应栅格中心
xmin_adjusted = xmin - resolution/2
xmax_adjusted = xmax + resolution/2
ymin_adjusted = ymin - resolution/2
ymax_adjusted = ymax + resolution/2

width = int((xmax_adjusted - xmin_adjusted) / resolution)
height = int((ymax_adjusted - ymin_adjusted) / resolution)

4. 创建栅格变换（Transform） #

# 构建栅格变换（左上角坐标 + 分辨率）
transform = from_origin(xmin_adjusted, ymax_adjusted, resolution, resolution)

5. 初始化栅格数组并赋值 #

# 创建空栅格数组，用np.nan填充空值
raster_array = np.full((height, width), np.nan, dtype=np.float32)

# 遍历每个点，计算其在栅格中的行列位置并赋值
for idx, row in gdf.iterrows():
    # 获取点坐标
    x, y = row.geometry.x, row.geometry.y
    # 计算对应的列和行（rasterio的行列是逆序的，行对应y轴，列对应x轴）
    col = int((x - xmin_adjusted) / resolution)
    row_idx = int((ymax_adjusted - y) / resolution)
    # 赋值
    raster_array[row_idx, col] = row['Value']

6. 写入TIFF文件 #

# 设置TIFF文件的元数据
metadata = {
    'driver': 'GTiff',
    'height': height,
    'width': width,
    'count': 1,
    'dtype': 'float32',
    'crs': gdf.crs,
    'transform': transform,
    'nodata': np.nan
}

# 写入文件
output_path = 'output_raster.tif'
with rasterio.open(output_path, 'w', **metadata) as dst:
    dst.write(raster_array, 1)

关键说明 #

• 确保GeoDataFrame的CRS与栅格输出的CRS一致，否则会出现坐标偏移问题。
• 若点数据的坐标是栅格的左上角而非中心，需调整xmin_adjusted和ymax_adjusted的计算逻辑，直接使用原始的xmin和ymax即可。
• 若存在多个点落在同一个栅格单元格内，可根据需求修改赋值逻辑（如取均值、最大值等）。

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