使用GDAL按边界框裁剪栅格时报错：无法创建0x0数据集

问题描述 #

使用GDAL裁剪.tif格式栅格文件时触发错误：Attempt to create 0x0 dataset is illegal,sizes must be larger than zero，请求排查代码问题。

补充信息：
输入tif文件的投影信息：

PROJCS["WGS 84 / UTM zone 44N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",81],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32644"]]

使用的代码：

from osgeo import gdal, osr

fn_input = r'path\File.TIF'

fn_extent = 'path/clip_extent.tif'

utm11_srs = osr.SpatialReference()
utm11_srs.SetWellKnownGeogCS("WGS84")
utm11_srs.SetUTM(11,1)

bounds_utm = (78.54760607220511, 11.670880584988538, 78.58386873862227, 11.649485363749173)

gdal.Warp(fn_extent, fn_input, outputBounds=bounds_utm, outputBoundsSRS=utm11_srs)

问题排查与修复 #

1. 坐标系不匹配 #

• 输入tif的投影是UTM 44N（EPSG:32644），但代码里指定的裁剪边界坐标系是UTM 11N，两者区域完全不重叠，GDAL转换边界后无法匹配原栅格范围，导致生成0x0数据集。
• 修复：你提供的边界值（78.xx、11.xx）是WGS84地理坐标系的经纬度，应将outputBoundsSRS设置为WGS84地理坐标系，而非UTM 11N。

2. 边界值顺序错误 #

• GDAL的outputBounds参数要求顺序为**(minX, minY, maxX, maxY)**，但你提供的bounds_utm中minY（11.67）大于maxY（11.64），计算出的高度为负数，最终触发0x0数据集错误。
• 修复：调整边界顺序为正确的(minX, minY, maxX, maxY)，即交换后两个值。

3. 修正后的代码示例 #

from osgeo import gdal, osr

fn_input = r'path\File.TIF'
fn_extent = 'path/clip_extent.tif'

# 设置为WGS84地理坐标系（EPSG:4326）
wgs84_srs = osr.SpatialReference()
wgs84_srs.SetWellKnownGeogCS("WGS84")

# 修正边界顺序为(minX, minY, maxX, maxY)
bounds_wgs84 = (78.54760607220511, 11.649485363749173, 78.58386873862227, 11.670880584988538)

gdal.Warp(fn_extent, fn_input, outputBounds=bounds_wgs84, outputBoundsSRS=wgs84_srs)

额外验证步骤 #

可以先打印原栅格的地理范围，确认裁剪边界在原栅格覆盖区内：

ds = gdal.Open(fn_input)
gt = ds.GetGeoTransform()
minX = gt[0]
maxY = gt[3]
maxX = minX + gt[1] * ds.RasterXSize
minY = maxY + gt[5] * ds.RasterYSize
print(f"原栅格范围：({minX}, {minY}, {maxX}, {maxY})")
ds = None

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