WorldWind学习系列十四：DEM数据加载和应用——以SRTM为例

it2024-11-24 26

　　今天是2010的第一天，总想把它过得充实点，为我自己新的一年开个好头吧！首先，向关注我博客的网友道声：“元旦快乐！”，其次，接着跟大家分享一下我学习WW中DEM数据的加载和应用心得，希望大家从中有所收获！

　　DEM应用在WW的三维表现中占有很重要的位置，跟影像数据同等重要！幸好影像和DEM的加载和处理原理上几乎一致，对基于WW搞GIS三维开发来说是件好事，理解好任何一种，另一种触类旁通！前一篇，主要从功能上做了简单入门介绍，该篇将从代码级别分析WW内置的SRTM的DEM数据加载和应用，下一篇讲从二次开发角度上讲解如何处理、配置自己的影像和DEM数据。呵呵，因为DEM部分很重要，且是放假期间我也有时间，争取篇篇精彩！

　　两个缩写词介绍：因为这两个缩写词常出现，知道是什么缩写，就不觉得神秘啦！

　　SRTM：The Shuttle Radar Topography Mission (SRTM) obtained elevation data on a near-global scale to generate the most complete high-resolution digital topographic database of Earth. SRTM consisted of a specially modified radar system that flew onboard the Space Shuttle Endeavour during an 11-day mission in February of 2000.

　　NLT:NASA Learning Technologies.

　　我从BMNG.cs为例入手研究DEM的使用，当然研究瓦片影像也该从此入手，但，今天影像不是我们关注的重点。现在正式步入主题，跟我一起分析和学习代码吧！

　　BMNG.cs类144行构造函数中代码，

WorldWind. NltImageStore imageStore = new WorldWind.NltImageStore(String.Format( " bmng.topo.2004{0:D2} " , i + 1 ), " http://worldwind25.arc.nasa.gov/tile/tile.aspx " ); imageStore.DataDirectory = null ; imageStore.LevelZeroTileSizeDegrees = 36.0 ; imageStore.LevelCount = 5 ; imageStore.ImageExtension = " jpg " ; imageStore.CacheDirectory = String.Format( " {0}\\BMNG\\{1} " , m_WorldWindow.Cache.CacheDirectory, String.Format( " BMNG (Shaded) Tiled - {0}.2004 " , i + 1 )); ias = new WorldWind.ImageStore[ 1 ]; ias[ 0 ] = imageStore; m_QuadTileLayers[ 0 , i] = new WorldWind.Renderable. QuadTileSet( String.Format( " Tiled - {0}.2004 " , i + 1 ), m_WorldWindow.CurrentWorld, 0 , 90 , - 90 , - 180 , 180 , true , ias);

BMNG中的NltImageStore.cs、QuadTileSet类。这是我们关注的对象。QuadTileSet继承自RenderableObject，是要绘制渲染的对象类。关注它的562行Update（）方法、517行Initialize（）方法、 701行Render()方法。Update（）方法

QuadTileSet的Update()方法 public override void Update(DrawArgs drawArgs) { if ( ! isInitialized) Initialize(drawArgs); if (m_effectPath != null && m_effect == null ) { string errs = string .Empty; m_effect = Effect.FromFile(DrawArgs.Device, m_effectPath, null , "" , ShaderFlags.None, m_effectPool, out errs); if (errs != null && errs != string .Empty) { Log.Write(Log.Levels.Warning, " Could not load effect " + m_effectPath + " : " + errs); Log.Write(Log.Levels.Warning, " Effect has been disabled. " ); m_effectPath = null ; m_effect = null ; } } if (ImageStores[ 0 ].LevelZeroTileSizeDegrees < 180 ) { // Check for layer outside view double vrd = DrawArgs.Camera.ViewRange.Degrees; double latitudeMax = DrawArgs.Camera.Latitude.Degrees + vrd; double latitudeMin = DrawArgs.Camera.Latitude.Degrees - vrd; double longitudeMax = DrawArgs.Camera.Longitude.Degrees + vrd; double longitudeMin = DrawArgs.Camera.Longitude.Degrees - vrd; if (latitudeMax < m_south || latitudeMin > m_north || longitudeMax < m_west || longitudeMin > m_east) return ; } if (DrawArgs.Camera.ViewRange * 0.5f > Angle.FromDegrees(TileDrawDistance * ImageStores[ 0 ].LevelZeroTileSizeDegrees)) { lock (m_topmostTiles.SyncRoot) { foreach (QuadTile qt in m_topmostTiles.Values) qt.Dispose(); m_topmostTiles.Clear(); ClearDownloadRequests(); } return ; }　　　　　　　//知识点，可以看看，如何计算不可见瓦片的算法。 RemoveInvisibleTiles(DrawArgs.Camera); 下面主要是如何计算和加载瓦片式影像的，是重点，但不是这次的重点。 try { 　　　　　　　　//根据Camera所对的中心经纬度，计算中心点的行列号 int middleRow = MathEngine.GetRowFromLatitude(DrawArgs.Camera.Latitude, ImageStores[ 0 ].LevelZeroTileSizeDegrees); int middleCol = MathEngine.GetColFromLongitude(DrawArgs.Camera.Longitude, ImageStores[ 0 ].LevelZeroTileSizeDegrees);　　　　　　　　　//根据行列号，反推瓦片的四点对应的经度或纬度 double middleSouth = - 90.0f + middleRow * ImageStores[ 0 ].LevelZeroTileSizeDegrees; double middleNorth = - 90.0f + middleRow * ImageStores[ 0 ].LevelZeroTileSizeDegrees + ImageStores[ 0 ].LevelZeroTileSizeDegrees; double middleWest = - 180.0f + middleCol * ImageStores[ 0 ].LevelZeroTileSizeDegrees; double middleEast = - 180.0f + middleCol * ImageStores[ 0 ].LevelZeroTileSizeDegrees + ImageStores[ 0 ].LevelZeroTileSizeDegrees; double middleCenterLat = 0.5f * (middleNorth + middleSouth); double middleCenterLon = 0.5f * (middleWest + middleEast);　//这里存在一个算法，由中心瓦片框，向四周扩散地找相邻的瓦片矩形框。 //有兴趣的网友可以看一下，根据算法画出图来就好理解啦。（我感觉该算法对以后开发会有用的） int tileSpread = 4 ; for ( int i = 0 ; i < tileSpread; i ++ ) { for ( double j = middleCenterLat - i * ImageStores[ 0 ].LevelZeroTileSizeDegrees; j < middleCenterLat + i * ImageStores[ 0 ].LevelZeroTileSizeDegrees; j += ImageStores[ 0 ].LevelZeroTileSizeDegrees) { for ( double k = middleCenterLon - i * ImageStores[ 0 ].LevelZeroTileSizeDegrees; k < middleCenterLon + i * ImageStores[ 0 ].LevelZeroTileSizeDegrees; k += ImageStores[ 0 ].LevelZeroTileSizeDegrees) { //根据经\纬度和tileSize来计算行列号，这里LevelZeroTileSizeDegrees为第0层的瓦片大小为36度，瓦片总个数为50片 int curRow = MathEngine. GetRowFromLatitude (Angle.FromDegrees(j), ImageStores[ 0 ].LevelZeroTileSizeDegrees); int curCol = MathEngine. GetColFromLongitude (Angle.FromDegrees(k), ImageStores[ 0 ].LevelZeroTileSizeDegrees); long key = (( long )curRow << 32 ) + curCol;　　　　　　　　　　　　　　　　//如果集合m_topmostTiles已经存在QuadTile，则更新QuadTile QuadTile qt = (QuadTile)m_topmostTiles[key]; if (qt != null ) { qt.Update(drawArgs); continue ; } // Check for tile outside layer boundaries，获取外边框四点经度或纬度坐标 double west = - 180.0f + curCol * ImageStores[ 0 ].LevelZeroTileSizeDegrees; if (west > m_east) continue ; double east = west + ImageStores[ 0 ].LevelZeroTileSizeDegrees; if (east < m_west) continue ; double south = - 90.0f + curRow * ImageStores[ 0 ].LevelZeroTileSizeDegrees; if (south > m_north) continue ; double north = south + ImageStores[ 0 ].LevelZeroTileSizeDegrees; if (north < m_south) continue ;　　　　　　　　　　　　　　　　//结合中不存在，创建新的QuadTile qt = new QuadTile(south, north, west, east, 0 , this ); //判断新的QuadTile是否在可视区域中。（可以关注一下：Intersects（）方法判断矩形框相交）　　　　　　　　　　　　　　　　if (DrawArgs.Camera.ViewFrustum.Intersects(qt.BoundingBox)) { lock (m_topmostTiles.SyncRoot) m_topmostTiles.Add(key, qt); 　　　　　　　　　　　　　　　　　　//调用QuadTile的Update()方法 qt.Update(drawArgs); } } } } } catch (System.Threading.ThreadAbortException) { } catch (Exception caught) { Log.Write(caught); } }

Render()方法的关键代码为：

device.VertexFormat = CustomVertex.PositionNormalTextured.Format; foreach (QuadTile qt in m_topmostTiles.Values) qt.Render(drawArgs);

从上面可以看出，QuadTileSet可看作是QuadTile的集合，真正实现更新和渲染的是QuadTile对象。里面有影像的加载和渲染绘制，也有DEM的渲染绘制。

　　我们先看看QuadTile.cs 中Update()方法：

QuadTile的Update()代码 public virtual void Update(DrawArgs drawArgs) { if (m_isResetingCache) return ; try { double tileSize = North - South; if ( ! isInitialized) { if (DrawArgs.Camera.ViewRange * 0.5f < Angle.FromDegrees(QuadTileSet.TileDrawDistance * tileSize) && MathEngine.SphericalDistance(CenterLatitude, CenterLongitude, DrawArgs.Camera.Latitude, DrawArgs.Camera.Longitude) < Angle.FromDegrees(QuadTileSet.TileDrawSpread * tileSize * 1.25f ) && DrawArgs.Camera.ViewFrustum.Intersects(BoundingBox) ) Initialize(); } if (isInitialized && World.Settings.VerticalExaggeration != verticalExaggeration || m_CurrentOpacity != QuadTileSet.Opacity || QuadTileSet.RenderStruts != renderStruts) { 　　　　　　　　　　　//创建瓦片网格(重点) CreateTileMesh(); } if (isInitialized) { 　　　　　　　　　　//判断进入下一层的条件（ViewRange角度、球面距离、可视区域） if (DrawArgs.Camera.ViewRange < Angle.FromDegrees(QuadTileSet.TileDrawDistance * tileSize) && MathEngine.SphericalDistance(CenterLatitude, CenterLongitude, DrawArgs.Camera.Latitude, DrawArgs.Camera.Longitude) < Angle.FromDegrees(QuadTileSet.TileDrawSpread * tileSize) && DrawArgs.Camera.ViewFrustum.Intersects(BoundingBox) ) { if (northEastChild == null || northWestChild == null || southEastChild == null || southWestChild == null ) {　　　　　　　　　　　　　　　　　　//计算下一级别的四个子瓦片（重点，稍后一起看看） ComputeChildren(drawArgs); } if (northEastChild != null ) { northEastChild.Update(drawArgs); } if (northWestChild != null ) { northWestChild.Update(drawArgs); } if (southEastChild != null ) { southEastChild.Update(drawArgs); } if (southWestChild != null ) { southWestChild.Update(drawArgs); } } else { if (northWestChild != null ) { northWestChild.Dispose(); northWestChild = null ; } if (northEastChild != null ) { northEastChild.Dispose(); northEastChild = null ; } if (southEastChild != null ) { southEastChild.Dispose(); southEastChild = null ; } if (southWestChild != null ) { southWestChild.Dispose(); southWestChild = null ; } } } if (isInitialized) { if (DrawArgs.Camera.ViewRange / 2 > Angle.FromDegrees(QuadTileSet.TileDrawDistance * tileSize * 1.5f ) || MathEngine.SphericalDistance(CenterLatitude, CenterLongitude, DrawArgs.Camera.Latitude, DrawArgs.Camera.Longitude) > Angle.FromDegrees(QuadTileSet.TileDrawSpread * tileSize * 1.5f )) { if (Level != 0 || (Level == 0 && ! QuadTileSet.AlwaysRenderBaseTiles)) this .Dispose(); } } } catch { } }

创建下一级的四个瓦片的方法：（可以被我们重用）

ComputeChildren(DrawArgs drawArgs) public virtual void ComputeChildren(DrawArgs drawArgs) { 　　　　　　　//判断是否是最高级别 if (Level + 1 >= QuadTileSet.ImageStores[ 0 ].LevelCount) return ;　　　　　　　//计算瓦片的中点经纬度 double CenterLat = 0.5f * (South + North); double CenterLon = 0.5f * (East + West); if (northWestChild == null ) northWestChild = ComputeChild(CenterLat, North, West, CenterLon); if (northEastChild == null ) northEastChild = ComputeChild(CenterLat, North, CenterLon, East); if (southWestChild == null ) southWestChild = ComputeChild(South, CenterLat, West, CenterLon); if (southEastChild == null ) southEastChild = ComputeChild(South, CenterLat, CenterLon, East); }

ComputeChild(double childSouth, double childNorth, double childWest, double childEast) /// <summary> /// Returns the QuadTile for specified location if available. /// Tries to queue a download if not available. /// </summary> /// <returns> Initialized QuadTile if available locally, else null. </returns> private QuadTile ComputeChild( double childSouth, double childNorth, double childWest, double childEast) { QuadTile child = new QuadTile( childSouth, childNorth, childWest, childEast, this .Level + 1 , QuadTileSet); return child; }

QuadTile.cs 中的CreateTileMesh()方法用来创建瓦片格网的，分别在Initialize() 、Update()方法中调用。

410行这里调用的CreateElevatedMesh();是添加DEM数据创建高程格网的。

/// <summary> /// Builds flat or terrain mesh for current tile /// </summary> public virtual void CreateTileMesh() { verticalExaggeration = World.Settings.VerticalExaggeration; m_CurrentOpacity = QuadTileSet.Opacity; renderStruts = QuadTileSet.RenderStruts; if (QuadTileSet.TerrainMapped && Math.Abs(verticalExaggeration) > 1e - 3 ) 　　　　　　　　　//创建具有高程值的格网（今天要关注的） CreateElevatedMesh(); else //创建没有高程值的平面格网　　　　　　　　 CreateFlatMesh(); }

591行CreateElevatedMesh()

创建具有高程值的格网 /// <summary> /// Build the elevated terrain mesh /// </summary> protected virtual void CreateElevatedMesh() { isDownloadingTerrain = true ; //vertexCountElevated值为40；向四周分别扩充一个样本点 // Get height data with one extra sample around the tile double degreePerSample = LatitudeSpan / vertexCountElevated; //获取具有高程值的TerrainTile对象（这是最关键部分，深入分析） TerrainTile tile = QuadTileSet.World.TerrainAccessor. GetElevationArray (North + degreePerSample, South - degreePerSample, West - degreePerSample, East + degreePerSample, vertexCountElevated + 3 ); float [,] heightData = tile.ElevationData; int vertexCountElevatedPlus3 = vertexCountElevated / 2 + 3 ; int totalVertexCount = vertexCountElevatedPlus3 * vertexCountElevatedPlus3; northWestVertices = new CustomVertex.PositionNormalTextured[totalVertexCount]; southWestVertices = new CustomVertex.PositionNormalTextured[totalVertexCount]; northEastVertices = new CustomVertex.PositionNormalTextured[totalVertexCount]; southEastVertices = new CustomVertex.PositionNormalTextured[totalVertexCount]; double layerRadius = ( double )QuadTileSet.LayerRadius; // Calculate mesh base radius (bottom vertices) // Find minimum elevation to account for possible bathymetry float minimumElevation = float .MaxValue; float maximumElevation = float .MinValue; foreach ( float height in heightData) { if (height < minimumElevation) minimumElevation = height; if (height > maximumElevation) maximumElevation = height; } minimumElevation *= verticalExaggeration; maximumElevation *= verticalExaggeration; if (minimumElevation > maximumElevation) { // Compensate for negative vertical exaggeration minimumElevation = maximumElevation; maximumElevation = minimumElevation; } double overlap = 500 * verticalExaggeration; // 500m high tiles // Radius of mesh bottom grid meshBaseRadius = layerRadius + minimumElevation - overlap; CreateElevatedMesh(ChildLocation.NorthWest, northWestVertices, meshBaseRadius, heightData); CreateElevatedMesh(ChildLocation.SouthWest, southWestVertices, meshBaseRadius, heightData); CreateElevatedMesh(ChildLocation.NorthEast, northEastVertices, meshBaseRadius, heightData); CreateElevatedMesh(ChildLocation.SouthEast, southEastVertices, meshBaseRadius, heightData); BoundingBox = new BoundingBox(( float )South, ( float )North, ( float )West, ( float )East, ( float )layerRadius, ( float )layerRadius + 10000 * this .verticalExaggeration); QuadTileSet.IsDownloadingElevation = false ; // Build common set of indexes for the 4 child meshes int vertexCountElevatedPlus2 = vertexCountElevated / 2 + 2 ; vertexIndexes = new short [ 2 * vertexCountElevatedPlus2 * vertexCountElevatedPlus2 * 3 ]; int elevated_idx = 0 ; for ( int i = 0 ; i < vertexCountElevatedPlus2; i ++ ) { for ( int j = 0 ; j < vertexCountElevatedPlus2; j ++ ) { vertexIndexes[elevated_idx ++ ] = ( short )(i * vertexCountElevatedPlus3 + j); vertexIndexes[elevated_idx ++ ] = ( short )((i + 1 ) * vertexCountElevatedPlus3 + j); vertexIndexes[elevated_idx ++ ] = ( short )(i * vertexCountElevatedPlus3 + j + 1 ); vertexIndexes[elevated_idx ++ ] = ( short )(i * vertexCountElevatedPlus3 + j + 1 ); vertexIndexes[elevated_idx ++ ] = ( short )((i + 1 ) * vertexCountElevatedPlus3 + j); vertexIndexes[elevated_idx ++ ] = ( short )((i + 1 ) * vertexCountElevatedPlus3 + j + 1 ); } } calculate_normals( ref northWestVertices, vertexIndexes); calculate_normals( ref southWestVertices, vertexIndexes); calculate_normals( ref northEastVertices, vertexIndexes); calculate_normals( ref southEastVertices, vertexIndexes); isDownloadingTerrain = false ; }

596行TerrainTile tile = QuadTileSet.World.TerrainAccessor.GetElevationArray(North + degreePerSample, South - degreePerSample, West - degreePerSample, East + degreePerSample, vertexCountElevated + 3);获取样本点的高程值数组。

使用了TerrainAccessor.cs类120行代码

public virtual TerrainTile GetElevationArray( double north, double south, double west, double east, int samples) { TerrainTile res = null ; res = new TerrainTile( null ); res.North = north; res.South = south; res.West = west; res.East = east; res.SamplesPerTile = samples; res.IsInitialized = true ; res.IsValid = true ; double latrange = Math.Abs(north - south); double lonrange = Math.Abs(east - west); float [,] data = new float [samples,samples]; float scaleFactor = ( float ) 1 / (samples - 1 ); for ( int x = 0 ; x < samples; x ++ ) { for ( int y = 0 ; y < samples; y ++ ) { double curLat = north - scaleFactor * latrange * x; double curLon = west + scaleFactor * lonrange * y; //关键，获取瓦片所有样本点的高程值 data[x,y] = GetElevationAt(curLat, curLon, 0 ); } } res.ElevationData = data; return res; }

关键代码：data[x,y] = GetElevationAt(curLat, curLon, 0);

GetElevationAt();在TerrainAccessor.cs是抽象方法。真正实现是在TerrainAccessor的子类NltTerrainAccessor中重载实现的。120行 public override float GetElevationAt(double latitude, double longitude) { return GetElevationAt(latitude, longitude, m_terrainTileService.SamplesPerTile / m_terrainTileService.LevelZeroTileSizeDegrees); }

TerrainAccessor对象哪里来的（即：在哪完成初始化和传入的？）

ConfigurationLoader.cs的Load()方法的97行代码：

TerrainAccessor[] terrainAccessor = getTerrainAccessorsFromXPathNodeIterator(worldIter.Current.Select("TerrainAccessor"), System.IO.Path.Combine(cache.CacheDirectory, worldName));

World newWorld = new World( worldName, new Microsoft.DirectX.Vector3(0, 0, 0), new Microsoft.DirectX.Quaternion(0, 0, 0, 0), equatorialRadius, cache.CacheDirectory, (terrainAccessor != null ? terrainAccessor[0] : null)//TODO: Oops, World should be able to handle an array of terrainAccessors );

然后通过World对象传入到QuadTileSet类中的。

我们看看getTerrainAccessorsFromXPathNodeIterator方法如何完成完成TerrainAccessor对象。注意：该方法返回值为TerrainAccessor[]，是个数组,为什么呢？？（请关注我下一篇文章）

2078行

getTerrainAccessorsFromXPathNodeIterator(XPathNodeIterator iter, string cacheDirectory)代码 private static TerrainAccessor[] getTerrainAccessorsFromXPathNodeIterator(XPathNodeIterator iter, string cacheDirectory) { System.Collections.ArrayList terrainAccessorList = new System.Collections.ArrayList(); //下面是读取DEM配置XML，并根据配置信息创建TerrainTileService对象和TerrainAccessor对象 while (iter.MoveNext()) { string terrainAccessorName = iter.Current.GetAttribute( " Name " , "" ); if (terrainAccessorName == null ) { // TODO: Throw exception? continue ; } XPathNodeIterator latLonBoxIter = iter.Current.Select( " LatLonBoundingBox " ); if (latLonBoxIter.Count != 1 ) { // TODO: Throw exception? continue ; } double north = 0 ; double south = 0 ; double west = 0 ; double east = 0 ; latLonBoxIter.MoveNext(); north = ParseDouble(getInnerTextFromFirstChild(latLonBoxIter.Current.Select( " North " ))); south = ParseDouble(getInnerTextFromFirstChild(latLonBoxIter.Current.Select( " South " ))); west = ParseDouble(getInnerTextFromFirstChild(latLonBoxIter.Current.Select( " West " ))); east = ParseDouble(getInnerTextFromFirstChild(latLonBoxIter.Current.Select( " East " ))); TerrainAccessor[] higerResolutionSubsets = getTerrainAccessorsFromXPathNodeIterator( iter.Current.Select( " HigherResolutionSubsets " ), Path.Combine(cacheDirectory, terrainAccessorName)); XPathNodeIterator tileServiceIter = iter.Current.Select( " TerrainTileService " ); if (tileServiceIter.Count == 1 ) { string serverUrl = null ; string dataSetName = null ; double levelZeroTileSizeDegrees = double .NaN; uint numberLevels = 0 ; uint samplesPerTile = 0 ; string dataFormat = null ; string fileExtension = null ; string compressionType = null ; tileServiceIter.MoveNext(); serverUrl = getInnerTextFromFirstChild(tileServiceIter.Current.Select( " ServerUrl " )); dataSetName = getInnerTextFromFirstChild(tileServiceIter.Current.Select( " DataSetName " )); levelZeroTileSizeDegrees = ParseDouble(getInnerTextFromFirstChild(tileServiceIter.Current.Select( " LevelZeroTileSizeDegrees " ))); numberLevels = uint .Parse(getInnerTextFromFirstChild(tileServiceIter.Current.Select( " NumberLevels " ))); samplesPerTile = uint .Parse(getInnerTextFromFirstChild(tileServiceIter.Current.Select( " SamplesPerTile " ))); dataFormat = getInnerTextFromFirstChild(tileServiceIter.Current.Select( " DataFormat " )); fileExtension = getInnerTextFromFirstChild(tileServiceIter.Current.Select( " FileExtension " )); compressionType = getInnerTextFromFirstChild(tileServiceIter.Current.Select( " CompressonType " )); //根据配置信息创建TerrainTileService对象和TerrainAccessor对象 TerrainTileService tts = new TerrainTileService( serverUrl, dataSetName, levelZeroTileSizeDegrees, ( int )samplesPerTile, fileExtension, ( int )numberLevels, Path.Combine(cacheDirectory, terrainAccessorName), World.Settings.TerrainTileRetryInterval, dataFormat); TerrainAccessor newTerrainAccessor = new NltTerrainAccessor( terrainAccessorName, west, south, east, north, tts, higerResolutionSubsets); terrainAccessorList.Add(newTerrainAccessor); } // TODO: Add Floating point terrain Accessor code // TODO: Add WMSAccessor code and make it work in TerrainAccessor (which it currently doesn't) } if (terrainAccessorList.Count > 0 ) { return (TerrainAccessor[])terrainAccessorList.ToArray( typeof (TerrainAccessor)); } else { return null ; } }

再来看看DEM的配置在哪里和XML内容吧路径：C:\Program Files\NASA\World Wind 1.4\Config\Earth.xml配置内容：

?xml version="1.0" encoding="UTF-8"?> < World Name ="Earth" EquatorialRadius ="6378137.0" LayerDirectory ="Earth" xmlns:xsi ="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation ="WorldXmlDescriptor.xsd" > < TerrainAccessor Name ="SRTM" > //黄色之间的XML就是一个TerrainAccessor配置 < TerrainTileService > < ServerUrl > http://worldwind25.arc.nasa.gov/wwelevation/wwelevation.aspx </ ServerUrl > < DataSetName > srtm30pluszip </ DataSetName > < LevelZeroTileSizeDegrees > 20.0 </ LevelZeroTileSizeDegrees > < NumberLevels > 12 </ NumberLevels > < SamplesPerTile > 150 </ SamplesPerTile > < DataFormat > Int16 </ DataFormat > < FileExtension > bil </ FileExtension > < CompressonType > zip </ CompressonType > </ TerrainTileService > < LatLonBoundingBox > < North > < Value > 90.0 </ Value > </ North > < South > < Value > -90.0 </ Value > </ South > < West > < Value > -180.0 </ Value > </ West > < East > < Value > 180.0 </ Value > </ East > </ LatLonBoundingBox > </ TerrainAccessor > </ World >

接着上面的讲NltTerrainAccessor类76行代码GetElevationAt(double latitude, double longitude, double targetSamplesPerDegree)方法。

获取特定点的高程值 /// <summary> /// Get terrain elevation at specified location. /// </summary> /// <param name="latitude"> Latitude in decimal degrees. </param> /// <param name="longitude"> Longitude in decimal degrees. </param> /// <param name="targetSamplesPerDegree"></param> /// <returns> Returns 0 if the tile is not available on disk. </returns> public override float GetElevationAt( double latitude, double longitude, double targetSamplesPerDegree) { try { if (m_terrainTileService == null || targetSamplesPerDegree < World.Settings.MinSamplesPerDegree) return 0 ; if (m_higherResolutionSubsets != null ) { foreach (TerrainAccessor higherResSub in m_higherResolutionSubsets) { if (latitude > higherResSub.South && latitude < higherResSub.North && longitude > higherResSub.West && longitude < higherResSub.East) { return higherResSub.GetElevationAt(latitude, longitude, targetSamplesPerDegree); } } } //自己可以看看如何完成TerrainTile的初始化构建的 TerrainTile tt = m_terrainTileService. GetTerrainTile(latitude, longitude, targetSamplesPerDegree); TerrainTileCacheEntry ttce = (TerrainTileCacheEntry)m_tileCache[tt.TerrainTileFilePath]; if (ttce == null ) { ttce = new TerrainTileCacheEntry(tt); AddToCache(ttce); } if ( ! ttce.TerrainTile.IsInitialized) ttce.TerrainTile.Initialize(); ttce.LastAccess = DateTime.Now; //获取高程值 return ttce.TerrainTile.GetElevationAt(latitude, longitude); } catch (Exception) { } return 0 ; }

上面获取高程值的关键：TerrainTile类的330行的GetElevationAt(double latitude, double longitude)方法

public float GetElevationAt( double latitude, double longitude) { try { double deltaLat = North - latitude; double deltaLon = longitude - West;//TileSizeDegrees为当前级别下瓦片的度数大小 //计算方法：158行tile.TileSizeDegrees = m_levelZeroTileSizeDegrees * Math.Pow(0.5, tile.TargetLevel); //注意思考：SamplesPerTile-1为什么是减去1？传进来的SamplesPerTile=43而不是44？ //如果传入的是44，这里应该减2 double df2 = (SamplesPerTile - 1 ) / TileSizeDegrees; float lat_pixel = ( float )(deltaLat * df2); float lon_pixel = ( float )(deltaLon * df2); //这里是将点，近似成包含点的最小正方形（经纬度取整） int lat_min = ( int )lat_pixel; int lat_max = ( int )Math.Ceiling(lat_pixel); int lon_min = ( int )lon_pixel; int lon_max = ( int )Math.Ceiling(lon_pixel); if (lat_min >= SamplesPerTile) lat_min = SamplesPerTile - 1 ; if (lat_max >= SamplesPerTile) lat_max = SamplesPerTile - 1 ; if (lon_min >= SamplesPerTile) lon_min = SamplesPerTile - 1 ; if (lon_max >= SamplesPerTile) lon_max = SamplesPerTile - 1 ; if (lat_min < 0 ) lat_min = 0 ; if (lat_max < 0 ) lat_max = 0 ; if (lon_min < 0 ) lon_min = 0 ; if (lon_max < 0 ) lon_max = 0 ; float delta = lat_pixel - lat_min; //根据外矩形四顶点的经纬度分别插值计算中间线的高程 float westElevation = ElevationData [lon_min, lat_min] * ( 1 - delta) + ElevationData [lon_min, lat_max] * delta; float eastElevation = ElevationData[lon_max, lat_min] * ( 1 - delta) + ElevationData[lon_max, lat_max] * delta; //利用插值计算点的高程值 delta = lon_pixel - lon_min; float interpolatedElevation = westElevation * ( 1 - delta) + eastElevation * delta; return interpolatedElevation; } catch { } return 0 ; }

public float[,] ElevationData就是存放当前瓦片所有样本点高程值的数值。这是通过Initialize()中读取DEM（.bil）文件来获取的。

读取BIL文件高程数据 /// <summary> /// This method initializes the terrain tile add switches to /// Initialize floating point/int 16 tiles /// </summary> public void Initialize() { if (IsInitialized) return ; if ( ! File.Exists(TerrainTileFilePath)) { // Download elevation if (request == null ) { using ( request = new TerrainDownloadRequest( this , m_owner, Row, Col, TargetLevel) ) { request.SaveFilePath = TerrainTileFilePath; request.DownloadInForeground(); } } } if (ElevationData == null ) ElevationData = new float [SamplesPerTile, SamplesPerTile]; if (File.Exists(TerrainTileFilePath)) { // Load elevation file try { // TerrainDownloadRequest's FlagBadTile() creates empty files // as a way to flag "bad" terrain tiles. // Remove the empty 'flag' files after preset time. try { FileInfo tileInfo = new FileInfo(TerrainTileFilePath); if (tileInfo.Length == 0 ) { TimeSpan age = DateTime.Now.Subtract( tileInfo.LastWriteTime ); if (age < m_owner.TerrainTileRetryInterval) { // This tile is still flagged bad IsInitialized = true ; } else { // remove the empty 'flag' file File.Delete(TerrainTileFilePath); } return ; } } catch { // Ignore any errors in the above block, and continue. // For example, if someone had the empty 'flag' file // open, the delete would fail. } //读取BIL文件数据的关键代码，可以被我们借鉴使用 using ( Stream s = File.OpenRead(TerrainTileFilePath)) { BinaryReader reader = new BinaryReader(s); if (m_owner.DataType == " Int16 " ) { for ( int y = 0 ; y < SamplesPerTile; y ++ ) for ( int x = 0 ; x < SamplesPerTile; x ++ ) ElevationData[x,y] = reader.ReadInt16(); } if (m_owner.DataType == " Float32 " ) { for ( int y = 0 ; y < SamplesPerTile; y ++ ) for ( int x = 0 ; x < SamplesPerTile; x ++ ) { ElevationData[x,y] = reader.ReadSingle(); } } IsInitialized = true ; IsValid = true ; } return ; } catch (IOException) { // If there is an IO exception when reading the terrain tile, // then either something is wrong with the file, or with // access to the file, so try and remove it. try { File.Delete(TerrainTileFilePath); } catch (Exception ex) { throw new ApplicationException(String.Format( " Error while trying to delete corrupt terrain tile {0} " , TerrainTileFilePath), ex); } } catch (Exception ex) { // Some other type of error when reading the terrain tile. throw new ApplicationException(String.Format( " Error while trying to read terrain tile {0} " , TerrainTileFilePath), ex); } } }

另注：SRTM的高程数据存放路径如下图：（DEM跟影像也是分级存放的，存放方式一致）

至此，如何加载DEM数据创建网格的过程已经分析完了。

接下来，继续分析QuadTile.cs中CreateElevatedMesh()和Render方法，内容主要是DirectX编程，稍后添加……

转载于:https://www.cnblogs.com/wuhenke/archive/2010/01/01/1637544.html

相关资源：数据结构—成绩单生成器

最新回复(0)