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

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

核心函數(shù)

GDAL中柵格轉(zhuǎn)矢量的函數(shù)主要是以下兩個，二者的參數(shù)沒有任何區(qū)別，只是功能有區(qū)別：

FPolygonize(*args, **kwargs)

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

將每個像元轉(zhuǎn)成一個矩形。

Polygonize(*args, **kwargs) **

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

將每個像元轉(zhuǎn)成一個矩形，然后將相似的像元進(jìn)行合并。

轉(zhuǎn)換代碼

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)  # 讀取路徑中的柵格數(shù)據(jù)
 inband = inraster.GetRasterBand(1)  # 這個波段就是最后想要轉(zhuǎn)為矢量的波段，如果是單波段數(shù)據(jù)的話那就都是1
 prj = osr.SpatialReference()
 prj.ImportFromWkt(inraster.GetProjection())  # 讀取柵格數(shù)據(jù)的投影信息，用來為后面生成的矢量做準(zhǔn)備
 outshp = path[:-4] + ".shp"  # 給后面生成的矢量準(zhǔn)備一個輸出文件名，這里就是把原柵格的文件名后綴名改成shp了
 drv = ogr.GetDriverByName("ESRI Shapefile")
 if os.path.exists(outshp):  # 若文件已經(jīng)存在，則刪除它繼續(xù)重新做一遍
  drv.DeleteDataSource(outshp)
 Polygon = drv.CreateDataSource(outshp)  # 創(chuàng)建一個目標(biāo)文件
 Poly_layer = Polygon.CreateLayer(path[:-4], srs=prj, geom_type=ogr.wkbMultiPolygon)  # 對shp文件創(chuàng)建一個圖層，定義為多個面類
 newField = ogr.FieldDefn('value', ogr.OFTReal)  # 給目標(biāo)shp文件添加一個字段，用來存儲原始柵格的pixel value,浮點(diǎn)型，
 Poly_layer.CreateField(newField)
 gdal.Polygonize(inband, None, Poly_layer, 0)  # 核心函數(shù)，執(zhí)行的就是柵格轉(zhuǎn)矢量操作
 # gdal.FPolygonize(inband, None, Poly_layer, 0)  # 只轉(zhuǎn)矩形，不合并
 Polygon.SyncToDisk()
 Polygon = None
 end_time = datetime.datetime.now()
 print("Succeeded at", end_time)
 print("Elapsed Time:", end_time - start_time)  # 輸出程序運(yùn)行所需時間

轉(zhuǎn)換效果

• 使用FPolygonize
  

轉(zhuǎn)換之后的矢量數(shù)據(jù)有270MB，非常大，打開非?？?/p>

• 使用Polygonize
  

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

到此這篇關(guān)于python 使用GDAL實(shí)現(xiàn)柵格tif轉(zhuǎn)矢量shp的文章就介紹到這了,更多相關(guān)python柵格tif轉(zhuǎn)矢量shp內(nèi)容請搜索本站以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持本站！

香港穩(wěn)定服務(wù)器

版權(quán)聲明：本站文章來源標(biāo)注為YINGSOO的內(nèi)容版權(quán)均為本站所有，歡迎引用、轉(zhuǎn)載，請保持原文完整并注明來源及原文鏈接。禁止復(fù)制或仿造本網(wǎng)站，禁止在非maisonbaluchon.cn所屬的服務(wù)器上建立鏡像，否則將依法追究法律責(zé)任。本站部分內(nèi)容來源于網(wǎng)友推薦、互聯(lián)網(wǎng)收集整理而來，僅供學(xué)習(xí)參考，不代表本站立場，如有內(nèi)容涉嫌侵權(quán)，請聯(lián)系alex-e#qq.com處理。