一个从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*
干脆贴出来吧,反正也不多,贴出来方便大家,你好,我也好。大家好才是真的好。*_*
1 #-*- coding: utf8 -*-
2 """
3 该代码是进行从DEM到VRML的转换
4 $ writer:lilin; create:2006.9.19; version:0.1 $
5 """
6 import gdal
7 import math,os
8 from Numeric import *
9 #打开dem
10 #ds = gdal.Open("J:/arcgis/ArcTutor/Catalog/Yellowstone/dem30")
11 ds = gdal.Open("J:/gisdata/gtif/dem.tif")
12 #读取波段信息
13 band = ds.GetRasterBand(1)
14 print band.XSize,band.YSize
15 trans = ds.GetGeoTransform()
16 bandXSize = band.XSize
17 bandYSize = band.YSize
18 if bandXSize>300:#压缩像素到最大300个像素,太大了浏览器跑不动
19 scale = 300.0/bandXSize
20 bandXSize = int(bandXSize*scale)
21 bandYSize = int(bandYSize*scale)
22 else:scale=1
23 print bandXSize,bandYSize
24 xres = trans[1]/scale
25 yres = math.fabs(trans[5])/scale
26 print xres,yres
27 #读取波段数据
28 elevs = band.ReadAsArray(0,0,band.XSize,band.YSize,bandXSize,bandYSize)
29 #包装成Shape节点
30 def getShape(shape):
31 return "Shape{\n%s\n}" % shape
32 #包装成ElevationGrid节点
33 def 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节点的纹理映射坐标
41 def 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节点
47 def getTexCoord(w,h):
48 points = getCoorArray(w,h)
49 return "TextureCoordinate{\n%s\n}" % points
50 #包装材质节点
51 def getAppearance(imgname):
52 return """
53 appearance DEF DEMCOLOR Appearance {
54 material Material {
55 diffuseColor 0.1843 0.5 0.2
56 }
57 texture ImageTexture {
58 url ["%s"]
59 repeatS TRUE
60 repeatT TRUE
61 }
62 }
63 """ % imgname
64 #开始进行转换
65 imgname = "dem.jpg"
66 vrname = "dem.wrl"
67 if os.access(vrname,os.R_OK):os.remove(vrname)
68 #打开文件输入vrml文本
69 file = open(vrname,'w')
70 file.write("#VRML V2.0 utf8\n")
71 appeartext = getAppearance(imgname)
72 elevtext = 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))
80 shapetext = "\n".join([appeartext,elevtext])
81 context = getShape(shapetext)
82 file.write(context)
83 file.close()
84 #拷贝图像成jpg
85 #这里有一些问题,如果不是Byte的数据类型就会出错。
86 #若出错,可以用其他的软件手动改,或者再写一段代码
87 #其实这些都没有关系,有了映射坐标,随便什么图都好,要的就是那张皮。
88 #不过支持的格式有限,只有jpg,gif等少数格式
89 if os.access(imgname,os.R_OK):os.remove(imgname)
90 format = "jpeg"
91 driver = gdal.GetDriverByName( format )
92 dst_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进来(这就是下下策了)。
反馈
如果您发现我写的东西中有问题,或者您对我写的东西有意见,请登陆这个论坛。