python gdal 经纬度转地理坐标_【Python】GDAL批量栅格转矢量

最近手头有个工作需要实现批量的栅格转矢量操作。一开始我是用的ArcMap 10.2，想的是写个循环直接使用软件提供的栅格转矢量函数就可以，但不知道是软件自身的bug还是因为非正版的原因，每个脚本里无法运行多次栅格转矢量函数，运行多次程序就会直接崩溃，或Python Kernel死掉。调试了很久也未解决于是就打算试试GDAL，没想到非常的方便快捷，开源GIS真的香~

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

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

start_time = datetime.datetime.now()
folder = 'E:/test' #这里就是你的批量栅格存储的文件夹。文件夹里最好除了你的目标栅格数据不要有其他文件了。
os.chdir(folder) #设置默认路径
for raster in os.listdir(): #遍历路径中每一个文件，如果存在gdal不能打开的文件类型，则后续代码可能会报错。
    inraster = gdal.Open(raster)  #读取路径中的栅格数据
    inband = inraster.GetRasterBand(1)  #这个波段就是最后想要转为矢量的波段，如果是单波段数据的话那就都是1
    prj = osr.SpatialReference()  
    prj.ImportFromWkt(inraster.GetProjection())   #读取栅格数据的投影信息，用来为后面生成的矢量做准备
    
    outshp = raster[:-4] + ".shp" #给后面生成的矢量准备一个输出文件名，这里就是把原栅格的文件名后缀名改成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() 
    Polygon = None
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'字段。
• 后面的三个参数目前还不清楚有什么作用，只使用默认值就可以正常运行得到结果。