== 一个从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*


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

{{{
#!python
#-*- coding: utf8 -*-
"""
该代码是进行从DEM到VRML的转换
$ writer:lilin; create:2006.9.19; version:0.1 $
"""
import gdal
import math,os
from Numeric import *
#打开dem
#ds = gdal.Open("J:/arcgis/ArcTutor/Catalog/Yellowstone/dem30")
ds = gdal.Open("J:/gisdata/gtif/dem.tif")
#读取波段信息
band = ds.GetRasterBand(1)
print band.XSize,band.YSize
trans = ds.GetGeoTransform()
bandXSize = band.XSize
bandYSize = band.YSize
if bandXSize>300:#压缩像素到最大300个像素，太大了浏览器跑不动
    scale = 300.0/bandXSize
    bandXSize = int(bandXSize*scale)
    bandYSize = int(bandYSize*scale)
else:scale=1
print bandXSize,bandYSize
xres = trans[1]/scale
yres = math.fabs(trans[5])/scale
print xres,yres
#读取波段数据
elevs = band.ReadAsArray(0,0,band.XSize,band.YSize,bandXSize,bandYSize)
#包装成Shape节点
def getShape(shape):
    return "Shape{\n%s\n}" % shape
#包装成ElevationGrid节点
def getElevationGrid(heights,**keymap):
    pstr = ["%s %s" % (k,str(v)) for k,v in keymap.items()]
    prep = "\n".join(pstr)+"\n"
    elevs = "\n\t".join([", ".join([str(j) for j in i]) for i in heights])
    h = "height [\n\t%s\n]" % elevs
    context = prep+h
    return "geometry ElevationGrid{\n%s\n}" % context
#包装成TexCoord节点的纹理映射坐标
def getCoorArray(w,h):
    warr = arange(0.0, 1.0, 1.0/(w+1))[1:]
    harr = arange(1.0, 0.0, -1.0/(h+1))[1:]
    points = "\n\t".join([",".join(["%s %s" % (str(i),str(j)) for i in warr])for j in harr])
    return "point[\n\t%s\n]" % points
#包装成TexCoord节点
def getTexCoord(w,h):
    points = getCoorArray(w,h)
    return "TextureCoordinate{\n%s\n}" % points
#包装材质节点
def getAppearance(imgname):
    return """
    appearance DEF DEMCOLOR Appearance {
                material Material {
                        diffuseColor 0.1843 0.5 0.2
                }
                texture ImageTexture {
                        url ["%s"]
                        repeatS TRUE
                        repeatT TRUE
                }
        }
    """ % imgname
#开始进行转换
imgname = "dem.jpg"
vrname = "dem.wrl"
if os.access(vrname,os.R_OK):os.remove(vrname)
#打开文件输入vrml文本
file = open(vrname,'w')
file.write("#VRML V2.0 utf8\n")
appeartext = getAppearance(imgname)
elevtext = getElevationGrid(elevs,
                xDimension=bandXSize,
                zDimension= bandYSize,
                xSpacing=xres,
                zSpacing=yres,
                solid='TRUE',
                creaseAngle=0.0,
        texCoord = getTexCoord(bandXSize,bandYSize))
shapetext = "\n".join([appeartext,elevtext])
context = getShape(shapetext)
file.write(context)
file.close()
#拷贝图像成jpg
#这里有一些问题，如果不是Byte的数据类型就会出错。
#若出错，可以用其他的软件手动改，或者再写一段代码
#其实这些都没有关系，有了映射坐标，随便什么图都好，要的就是那张皮。
#不过支持的格式有限，只有jpg，gif等少数格式
if os.access(imgname,os.R_OK):os.remove(imgname)
format = "jpeg"
driver = gdal.GetDriverByName( format )
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了。

 . {{attachment:dem.jpg}}
不是很突兀啊！~~~

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

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

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

== 反馈 ==
如果您发现我写的东西中有问题，或者您对我写的东西有意见，请登陆[[http://groups.google.com/group/geosings?hl=zh-CN|这个论坛]]。