【Python】GDAL批量栅格转矢量

最近需要实现批量栅格转矢量操作。一开始用Arcmap操作，全国34各省太多，就写了一个脚本，从文件夹里读取循环读取后缀是“.tif”的文件，使用GDAL读取，非常方便快捷，手点到抽筋的步骤，一会就自动读取转换，并自动分类存放。。。

话不多说，直接上代码吧，只需要修改文件夹路径到自己的目标路径，就可以直接运行下面代码（python3）：

import datetime
import os
from osgeo import gdal, osr, ogr

start_time = datetime.datetime.now()
input_path = r'G:\viirs_sys_snpp2\点密度\光伏板栅格2020'
out_path = r"G:\viirs_sys_snpp2\点密度\guanfuban_shp"
os.chdir(input_path) # 设置默认路径
for province in os.listdir(input_path):
    for raster in os.listdir(province):
        filename, txt = os.path.splitext(raster)
        if txt == ".tif":
            raster_path = input_path + "\\" + province + "\\" + raster
            inraster = gdal.Open(raster_path)
            inband = inraster.GetRasterBand(1)  # 要转为矢量的波段
            prj = osr.SpatialReference()
            prj.ImportFromWkt(inraster.GetProjection())  #读取栅格数据结构投影信息

            outshp = out_path + '//' + raster[:-4] + ".shp"
            drv =ogr.GetDriverByName("ESRI Shapefile")
            if os.path.exists(outshp):  # 若文件已经存在，则删除它继续重新做一遍
                drv.DeleteDataSource(outshp)
            Polygon = drv.CreateDataSource(outshp)  #创建一个目标文件夹
            Poly_layer = Polygon.CreateLayer(raster[:-4], srs = prj, geom_type = ogr.wkbMultiPolygon)  #对shp文件创建一个图层，定义为多个面类
            newField = ogr.FieldDefn('value',ogr.OFTReal) # 给目标shp文件添加一个字段，用来存储原始栅格的pixel value
            Poly_layer.CreateField(newField)

            gdal.FPolygonize(inband, None, Poly_layer, 0) # 核心函数，执行的就是栅格转矢量操作
            Polygon.SyncToDisk()
            print(Polygon)
            Polygon = None
            print(raster)
end_time = datetime.datetime.now()
print ("Succeeded at",end_time)
print ("Elapsed Time:",end_time-start_time)  #输出程序运行所需时间

从上面的代码中可以看到，最核心的代码就是gdal.FPolygonize(）这个函数。这个函数的参数有以下几：

• srcBand: 想要转成矢量的那个波段
• maskBand: 转换的范围。这里因为我的是单波段图像，无需掩膜处理，所以给的None。
• outLayer：这个参数是指目标shp文件的图层名
• iPixValField: 栅格值存储在目标shp文件属性表的第几个字段。0代表存储在第一个字段，即我们在代码中建立的’value’字段。
• 后面的三个参数目前还不清楚有什么作用，只使用默认值就可以正常运行得到结果。