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

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值