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

Z值

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

步骤

1. 读取SHP格式，获得Geometry，再根据几何类型读出相应的点
2. 读取栅格数据，跟据已获得相应的点，读取高程值
3. 新建输出数据源、层。将原图层中的图层定义赋予新建的图层中
4. 新建特定的Geometry对象，将有z值的点添加到Geometry对象中
5. 新建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值