python如何使用GDAL實現柵格tif轉矢量shp

這篇文章給大家分享的是有關python如何使用GDAL實現柵格tif轉矢量shp的內容。小編覺得挺實用的，因此分享給大家做個參考，一起跟隨小編過來看看吧。

前言

目前有一張tif格式的柵格影像，需要在web地圖上進行展示，使用動態切片WMS的方式，渲染速度比較慢，而且大的時候會出現模糊的問題。并且后面需要做多期影像的切換，渲染與加載效率也值得關注。

計劃是使用柵格轉矢量的方式，將柵格數據轉為矢量shp文件，然后進行矢量切片，使用Mapbox進行前端動態渲染。在網上查詢了很多資料，有人說使用d3-contour在node.js中生成或者使用rasterio在python中進行轉換，整體過程都比較麻煩，很不易實現。最終選定了使用GDAL進行柵格轉矢量的方法，代碼比較簡單。
原始tif影像（12.8MB）如下：

核心函數

GDAL中柵格轉矢量的函數主要是以下兩個，二者的參數沒有任何區別，只是功能有區別：

FPolygonize(*args, **kwargs)

FPolygonize(Band srcBand, Band maskBand, Layer outLayer, int iPixValField, char options=None, GDALProgressFunc callback=0, void * callback_data=None) -> int

將每個像元轉成一個矩形。

Polygonize(*args, **kwargs) **

Polygonize(Band srcBand, Band maskBand, Layer outLayer, int iPixValField, char ** options=None, GDALProgressFunc callback=0, void * callback_data=None) -> int

將每個像元轉成一個矩形，然后將相似的像元進行合并。

轉換代碼

from osgeo import gdal, ogr, osr
import os
import datetime
import numpy as np

path = "Z_NAFP20210727.tif"


if __name__ == '__main__':
    start_time = datetime.datetime.now()

    inraster = gdal.Open(path)  # 讀取路徑中的柵格數據
    inband = inraster.GetRasterBand(1)  # 這個波段就是最后想要轉為矢量的波段，如果是單波段數據的話那就都是1
    prj = osr.SpatialReference()
    prj.ImportFromWkt(inraster.GetProjection())  # 讀取柵格數據的投影信息，用來為后面生成的矢量做準備

    outshp = path[:-4] + ".shp"  # 給后面生成的矢量準備一個輸出文件名，這里就是把原柵格的文件名后綴名改成shp了
    drv = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(outshp):  # 若文件已經存在，則刪除它繼續重新做一遍
        drv.DeleteDataSource(outshp)
    Polygon = drv.CreateDataSource(outshp)  # 創建一個目標文件
    Poly_layer = Polygon.CreateLayer(path[:-4], srs=prj, geom_type=ogr.wkbMultiPolygon)  # 對shp文件創建一個圖層，定義為多個面類
    newField = ogr.FieldDefn('value', ogr.OFTReal)  # 給目標shp文件添加一個字段，用來存儲原始柵格的pixel value,浮點型，
    Poly_layer.CreateField(newField)

    gdal.Polygonize(inband, None, Poly_layer, 0)  # 核心函數，執行的就是柵格轉矢量操作
    # 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)  # 輸出程序運行所需時間

轉換效果

• 使用FPolygonize
  

轉換之后的矢量數據有270MB，非常大，打開非常卡

• 使用Polygonize
  

合并之后的矢量數據有48MB，相對第一種方法數據量大大減少

感謝各位的閱讀！關于“python如何使用GDAL實現柵格tif轉矢量shp”這篇文章就分享到這里了，希望以上內容可以對大家有一定的幫助，讓大家可以學到更多知識，如果覺得文章不錯，可以把它分享出去讓更多的人看到吧！