利用GDAL-python库为SHP格式添加Z值

最新推荐文章于 2024-04-01 17:36:15 发布

原创最新推荐文章于 2024-04-01 17:36:15 发布 · 2.1k 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#float #api

Python 专栏收录该内容

8 篇文章

订阅专栏

本文介绍如何从栅格数据中读取点的Z值，并将其写入新的SHP文件中。通过使用OGR和GDAL库，文章详细说明了处理步骤，包括读取SHP文件、获取几何信息、读取栅格数据、创建带有Z值的输出文件等。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Z值

z值通常表示一个点的高程值，高程值一般都包含在栅格数据中。本文偿试从栅格数据读出对应点的z值，再写入相应点中。

步骤

读取SHP格式，获得Geometry，再根据几何类型读出相应的点
读取栅格数据，跟据已获得相应的点，读取高程值
新建输出数据源、层。将原图层中的图层定义赋予新建的图层中
新建特定的Geometry对象，将有z值的点添加到Geometry对象中
新建Feature对象，并将原有图层的属性表、新建Geometry对象赋予新建Feature对象，

写得很啰嗦，要是有耐心看完wkb的格式，也可以直接操作wkb。

代码

driver = ogr.GetDriverByName("ESRI Shapefile") inDs = driver.Open("boundary.shp",0) inLayer = inDs.GetLayer(0) outDs = driver.CreateDataSource("boundary_z.shp") outLayer = outDs.CreateLayer('boundary',None, ogr.wkbPolygon) layerDfn = inLayer.GetLayerDefn() for i in range(layerDfn.GetFieldCount()): fieldDfn = layerDfn.GetFieldDefn(i) fieldName = fieldDfn.GetNameRef() outLayer.CreateField(fieldDfn) #得到所有的feature #得到raster中的z oRing.AddPoint(float(p[0]),float(p[1]),float(p[2])) ploygon = ogr.Geometry(ogr.wkbPolygon) ploygon.AddGeometry(oRing) newFeature.SetGeometry(ploygon) outLayer.CreateFeature(newFeature)

* 在QGIS库中的QgsRasterLayer有个很简单的方法就可以提出z值—— rasterLayer.identify(point)，也可以用GDAL-python的api读取z值