欢迎访问 生活随笔!

生活随笔

当前位置: 首页 > 编程资源 > 编程问答 >内容正文

编程问答

6. GDAL进行栅格转矢量

发布时间:2024/3/26 编程问答 44 豆豆
生活随笔 收集整理的这篇文章主要介绍了 6. GDAL进行栅格转矢量 小编觉得挺不错的,现在分享给大家,帮大家做个参考.

1 代码实现

参考下面代码实现
python 使用GDAL实现栅格tif转矢量shp的方式小结_python_脚本之家
https://www.jb51.net/article/219157.htm

from osgeo import gdal, ogr, osr import os import datetime import numpy as nppath = "Z_NAFP20210727.tif"if __name__ == '__main__':start_time = datetime.datetime.now()inraster = gdal.Open(path) # 读取路径中的栅格数据inband = inraster.GetRasterBand(1) # 这个波段就是最后想要转为矢量的波段,如果是单波段数据的话那就都是1prj = 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 = Noneend_time = datetime.datetime.now()print("Succeeded at", end_time)print("Elapsed Time:", end_time - start_time) # 输出程序运行所需时间

2 栅格转矢量并删除指定属性的要素

from osgeo import gdal, ogr, osr import os import datetime import numpy as npdef raster_shp(path):inraster = gdal.Open(path) # 读取路径中的栅格数据inband = inraster.GetRasterBand(1) # 这个波段就是最后想要转为矢量的波段,如果是单波段数据的话那就都是1prj = 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 = Nonedef raster_shp_delete_0(path):inraster = gdal.Open(path) # 读取路径中的栅格数据inband = inraster.GetRasterBand(1) # 这个波段就是最后想要转为矢量的波段,如果是单波段数据的话那就都是1prj = osr.SpatialReference()prj.ImportFromWkt(inraster.GetProjection()) # 读取栅格数据的投影信息,用来为后面生成的矢量做准备outshp = path[:-4] + "d0.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) # 只转矩形,不合并strValue = 0strFilter = "Value = '" + str(strValue) + "'"Poly_layer.SetAttributeFilter(strFilter)pFeatureDef = Poly_layer.GetLayerDefn()pLayerName = Poly_layer.GetName()pFieldName = "Value"pFieldIndex = pFeatureDef.GetFieldIndex(pFieldName)for pFeature in Poly_layer:pFeatureFID = pFeature.GetFID()Poly_layer.DeleteFeature(int(pFeatureFID))strSQL = "REPACK " + str(Poly_layer.GetName())Polygon.ExecuteSQL(strSQL, None, "");Poly_layer = NonePolygon.SyncToDisk()Polygon = Noneif __name__ == '__main__':path = r"I:\5GF2\1ALLimg20211130\8波段4预测结果\GF7_DLC_E116.6_N31.1_20210911_L1A0000560021-BWDPAN_ORTHO_PSH_preR.tif"start_time = datetime.datetime.now()raster_shp_delete_0(path)end_time = datetime.datetime.now()print("Succeeded at", end_time)print("Elapsed Time:", end_time - start_time) # 输出程序运行所需时间

参考文献:

看下面两个就可以了

[1] 【Python】GDAL批量栅格转矢量 - 知乎
https://zhuanlan.zhihu.com/p/210939809

[2] (39条消息) GDAL–栅格转矢量_secyb的博客-CSDN博客_gdal栅格转矢量
https://blog.csdn.net/secyb/article/details/80245157

【3】栅格结果不保存到磁盘中,Python ogr.GetDriverByName方法代碼示例 - 純淨天空
https://vimsky.com/zh-tw/examples/detail/python-method-osgeo.ogr.GetDriverByName.html

总结

以上是生活随笔为你收集整理的6. GDAL进行栅格转矢量的全部内容,希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错,欢迎将生活随笔推荐给好友。