超大影像栅格转矢量快速实现

      距离上一次博客更新，起码又是大半年，时光飞逝，我也已经老了。。。这一次，我解决了一个工程上的小问题，可能在行家看来简单，但是呢，它好像又没那么简单，就是我们通常用的栅格转矢量，

我们知道栅格转矢量，通常有以下方法：采用Arcgis进行栅格转矢量，然后工程化呢，就用arcpy实现，就可以了，或者用qgis，原理也差不多，编程的话，绕不过去的，当然是GDAL，这个大家都是直接

调用GDAL的栅格转矢量API就完事了，对于小数据，或者几百MB这种数据效率也不错，问题不大，但是当你的数据规模达到几个G，或者几百个G之后，在涉密的单机上，你怎么解决，只能重写，我这里

采用GDAL从源码层实现了超大影像的栅格转矢量，相比于直接调用GDAL的函数接口实现，效率要高太多了，我们看一下栅格转矢量的实现代码：

def RasterToPoly(rasterName, shpName):
    """
        栅格转矢量
        :param rasterName: 输入分类后的栅格名称
        :param shpName: 输出矢量名称
        :return:
   """
    inraster = gdal.Open(rasterName)  # 读取路径中的栅格数据
    inband = inraster.GetRasterBand(1)  # 这个波段就是最后想要转为矢量的波段，如果是单波段数据的话那就都是1
    prj = osr.SpatialReference()
    prj.ImportFromWkt(inraster.GetProjection())  # 读取栅格数据的投影信息，用来为后面生成的矢量做准备

    outshp = shpName
    drv = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(outshp):  # 若文件已经存在，则删除它继续重新做一遍
        drv.DeleteDataSource(outshp)
    Polygon = drv.CreateDataSource(outshp)  # 创建一个目标文件
    Poly_layer = Polygon.CreateLayer(shpName[:-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

    deleteBackground(shpName, 0)  # 删除背景　　　　

 然后转成矢量后，有一些Nodata图层需要删掉，调用如下代码解决：

def deleteBackground(shpName, backGroundValue):
    """
    删除背景,一般背景的像素值为0
    """
    driver = ogr.GetDriverByName('ESRI Shapefile')
    pFeatureDataset = driver.Open(shpName, 1)
    pFeaturelayer = pFeatureDataset.GetLayer(0)
    strValue = backGroundValue

    strFilter = "Value = '" + str(strValue) + "'"
    pFeaturelayer.SetAttributeFilter(strFilter)
    pFeatureDef = pFeaturelayer.GetLayerDefn()
    pLayerName = pFeaturelayer.GetName()
    pFieldName = "Value"
    pFieldIndex = pFeatureDef.GetFieldIndex(pFieldName)

    for pFeature in pFeaturelayer:
        pFeatureFID = pFeature.GetFID()
        pFeaturelayer.DeleteFeature(int(pFeatureFID))

    strSQL = "REPACK " + str(pFeaturelayer.GetName())
    pFeatureDataset.ExecuteSQL(strSQL, None, "")
    pFeatureLayer = None
    pFeatureDataset = None

 我们来看一下时间对比情况（影像大小为9G）：

实现平台	时间
GDAL 栅格转矢量API	1小时43分钟
改进后API	3分45秒

 当然，我没有用高性能服务器跑，用的是家里的一台普通台式机，CPU为AMD 3600，（这里我要插一句，AMD YES！！！干死intel）。12个核心全部跑起来，执行效率真的是相当惊人。

好吧，没图我说什么(这是南京地区的高分辨率二值图影像，感谢南京朋友提供，我也忘记他的名字了。。)：

 

            图 1  转换后的矢量

从这个来看，结果还是不错的。如果想技术交流或源码，QQ：1044625113，请备注超大影像处理。

最后，我想说几句心里话，我其实挺讨厌写那种狗屁论文的，啥用没有，就在哪里讲故事，吹牛逼，我之前也干过这种事，包括现在也在干这种事，唉，没办法，其实很多其他东西远远比

发论文重要的多，比如这种超大影像转矢量问题，其实还有很多其他工程问题，有待我们去解决的。比如沈阳的李国春老师做的事情，其实我很欣赏的，一句话来总结我的观点：“talk is cheap， show me your code.”