超大影像柵格轉向量快速實現

我愛木葉123qq發表於2020-11-26

      距離上一次部落格更新，起碼又是大半年，時光飛逝，我也已經老了。。。這一次，我解決了一個工程上的小問題，可能在行家看來簡單，但是呢，它好像又沒那麼簡單，就是我們通常用的柵格轉向量，

我們知道柵格轉向量，通常有以下方法：採用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.”

 

相關文章