python脚本GIS中批量栅格转面

在Python中，可以使用GDAL库（Geospatial Data Abstraction Library）来实现批量栅格转面的操作。GDAL提供了一套丰富的函数和工具，能够方便地处理地理空间数据，包括栅格到矢量的转换。

以下是一个简单的示例，展示如何使用Python和GDAL进行批量栅格转面：

from osgeo import gdal, ogr

# 定义源目录和目标输出目录
src_dir = 'path/to/source/rasters'
dst_dir = 'path/to/output/polygons'

# 获取源目录下所有栅格文件
raster_files = [f for f in os.listdir(src_dir) if f.endswith('.tif')]

for raster_file in raster_files:
    # 打开栅格文件
    src_dataset = gdal.Open(os.path.join(src_dir, raster_file))
    
    # 创建输出矢量文件
    driver = ogr.GetDriverByName('ESRI Shapefile')
    dst_filename = os.path.splitext(raster_file)[0] + '.shp'
    dst_dataset = driver.CreateDataSource(os.path.join(dst_dir, dst_filename))
    
    # 获取栅格几何信息，创建字段
    srs = src_dataset.GetProjection()
    geom_type = gdal.GetDriverByName('GTiff').GetMetadataItem('GEOTIFF_GEOMETRY_TYPE', 'IMAGE_STRUCTURE')
    
    # 创建数据层
    dst_layer = dst_dataset.CreateLayer(dst_filename, srs=srs)
    dst_layer.CreateField(ogr.FieldDefn("FID", ogr.OFTInteger))  # 创建FID字段
    
    # 遍历栅格的每个像素
    for i in range(src_dataset.RasterYSize):
        for j in range(src_dataset.RasterXSize):
            # 根据栅格值创建面
            pixel_value = src_dataset.ReadAsArray(j, i, 1, 1)
            geom = ogr.Geometry(geom_type)
            if pixel_value != 0:  # 只处理非零值（假设非零代表有地物）
                geom.AddPoint(j, i)  # 添加坐标
                feat = ogr.Feature(dst_layer.GetLayerDefn())
                feat.SetGeometry(geom)
                feat.SetField('FID', feat.GetFID() + 1)  # 增加FID
                dst_layer.CreateFeature(feat)

    # 保存并关闭
    dst_layer.SyncToDisk()
    src_dataset = None
    dst_dataset = None