Differences between revisions 2 and 3
 ⇤ ← Revision 2 as of 2006-09-19 03:14:32 → Size: 7701 Editor: lilin Comment: ← Revision 3 as of 2006-09-19 08:54:27 → ⇥ Size: 7940 Editor: lilin Comment: Deletions are marked like this. Additions are marked like this. Line 3: Line 3: Line 116: Line 114: * 还有，我没有考虑NoDataValue的情况，NoDataValue是没有数值的地区。也就是透明地区，这些值一般都是用一些大的不能再大的值来填充的。所以有NoDataValue的时候就会出现非常壮观的景象！:-) Line 120: Line 119: Line 123: Line 120:

## 一个从DEM生成VRML的例子

```   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 #读取波段数据
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吧。

• attachment:dem.jpg

• attachment:dem2.jpg

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