前言
目前有一張tif格式的柵格影像,需要在web地圖上進行展示,使用動態(tài)切片WMS的方式,渲染速度比較慢,而且大的時候會出現(xiàn)模糊的問題。并且后面需要做多期影像的切換,渲染與加載效率也值得關注。
計劃是使用柵格轉矢量的方式,將柵格數(shù)據(jù)轉為矢量shp文件,然后進行矢量切片,使用Mapbox進行前端動態(tài)渲染。在網(wǎng)上查詢了很多資料,有人說使用d3-contour
在node.js中生成或者使用rasterio
在python中進行轉換,整體過程都比較麻煩,很不易實現(xiàn)。最終選定了使用GDAL進行柵格轉矢量的方法,代碼比較簡單。
原始tif影像(12.8MB)如下:
![](/d/20211017/02696d3abb062aea40661efc45a71af3.gif)
核心函數(shù)
GDAL中柵格轉矢量的函數(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
將每個像元轉成一個矩形。
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) # 讀取路徑中的柵格數(shù)據(jù)
inband = inraster.GetRasterBand(1) # 這個波段就是最后想要轉為矢量的波段,如果是單波段數(shù)據(jù)的話那就都是1
prj = osr.SpatialReference()
prj.ImportFromWkt(inraster.GetProjection()) # 讀取柵格數(shù)據(jù)的投影信息,用來為后面生成的矢量做準備
outshp = path[:-4] + ".shp" # 給后面生成的矢量準備一個輸出文件名,這里就是把原柵格的文件名后綴名改成shp了
drv = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(outshp): # 若文件已經(jīng)存在,則刪除它繼續(xù)重新做一遍
drv.DeleteDataSource(outshp)
Polygon = drv.CreateDataSource(outshp) # 創(chuàng)建一個目標文件
Poly_layer = Polygon.CreateLayer(path[:-4], srs=prj, geom_type=ogr.wkbMultiPolygon) # 對shp文件創(chuàng)建一個圖層,定義為多個面類
newField = ogr.FieldDefn('value', ogr.OFTReal) # 給目標shp文件添加一個字段,用來存儲原始柵格的pixel value,浮點型,
Poly_layer.CreateField(newField)
gdal.Polygonize(inband, None, Poly_layer, 0) # 核心函數(shù),執(zhí)行的就是柵格轉矢量操作
# 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) # 輸出程序運行所需時間
轉換效果
轉換之后的矢量數(shù)據(jù)有270MB,非常大,打開非???/p>
![](/d/20211017/d9c1bc1e326a03164e1de62210a55732.gif)
合并之后的矢量數(shù)據(jù)有48MB,相對第一種方法數(shù)據(jù)量大大減少
![](/d/20211017/8035afa7a3ab6450129d4c397aa09b1e.gif)
到此這篇關于python 使用GDAL實現(xiàn)柵格tif轉矢量shp的文章就介紹到這了,更多相關python柵格tif轉矢量shp內容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關文章希望大家以后多多支持腳本之家!
您可能感興趣的文章:- 在Python中用GDAL實現(xiàn)矢量對柵格的切割實例
- python使用gdal對shp讀取,新建和更新的實例
- python安裝gdal的兩種方法
- 利用python GDAL庫讀寫geotiff格式的遙感影像方法
- python 矢量數(shù)據(jù)轉柵格數(shù)據(jù)代碼實例