使用Python的make_geocube栅格化Shapefile生成空GeoTIFF的原因

问题现象 #

• 对属性为CA的海冰Shapefile执行栅格化后，生成的GeoTIFF仅1KB，QGIS中打开为空文件
• 触发警告信息：

C:\Users\ccw\anaconda3\envs\geo_env\lib\site-packages\rasterio_init_.py:230: NotGeoreferencedWarning: The given matrix is equal to Affine.identity or its flipped counterpart. GDAL may ignore this matrix and save no geotransform without raising an error. This behavior is somewhat driver-specific.
s = writer(path, mode, driver=driver,

原使用代码 #

import geopandas as gpd
from geocube.api.core import make_geocube

ds = gpd.read_file('cis_SGRDAEA_20110201_pl_a.shp')
ds['CA'] = ds['CA'].astype(float)
grid = make_geocube(vector_data=ds, measurements=['CA'], resolution=(5000,-5000))
grid.CA.rio.to_raster('test.tif')

问题根源 #

1. 地理参考绑定失败：警告提示仿射变换矩阵为单位矩阵，说明栅格未正确关联原始矢量的地理范围，导致生成的文件无有效空间信息，QGIS无法识别内容
2. make_geocube参数缺失：仅指定分辨率，未明确栅格的空间范围和坐标系，工具无法计算正确的地理变换规则

修复方案 #

步骤1：验证矢量地理信息 #

先确认原始Shapefile的坐标系和范围是否有效：

import geopandas as gpd

ds = gpd.read_file('cis_SGRDAEA_20110201_pl_a.shp')
print("矢量坐标系：", ds.crs)
print("矢量范围：", ds.total_bounds)

步骤2：完善栅格化参数 #

调用make_geocube时强制指定矢量的范围和坐标系：

import geopandas as gpd
from geocube.api.core import make_geocube

ds = gpd.read_file('cis_SGRDAEA_20110201_pl_a.shp')
ds['CA'] = ds['CA'].astype(float)

# 提取矢量的地理参数
vector_bounds = ds.total_bounds
vector_crs = ds.crs

# 生成带正确地理参考的栅格
grid = make_geocube(
    vector_data=ds,
    measurements=['CA'],
    resolution=(5000, -5000),
    bounds=vector_bounds,
    crs=vector_crs
)

grid.CA.rio.to_raster('test_fixed.tif')

步骤3：验证修复结果 #

用Rasterio检查生成的栅格：

import rasterio
with rasterio.open('test_fixed.tif') as src:
    print("栅格坐标系：", src.crs)
    print("地理变换矩阵：", src.transform)
    print("栅格尺寸：", src.shape)

额外注意事项 #

• 若原始PRJ文件异常，可手动指定正确的EPSG坐标系代码（例如crs="EPSG:3411"）替换vector_crs
• 分辨率的正负号需匹配坐标系的轴方向，确保栅格行/列方向正确

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