python 使用GDAL实现栅格tif转矢量shp的方式小结
作者:GIS开发者 时间:2021-10-02 07:13:56
前言
目前有一张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,相对第一种方法数据量大大减少
来源:https://blog.csdn.net/GISuuser/article/details/119452042
![](/images/zang.png)
![](/images/jiucuo.png)
猜你喜欢
shell命令行,一键创建 python 模板文件脚本方法
网页设计十大诀窍
Oracle数据库集复制方法浅议
Go语言流程控制详情
说说回车键触发表单提交的问题
Python configparser模块操作代码实例
Python的几个高级语法概念浅析(lambda表达式闭包装饰器)
Laravel框架路由管理简单示例
URL编码“陷阱”
关于Keras模型可视化教程及关键问题的解决
![](https://img.aspxhome.com/file/2023/1/78771_0s.jpg)
python实现根据图标提取分类应用程序实例
基于JS实现动态跟随特效的示例代码
![](https://img.aspxhome.com/file/2023/5/55865_0s.png)
python交互式图形编程实例(一)
Python面向对象之类和对象实例详解
![](https://img.aspxhome.com/file/2023/2/86672_0s.png)
巧用正则表达式获取新闻中图片地址
可爱动态背景输入框
python 3.10上如何安装pyqt5
![](https://img.aspxhome.com/file/2023/6/81866_0s.png)
python XlsxWriter模块创建aexcel表格的实例讲解
![](https://img.aspxhome.com/file/2023/3/65093_0s.jpg)