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

我愛木葉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.”

 

相關文章