1. 一个从DEM生成VRML的例子

现在我们牛刀小试一下，玩一个例子吧，用gdal可以读取那么多格式，不拿出来做点东西玩玩岂不是太浪费了？想当年我玩VR的时候，有过一段痛苦的经历，当时我为我们师大建VR的地形模型。首先把等高线描好，然后生成TIN，通过TIN转成GRID，然后通过GRID转成VRML。当时用的是ArcGIS8，然后转出的VRML居然是1.0，晕了。搞来搞去不行。最后花了几百兆空间装了个ERDAS转出来(其中那个版本ERDAS又和那台机器上的ArcGIS冲突，装了卸，卸了装，反复好几次，花了好几天)，最后勉强可以用了，但是用ERDAS转出的VR没有纹理映射(据说有个同学在后来研究了一下ERDAS，可以生成纹理映射-_-!，老天啊，为什么这么玩我!)，转出的东东又大，又难看。我又只好研究了几天那个文件，然后自己开C＋＋写了个映射转出程序，生成映射表，然后拷贝到wrl文件中，终于能用了。最后看效果，仿佛贴图贴得有点歪……

真是艰苦卓绝啊……(这是两年前的事了，具体过程也有点忘了，可能有些生成步骤错了，但是艰苦卓绝是一定的)

咳，要是当时知道有python，有gdal，那该多好啊！昨天晚上只花了两个小时，就搞定。这下好了，头不晕了，眼不花了，睡觉也踏实了，python＋gdal，效果好，价格便宜还实惠，我一直用它……那是相当高效！*o*

干脆贴出来吧，反正也不多，贴出来方便大家，你好，我也好。大家好才是真的好。*_*

Toggle line numbers 1#-*- coding: utf8 -*-

2"""

3该代码是进行从DEM到VRML的转换

4$ writer:lilin; create:2006.9.19; version:0.1 $

5"""

6import gdal

7import math,os

8from Numeric import *

9#打开dem

10#ds = gdal.Open("J:/arcgis/ArcTutor/Catalog/Yellowstone/dem30")

11ds = gdal.Open("J:/gisdata/gtif/dem.tif")

12#读取波段信息

13band = ds.GetRasterBand(1)

14print band.XSize,band.YSize

15trans = ds.GetGeoTransform()

16bandXSize = band.XSize

17bandYSize = band.YSize

18if bandXSize>300:#压缩像素到最大300个像素，太大了浏览器跑不动

19 scale = 300.0/bandXSize

20 bandXSize = int(bandXSize*scale)

21 bandYSize = int(bandYSize*scale)

22else:scale=1

23print bandXSize,bandYSize

24xres = trans[1]/scale

25yres = math.fabs(trans[5])/scale

26print xres,yres

27#读取波段数据

28elevs = band.ReadAsArray(0,0,band.XSize,band.YSize,bandXSize,bandYSize)

29#包装成Shape节点

30def getShape(shape):

31 return "Shape{\n%s\n}" % shape

32#包装成ElevationGrid节点

33def getElevationGrid(heights,**keymap):

34 pstr = ["%s %s" % (k,str(v)) for k,v in keymap.items()]

35 prep = "\n".join(pstr)+"\n"

36 elevs = "\n\t".join([", ".join([str(j) for j in i]) for i in heights])

37 h = "height [\n\t%s\n]" % elevs

38 context = prep+h

39 return "geometry ElevationGrid{\n%s\n}" % context

40#包装成TexCoord节点的纹理映射坐标

41def getCoorArray(w,h):

42 warr = arange(0.0, 1.0, 1.0/(w+1))[1:]

43 harr = arange(1.0, 0.0, -1.0/(h+1))[1:]

44 points = "\n\t".join([",".join(["%s %s" % (str(i),str(j)) for i in warr])for j in harr])

45 return "point[\n\t%s\n]" % points

46#包装成TexCoord节点

47def getTexCoord(w,h):

48 points = getCoorArray(w,h)

49 return "TextureCoordinate{\n%s\n}" % points

50#包装材质节点

51def getAppearance(imgname):

52 return """

53appearance DEF DEMCOLOR Appearance {

54material Material {

55diffuseColor 0.1843 0.5 0.2

56}

57texture ImageTexture {

58url ["%s"]

59repeatS TRUE

60repeatT TRUE

61}

62}

63""" % imgname

64#开始进行转换

65imgname = "dem.jpg"

66vrname = "dem.wrl"

67if os.access(vrname,os.R_OK):os.remove(vrname)

68#打开文件输入vrml文本

69file = open(vrname,'w')

70file.write("#VRML V2.0 utf8\n")

71appeartext = getAppearance(imgname)

72elevtext = getElevationGrid(elevs,

73 xDimension=bandXSize,

74 zDimension= bandYSize,

75 xSpacing=xres,

76 zSpacing=yres,

77 solid='TRUE',

78 creaseAngle=0.0,

79 texCoord = getTexCoord(bandXSize,bandYSize))

80shapetext = "\n".join([appeartext,elevtext])

81context = getShape(shapetext)

82file.write(context)

83file.close()

84#拷贝图像成jpg

85#这里有一些问题，如果不是Byte的数据类型就会出错。

86#若出错，可以用其他的软件手动改，或者再写一段代码

87#其实这些都没有关系，有了映射坐标，随便什么图都好，要的就是那张皮。

88#不过支持的格式有限，只有jpg，gif等少数格式

89if os.access(imgname,os.R_OK):os.remove(imgname)

90format = "jpeg"

91driver = gdal.GetDriverByName( format )

92dst_ds = driver.CreateCopy(imgname,ds)

哈，很“赶当”吧！

要说明的是，这断代码其实比较乱，我还没有整。而且有些地方会出一些小问题。

有几点需要解释:第18到21是万不得已的代码，太大了浏览器实在跑不动，那不仅是会卡的问题了。宽500像素的文件大小已经到3M的水平，wrl毕竟是文本，大小是呈现几何增长的，并且除了地形，还有纹理映射，又是一堆数值，大小要乘以2的，这么大，谁受得了！所以我将一张图的像素压缩到300像素内(X方向，Y方向的应该和X方向宽度差不多，如果你的图细长细长的，那你自己改代码吧，看看压缩到什么程度比较合适)，普通的机器性能上看应该还行，但这样的文件规模已经到900k了。

另外，最后一段是图像装出代码，一般情况下会出错(dem很少把数据限制在Byte类型，因为是高程数据，所以一般不会只在256范围内，比如被我注释掉的那个ArcGIS数据dem30就会出错，当时地形转出是可以的)。但幸运的是：我们在大多数情况下并不需要dem转出的图(dem转出的图一般都是黑乎乎不太美观的)。我们一般是找张遥感图或者纸质扫描图来蒙在地形上。所以，只要把那些图片转化成jpeg(vr似乎只支持jpg，gif这几种有限的格式)，命名为dem.jpg(不改名，就要开wrl改材质的贴图路径，路径设置代码在58行，“url”节点)放在wrl所在目录下就好。89~92如果总是有问题，那不如不要的好。如果你看它不顺眼，注释掉好了。

还有，我没有考虑NoDataValue的情况，NoDataValue是没有数值的地区。也就是透明地区，这些值一般都是用一些大的不能再大的值来填充的。所以有NoDataValue的时候就会出现非常壮观的景象！:-)

最后，生成的wrl文件没有缩进(要缩进也可以，代码就要多写了)，还好，现在的vr编写工具很多都可以format。看不惯？自己format吧。

效果图如下：这dem.tif是grass的demo数据集中的一个dem。我把它转出来成GeoTiff了。

不是很突兀啊！~~~

你以为，又不是青藏高原的边缘地区，哪里会那么突兀！ 一定要看突兀的地形？好吧，怕了你了。打开wrl，找到zSpacing，xSpacing 两个节点，改小点，我的是100改20。

怎么样？该凸的凸该凹的凹了吧。你也可以改代码86，87行，乘个百分几，地形就凸起来了。

当然，这样的地形，要walk是不行的。fly也够呛(单元太大)。要置身地形中，还需要建立另一个vr文件来包含这个地形文件。然后设视点，设导航属性……这些事就放在cosmo world里面弄吧。怕花钱？那就用White Dune调整也好。或者自己写，然后用inline包含dem.wrl进来(这就是下下策了)。

2. 反馈

如果您发现我写的东西中有问题，或者您对我写的东西有意见，请一定要发邮件跟我讲，Email(linux_23@163.com)