GDAL 计算栅格影像覆盖范围 (不含 0 值区域)

可以使用 GDAL 提供的 Polygonize 函数来计算栅格影像的覆盖范围，不包括 0 值区域。具体步骤如下：

导入 GDAL 库

from osgeo import gdal, ogr

打开栅格影像

raster_ds = gdal.Open('your_raster.tif')

获取栅格影像的波段数和范围

band_count = raster_ds.RasterCount
x_min, x_pixel_size, _, y_max, _, y_pixel_size = raster_ds.GetGeoTransform()
x_max = x_min + (raster_ds.RasterXSize * x_pixel_size)
y_min = y_max - (raster_ds.RasterYSize * abs(y_pixel_size))

创建空间参考

srs = osr.SpatialReference()
srs.ImportFromWkt(raster_ds.GetProjectionRef())

创建矢量数据源

driver = ogr.GetDriverByName('ESRI Shapefile')
out_shp = driver.CreateDataSource('your_output.shp')

创建图层

out_lyr = out_shp.CreateLayer('your_layer', srs=srs, geom_type=ogr.wkbPolygon)

添加字段

field_defn = ogr.FieldDefn('id', ogr.OFTInteger)
out_lyr.CreateField(field_defn)

使用 Polygonize 函数生成多边形

mask_band = raster_ds.GetRasterBand(1).ReadAsArray()
mask = (mask_band != 0).astype(int)
mask = np.flip(mask, axis=0)
mask = np.rot90(mask, k=3)
mask_ds = gdal.GetDriverByName('MEM').Create('', mask.shape[1], mask.shape[0], 1, gdal.GDT_Byte)
mask_ds.GetRasterBand(1).WriteArray(mask)
mask_ds.SetGeoTransform((x_min, x_pixel_size, 0, y_min, 0, y_pixel_size))
gdal.Polygonize(mask_ds.GetRasterBand(1), None, out_lyr, 0, [], callback=None)

关闭数据源

out_shp = None

通过以上步骤，我们就可以得到一个不包括 0 值区域的覆盖范围矢量文件。