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

　　今天是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编程，稍后添加……