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

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

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

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

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

 1 def RasterToPoly(rasterName, shpName):
 2     """
 3         栅格转矢量
 4         :param rasterName: 输入分类后的栅格名称
 5         :param shpName: 输出矢量名称
 6         :return:
 7    """
 8     inraster = gdal.Open(rasterName)  # 读取路径中的栅格数据
 9     inband = inraster.GetRasterBand(1)  # 这个波段就是最后想要转为矢量的波段,如果是单波段数据的话那就都是1
10     prj = osr.SpatialReference()
11     prj.ImportFromWkt(inraster.GetProjection())  # 读取栅格数据的投影信息,用来为后面生成的矢量做准备
12 
13     outshp = shpName
14     drv = ogr.GetDriverByName("ESRI Shapefile")
15     if os.path.exists(outshp):  # 若文件已经存在,则删除它继续重新做一遍
16         drv.DeleteDataSource(outshp)
17     Polygon = drv.CreateDataSource(outshp)  # 创建一个目标文件
18     Poly_layer = Polygon.CreateLayer(shpName[:-4], srs=prj, geom_type=ogr.wkbMultiPolygon)  # 对shp文件创建一个图层,定义为多个面类
19     newField = ogr.FieldDefn('Value', ogr.OFTReal)  # 给目标shp文件添加一个字段,用来存储原始栅格的pixel value
20     Poly_layer.CreateField(newField)
21 
22     gdal.FPolygonize(inband, None, Poly_layer, 0)  # 核心函数,执行的就是栅格转矢量操作
23     Polygon.SyncToDisk()
24     Polygon = None
25 
26     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.”