一个从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)

哈,很“赶当”吧!

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

有几点需要解释:

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

不是很突兀啊!~~~

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

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

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

反馈

如果您发现我写的东西中有问题,或者您对我写的东西有意见,请登陆这个论坛

lilin/raster2vr (last edited 2009-12-25 07:15:26 by localhost)