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

最新推荐文章于 2023-08-16 10:12:03 发布

aijingyi

最新推荐文章于 2023-08-16 10:12:03 发布

阅读量1.8k

点赞数 1

分类专栏： World Wind 文章标签： null float exception string qt download

World Wind 专栏收录该内容

24 篇文章 2 订阅

订阅专栏

原文转自：http://www.cnblogs.com/wuhenke/archive/2010/01/01/1637544.html

今天是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行构造函数中代码，

Code highlighting produced by Actipro CodeHighlighter (freeware)
http://www.CodeHighlighter.com/

-->   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,
,
, -90, -180, 180,
                    true,
                    ias);

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

QuadTileSet的Update()方法 




Code highlighting produced by Actipro CodeHighlighter (freeware)
http://www.CodeHighlighter.com/

-->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()代码 




Code highlighting produced by Actipro CodeHighlighter (freeware)
http://www.CodeHighlighter.com/

-->    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) 




Code highlighting produced by Actipro CodeHighlighter (freeware)
http://www.CodeHighlighter.com/

-->        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) 




Code highlighting produced by Actipro CodeHighlighter (freeware)
http://www.CodeHighlighter.com/

-->        /// <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数据创建高程格网的。

Code highlighting produced by Actipro CodeHighlighter (freeware)
http://www.CodeHighlighter.com/

-->        /// <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()

创建具有高程值的格网 




Code highlighting produced by Actipro CodeHighlighter (freeware)
http://www.CodeHighlighter.com/

-->        /// <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行代码

Code highlighting produced by Actipro CodeHighlighter (freeware)http://www.CodeHighlighter.com/-->  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)代码 


Code highlighting produced by Actipro CodeHighlighter (freeware)http://www.CodeHighlighter.com/-->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
配置内容：

Code highlighting produced by Actipro CodeHighlighter (freeware)http://www.CodeHighlighter.com/-->?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)方法。

获取特定点的高程值 


Code highlighting produced by Actipro CodeHighlighter (freeware)http://www.CodeHighlighter.com/-->    /// <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)方法。

Code highlighting produced by Actipro CodeHighlighter (freeware)http://www.CodeHighlighter.com/-->    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文件高程数据 


Code highlighting produced by Actipro CodeHighlighter (freeware)http://www.CodeHighlighter.com/-->        /// <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编程，稍后添加……

aijingyi

关注

1
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
WorldWind学习系列十四：DEM数据加载和应用——以SRTM为例

原文转自：http://www.cnblogs.com/wuhenke/archive/2010/01/01/1637544.html今天是2010的第一天，总想把它过得充实点，为我自己新的一年开个好头吧！首先，向关注我博客的网友道声：“元旦快乐！”，其次，接着跟大家分享一下我学习WW中DEM数据的加载和应用心得，希望大家从中有所收获！　　DEM应用在WW的三维表现中占有很重要的位置，跟影
复制链接

扫一扫

专栏目录