python 使用GDAL實(shí)現(xiàn)柵格tif轉(zhuǎn)矢量shp的方式小結(jié)
目前有一張tif格式的柵格影像,需要在web地圖上進(jìn)行展示,使用動(dòng)態(tài)切片WMS的方式,渲染速度比較慢,而且大的時(shí)候會(huì)出現(xiàn)模糊的問題。并且后面需要做多期影像的切換,渲染與加載效率也值得關(guān)注。
計(jì)劃是使用柵格轉(zhuǎn)矢量的方式,將柵格數(shù)據(jù)轉(zhuǎn)為矢量shp文件,然后進(jìn)行矢量切片,使用Mapbox進(jìn)行前端動(dòng)態(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)如下:
GDAL中柵格轉(zhuǎn)矢量的函數(shù)主要是以下兩個(gè),二者的參數(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
將每個(gè)像元轉(zhuǎn)成一個(gè)矩形。
Polygonize(*args, **kwargs) **
Polygonize(Band srcBand, Band maskBand, Layer outLayer, int iPixValField, char ** options=None, GDALProgressFunc callback=0, void * callback_data=None) -> int
將每個(gè)像元轉(zhuǎn)成一個(gè)矩形,然后將相似的像元進(jìn)行合并。
轉(zhuǎn)換代碼from osgeo import gdal, ogr, osrimport osimport datetimeimport numpy as nppath = 'Z_NAFP20210727.tif'if __name__ == ’__main__’: start_time = datetime.datetime.now() inraster = gdal.Open(path) # 讀取路徑中的柵格數(shù)據(jù) inband = inraster.GetRasterBand(1) # 這個(gè)波段就是最后想要轉(zhuǎn)為矢量的波段,如果是單波段數(shù)據(jù)的話那就都是1 prj = osr.SpatialReference() prj.ImportFromWkt(inraster.GetProjection()) # 讀取柵格數(shù)據(jù)的投影信息,用來為后面生成的矢量做準(zhǔn)備 outshp = path[:-4] + '.shp' # 給后面生成的矢量準(zhǔn)備一個(gè)輸出文件名,這里就是把原柵格的文件名后綴名改成shp了 drv = ogr.GetDriverByName('ESRI Shapefile') if os.path.exists(outshp): # 若文件已經(jīng)存在,則刪除它繼續(xù)重新做一遍drv.DeleteDataSource(outshp) Polygon = drv.CreateDataSource(outshp) # 創(chuàng)建一個(gè)目標(biāo)文件 Poly_layer = Polygon.CreateLayer(path[:-4], srs=prj, geom_type=ogr.wkbMultiPolygon) # 對shp文件創(chuàng)建一個(gè)圖層,定義為多個(gè)面類 newField = ogr.FieldDefn(’value’, ogr.OFTReal) # 給目標(biāo)shp文件添加一個(gè)字段,用來存儲(chǔ)原始柵格的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)行所需時(shí)間轉(zhuǎn)換效果 使用FPolygonize
轉(zhuǎn)換之后的矢量數(shù)據(jù)有270MB,非常大,打開非常卡
合并之后的矢量數(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)容請搜索好吧啦網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持好吧啦網(wǎng)!
相關(guān)文章:
1. IntelliJ IDEA設(shè)置默認(rèn)瀏覽器的方法2. 簡述JAVA同步、異步、阻塞和非阻塞之間的區(qū)別3. Python TestSuite生成測試報(bào)告過程解析4. 在JSP中使用formatNumber控制要顯示的小數(shù)位數(shù)方法5. SpringBoot項(xiàng)目優(yōu)雅的全局異常處理方式(全網(wǎng)最新)6. docker /var/lib/docker/aufs/mnt 目錄清理方法7. IntelliJ IDEA設(shè)置背景圖片的方法步驟8. 如何清空python的變量9. 解決python路徑錯(cuò)誤,運(yùn)行.py文件,找不到路徑的問題10. python操作數(shù)據(jù)庫獲取結(jié)果之fetchone和fetchall的區(qū)別說明
