用VTK绘制DEM之一

来源：互联网发布：数控车图片编程编辑：程序博客网时间：2024/06/04 13:28

DEM的显示是GIS的基本功能，利用GDAL和VTK强大的功能可以很容易做到这一点，而Python则使实现更为简单。

关于GDAL可以参考lilin的学习笔记http://wiki.woodpecker.org.cn/moin/lilin/gdal-index，VTK可以参考我博客里的其他文章。

先看看我的效果图，目前只显示地形，没有根据高程设置颜色，也没有叠加遥感影像。高程数据采用的是SRTM3，可以去http://www2.jpl.nasa.gov/srtm/cbanddataproducts.html下载。

先看代码再解释几个关键函数。

#!C:/Python24/python.exe

#coding=utf-8

#app.py

#author:liujunzhi

#date:20071210

import vtk

from vtk.wx.wxVTKRenderWindow import wxVTKRenderWindow

import sys, time

from gdal import gdal

from wxPython.wx import *

VTK_POLYGON = 7

def display_dem(renWin, demName):

#open dem file with gdal

dataset = gdal.Open(demName)

if dataset is None:

print 'open DEM file failed!'

sys.exit(0)

geoTrans = dataset.GetGeoTransform()

band = dataset.GetRasterBand(1)

#ReadAsArray(self, xoff=0, yoff=0, win_xsize=None, win_ysize=None, buf_xsize=None, buf_ysize=None, buf_obj=None)

array = band.ReadAsArray(0, 0, 500, 500)

width = len(array[0])

height = len(array)

polyData = vtk.vtkPolyData()

polyData.SetLines(vtk.vtkCellArray())

polyData.SetPolys(vtk.vtkCellArray())

points = vtk.vtkPoints()

for yPixel in range(height):

for xPixel in range(width):

xGeo = geoTrans[0] + xPixel*geoTrans[1] + yPixel*geoTrans[2]

yGeo = geoTrans[3] + xPixel*geoTrans[4] + yPixel*geoTrans[5]

#print "X:%f Y:%f Z:%f" % (xGeo, yGeo, array[yPixel][xPixel])

points.InsertNextPoint(xGeo, yGeo, array[yPixel][xPixel]/10000)

#set points

polyData.SetPoints(points)

#set polygon info

ids = vtk.vtkIdList()

for yPixel in range(height-1):

for xPixel in range(width-1):

indexTL = yPixel*width + xPixel

indexTR = yPixel*width + xPixel + 1

indexBL = indexTL + width

indexBR = indexTR + width

ids.Reset()

ids.InsertNextId(indexTL)

ids.InsertNextId(indexBL)

ids.InsertNextId(indexBR)

cell = polyData.InsertNextCell(VTK_POLYGON, ids)

ids.Reset()

ids.InsertNextId(indexTL)

ids.InsertNextId(indexTR)

ids.InsertNextId(indexBR)

cell = polyData.InsertNextCell(VTK_POLYGON, ids)

#end vtkPolyData

polyDataMapper = vtk.vtkPolyDataMapper()

polyDataMapper.SetInput(polyData)

actor = vtk.vtkActor()

actor.SetMapper(polyDataMapper)

actor.GetProperty().SetColor(1, 1, 0)

ren = vtk.vtkRenderer()

ren.AddActor(actor)

renWin.AddRenderer(ren)

#----------------------------------------------------------------------------

def main():

# every wx app needs an app

app = wxPySimpleApp()

# create the widget

frame = wxFrame(None, -1, "3d Demo", size=wxSize(400,400))

widget = wxVTKRenderWindow(frame, -1)

renWin = widget.GetRenderWindow()

demname = 'data/36074346'

display_dem(renWin, demname)

frame.Show(1)

app.MainLoop()

if __name__ == "__main__":

main()

解释几个关键的函数：

1.GDAL中的ReadAsArray

ReadAsArray(self, xoff=0, yoff=0, win_xsize=None, win_ysize=None, buf_xsize=None , buf_ysize=None, buf_obj=None)

这里引用lilin的解释：

>>> band.ReadAsArray(100,100,10,10)

头两个100是取值窗口的左上角在实际数据中所处象元的xy位置。后两个是取值窗口覆盖的区域大小，再后面是取值窗口取出数组进行缩放后数组的大小。这里需要注意的是这里的buffer大小是根据参数自动分配的，可以不指定，如果不指定，则和第3，4个参数一致。经过5，6两个参数的设置，可以进行缩放。比如真的取值窗口大小可以是20*20，而取出后数组就可以人工设置大小。让他称为10*10的数组，或者是40*40的数组。如果设置20*40则取出的数组对于真实图像来说有了变形。

2.VTK的vtkPolyData

这个类是用来组织数据的，分两部分：

1.vtkPoints

用来存放数据集中的点数据，可通过SetPoints函数指定

如：polyData.SetPoints(points)

2.vtkCellArray

用来指定这些点数据是怎么组成线、多边形的

比如points中有5个点，我指定前三个点组成一个三角形，后两个点组成一条线，代码如下：

import vtk

VTK_LINE = 3

VTK_POLYGON = 7

# initialize polydata object

polydata = vtk.vtkPolyData()

polydata.SetLines(vtk.vtkCellArray())

polydata.SetPolys(vtk.vtkCellArray())

points = vtk.vtkPoints()

# points for poly 1 (triangle)

points.InsertNextPoint(0.2, 0.3, 0)

points.InsertNextPoint(0.4, 0.55, 0)

points.InsertNextPoint(0.2, 0.8, 0)

# points for line 1 (line)

points.InsertNextPoint(0.5, 0.3, 0)

points.InsertNextPoint(0.5, 0.8, 0)

# create poly 1

ids = vtk.vtkIdList()

ids.InsertNextId(0)

ids.InsertNextId(1)

ids.InsertNextId(2)

cell = polydata.InsertNextCell(VTK_POLYGON, ids)

# create line 1

ids.Reset()

ids.InsertNextId(3)

ids.InsertNextId(4)

cell = polydata.InsertNextCell(VTK_LINE, ids)