java离散点生成dem_WorldWind学习系列十四：DEM数据加载和应用——以SRTM为例

最新推荐文章于 2023-02-03 10:18:36 发布

天天测评

最新推荐文章于 2023-02-03 10:18:36 发布

阅读量569

点赞数

文章标签： java离散点生成dem

本文链接：https://blog.csdn.net/weixin_35282700/article/details/114709261

版权

今天是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.NltImageStoreimageStore=newWorldWind.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=newWorldWind.ImageStore[1];

ias[0]=imageStore;

m_QuadTileLayers[0, i]=newWorldWind.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()方法

publicoverridevoidUpdate(DrawArgs drawArgs)

{if(!isInitialized)

Initialize(drawArgs);if(m_effectPath!=null&&m_effect==null)

{stringerrs=string.Empty;

m_effect=Effect.FromFile(DrawArgs.Device, m_effectPath,null,"", ShaderFlags.None, m_effectPool,outerrs);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 viewdoublevrd=DrawArgs.Camera.ViewRange.Degrees;doublelatitudeMax=DrawArgs.Camera.Latitude.Degrees+vrd;doublelatitudeMin=DrawArgs.Camera.Latitude.Degrees-vrd;doublelongitudeMax=DrawArgs.Camera.Longitude.Degrees+vrd;doublelongitudeMin=DrawArgs.Camera.Longitude.Degrees-vrd;if(latitudeMaxm_north||longitudeMaxm_east)return;

}if(DrawArgs.Camera.ViewRange*0.5f>Angle.FromDegrees(TileDrawDistance*ImageStores[0].LevelZeroTileSizeDegrees))

{lock(m_topmostTiles.SyncRoot)

{foreach(QuadTile qtinm_topmostTiles.Values)

qt.Dispose();

m_topmostTiles.Clear();

ClearDownloadRequests();

}return;

}

//知识点，可以看看，如何计算不可见瓦片的算法。

RemoveInvisibleTiles(DrawArgs.Camera);

下面主要是如何计算和加载瓦片式影像的，是重点，但不是这次的重点。try{

//根据Camera所对的中心经纬度，计算中心点的行列号intmiddleRow=MathEngine.GetRowFromLatitude(DrawArgs.Camera.Latitude, ImageStores[0].LevelZeroTileSizeDegrees);intmiddleCol=MathEngine.GetColFromLongitude(DrawArgs.Camera.Longitude, ImageStores[0].LevelZeroTileSizeDegrees);

//根据行列号，反推瓦片的四点对应的经度或纬度doublemiddleSouth=-90.0f+middleRow*ImageStores[0].LevelZeroTileSizeDegrees;doublemiddleNorth=-90.0f+middleRow*ImageStores[0].LevelZeroTileSizeDegrees+ImageStores[0].LevelZeroTileSizeDegrees;doublemiddleWest=-180.0f+middleCol*ImageStores[0].LevelZeroTileSizeDegrees;doublemiddleEast=-180.0f+middleCol*ImageStores[0].LevelZeroTileSizeDegrees+ImageStores[0].LevelZeroTileSizeDegrees;doublemiddleCenterLat=0.5f*(middleNorth+middleSouth);doublemiddleCenterLon=0.5f*(middleWest+middleEast);

//这里存在一个算法，由中心瓦片框，向四周扩散地找相邻的瓦片矩形框。

//有兴趣的网友可以看一下，根据算法画出图来就好理解啦。(我感觉该算法对以后开发会有用的)inttileSpread=4;for(inti=0; i

{for(doublej=middleCenterLat-i*ImageStores[0].LevelZeroTileSizeDegrees; j

{for(doublek=middleCenterLon-i*ImageStores[0].LevelZeroTileSizeDegrees; k

{

//根据经\纬度和tileSize来计算行列号，这里LevelZeroTileSizeDegrees为第0层的瓦片大小为36度，瓦片总个数为50片intcurRow=MathEngine.GetRowFromLatitude(Angle.FromDegrees(j), ImageStores[0].LevelZeroTileSizeDegrees);intcurCol=MathEngine.GetColFromLongitude(Angle.FromDegrees(k), ImageStores[0].LevelZeroTileSizeDegrees);longkey=((long)curRow<<32)+curCol;

//如果集合m_topmostTiles已经存在QuadTile，则更新QuadTileQuadTile qt=(QuadTile)m_topmostTiles[key];if(qt!=null)

{qt.Update(drawArgs);continue;

}//Check for tile outside layer boundaries，获取外边框四点经度或纬度坐标doublewest=-180.0f+curCol*ImageStores[0].LevelZeroTileSizeDegrees;if(west>m_east)continue;doubleeast=west+ImageStores[0].LevelZeroTileSizeDegrees;if(eastm_north)continue;doublenorth=south+ImageStores[0].LevelZeroTileSizeDegrees;if(north

//结合中不存在，创建新的QuadTileqt=newQuadTile(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()代码

publicvirtualvoidUpdate(DrawArgs drawArgs)

{if(m_isResetingCache)return;try{doubletileSize=North-South;if(!isInitialized)

{if(DrawArgs.Camera.ViewRange*0.5f

DrawArgs.Camera.Latitude, DrawArgs.Camera.Longitude)

)

Initialize();

}if(isInitialized&&World.Settings.VerticalExaggeration!=verticalExaggeration||m_CurrentOpacity!=QuadTileSet.Opacity||QuadTileSet.RenderStruts!=renderStruts)

{

//创建瓦片网格(重点)CreateTileMesh();}if(isInitialized)

{

//判断进入下一层的条件(ViewRange角度、球面距离、可视区域)if(DrawArgs.Camera.ViewRange

DrawArgs.Camera.Latitude, DrawArgs.Camera.Longitude)

)

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

publicvirtualvoidComputeChildren(DrawArgs drawArgs)

{

//判断是否是最高级别if(Level+1>=QuadTileSet.ImageStores[0].LevelCount)return;

//计算瓦片的中点经纬度doubleCenterLat=0.5f*(South+North);doubleCenterLon=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)

//Returns the QuadTile for specified location if available.///Tries to queue a download if not available.//Initialized QuadTile if available locally, else null.privateQuadTile ComputeChild(doublechildSouth,doublechildNorth,doublechildWest,doublechildEast)

{

QuadTile child=newQuadTile(

childSouth,

childNorth,

childWest,

childEast,this.Level+1,

QuadTileSet);returnchild;

}

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

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

//Builds flat or terrain mesh for current tile///publicvirtualvoidCreateTileMesh()

{

verticalExaggeration=World.Settings.VerticalExaggeration;

m_CurrentOpacity=QuadTileSet.Opacity;

renderStruts=QuadTileSet.RenderStruts;if(QuadTileSet.TerrainMapped&&Math.Abs(verticalExaggeration)>1e-3)

//创建具有高程值的格网(今天要关注的)CreateElevatedMesh();else

//创建没有高程值的平面格网CreateFlatMesh();

}

591行CreateElevatedMesh()

创建具有高程值的格网

//Build the elevated terrain mesh///protectedvirtualvoidCreateElevatedMesh()

{

isDownloadingTerrain=true;

//vertexCountElevated值为40；向四周分别扩充一个样本点//Get height data with one extra sample around the tiledoubledegreePerSample=LatitudeSpan/vertexCountElevated;

//获取具有高程值的TerrainTile对象(这是最关键部分，深入分析)

TerrainTile tile=QuadTileSet.World.TerrainAccessor.GetElevationArray(North+degreePerSample, South-degreePerSample, West-degreePerSample, East+degreePerSample, vertexCountElevated+3);float[,] heightData=tile.ElevationData;intvertexCountElevatedPlus3=vertexCountElevated/2+3;inttotalVertexCount=vertexCountElevatedPlus3*vertexCountElevatedPlus3;

northWestVertices=newCustomVertex.PositionNormalTextured[totalVertexCount];

southWestVertices=newCustomVertex.PositionNormalTextured[totalVertexCount];

northEastVertices=newCustomVertex.PositionNormalTextured[totalVertexCount];

southEastVertices=newCustomVertex.PositionNormalTextured[totalVertexCount];doublelayerRadius=(double)QuadTileSet.LayerRadius;//Calculate mesh base radius (bottom vertices)//Find minimum elevation to account for possible bathymetryfloatminimumElevation=float.MaxValue;floatmaximumElevation=float.MinValue;foreach(floatheightinheightData)

{if(height

minimumElevation=height;if(height>maximumElevation)

maximumElevation=height;

}

minimumElevation*=verticalExaggeration;

maximumElevation*=verticalExaggeration;if(minimumElevation>maximumElevation)

{//Compensate for negative vertical exaggerationminimumElevation=maximumElevation;

maximumElevation=minimumElevation;

}doubleoverlap=500*verticalExaggeration;//500m high tiles//Radius of mesh bottom gridmeshBaseRadius=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=newBoundingBox((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 meshesintvertexCountElevatedPlus2=vertexCountElevated/2+2;

vertexIndexes=newshort[2*vertexCountElevatedPlus2*vertexCountElevatedPlus2*3];intelevated_idx=0;for(inti=0; i

{for(intj=0; 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+1)*vertexCountElevatedPlus3+j);

vertexIndexes[elevated_idx++]=(short)((i+1)*vertexCountElevatedPlus3+j+1);

}

calculate_normals(refnorthWestVertices, vertexIndexes);

calculate_normals(refsouthWestVertices, vertexIndexes);

calculate_normals(refnorthEastVertices, vertexIndexes);

calculate_normals(refsouthEastVertices, vertexIndexes);

isDownloadingTerrain=false;

}

596行TerrainTile tile = QuadTileSet.World.TerrainAccessor.GetElevationArray(North + degreePerSample, South - degreePerSample, West - degreePerSample, East + degreePerSample, vertexCountElevated + 3);

获取样本点的高程值数组。

使用了TerrainAccessor.cs类120行代码

publicvirtualTerrainTile GetElevationArray(doublenorth,doublesouth,doublewest,doubleeast,intsamples)

{

TerrainTile res=null;

res=newTerrainTile(null);

res.North=north;

res.South=south;

res.West=west;

res.East=east;

res.SamplesPerTile=samples;

res.IsInitialized=true;

res.IsValid=true;doublelatrange=Math.Abs(north-south);doublelonrange=Math.Abs(east-west);float[,] data=newfloat[samples,samples];floatscaleFactor=(float)1/(samples-1);for(intx=0; x

{for(inty=0; y

{doublecurLat=north-scaleFactor*latrange*x;doublecurLon=west+scaleFactor*lonrange*y;

//关键，获取瓦片所有样本点的高程值data[x,y]=GetElevationAt(curLat, curLon,0);}

}

res.ElevationData=data;returnres;

}

关键代码：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)代码

privatestaticTerrainAccessor[] getTerrainAccessorsFromXPathNodeIterator(XPathNodeIterator iter,stringcacheDirectory)

{

System.Collections.ArrayList terrainAccessorList=newSystem.Collections.ArrayList();

//下面是读取DEM配置XML，并根据配置信息创建TerrainTileService对象和TerrainAccessor对象

while(iter.MoveNext())

{stringterrainAccessorName=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;

}doublenorth=0;doublesouth=0;doublewest=0;doubleeast=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)

{stringserverUrl=null;stringdataSetName=null;doublelevelZeroTileSizeDegrees=double.NaN;uintnumberLevels=0;uintsamplesPerTile=0;stringdataFormat=null;stringfileExtension=null;stringcompressionType=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=newTerrainTileService(

serverUrl,

dataSetName,

levelZeroTileSizeDegrees,

(int)samplesPerTile,

fileExtension,

(int)numberLevels,

Path.Combine(cacheDirectory, terrainAccessorName),

World.Settings.TerrainTileRetryInterval,

dataFormat);

TerrainAccessor newTerrainAccessor=newNltTerrainAccessor(

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{returnnull;

}

再来看看DEM的配置在哪里和XML内容吧

路径：

C:\Program Files\NASA\World Wind 1.4\Config\Earth.xml

配置内容：

?xml version="1.0" encoding="UTF-8"?> //黄色之间的XML就是一个TerrainAccessor配置http://worldwind25.arc.nasa.gov/wwelevation/wwelevation.aspxsrtm30pluszip20.012150Int16bilzip90.0-90.0-180.0180.0

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

获取特定点的高程值

//Get terrain elevation at specified location.//Latitude in decimal degrees.///Longitude in decimal degrees.//Returns 0 if the tile is not available on disk.publicoverridefloatGetElevationAt(doublelatitude,doublelongitude,doubletargetSamplesPerDegree)

{try{if(m_terrainTileService==null||targetSamplesPerDegree

{foreach(TerrainAccessor higherResSubinm_higherResolutionSubsets)

{if(latitude>higherResSub.South&&latitudehigherResSub.West&&longitude

{returnhigherResSub.GetElevationAt(latitude, longitude, targetSamplesPerDegree);

}

//自己可以看看如何完成TerrainTile的初始化构建的

TerrainTile tt=m_terrainTileService.GetTerrainTile(latitude, longitude, targetSamplesPerDegree);TerrainTileCacheEntry ttce=(TerrainTileCacheEntry)m_tileCache[tt.TerrainTileFilePath];if(ttce==null)

{

ttce=newTerrainTileCacheEntry(tt);

AddToCache(ttce);

}if(!ttce.TerrainTile.IsInitialized)

ttce.TerrainTile.Initialize();

ttce.LastAccess=DateTime.Now;

//获取高程值

returnttce.TerrainTile.GetElevationAt(latitude, longitude);}catch(Exception)

{

}return0;

}

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

publicfloatGetElevationAt(doublelatitude,doublelongitude)

{try{doubledeltaLat=North-latitude;doubledeltaLon=longitude-West;

//TileSizeDegrees为当前级别下瓦片的度数大小

//计算方法：158行tile.TileSizeDegrees = m_levelZeroTileSizeDegrees * Math.Pow(0.5, tile.TargetLevel);

//注意思考：SamplesPerTile-1为什么是减去1？传进来的SamplesPerTile=43而不是44？

//如果传入的是44，这里应该减2doubledf2=(SamplesPerTile-1)/TileSizeDegrees;floatlat_pixel=(float)(deltaLat*df2);floatlon_pixel=(float)(deltaLon*df2);

//这里是将点，近似成包含点的最小正方形(经纬度取整)intlat_min=(int)lat_pixel;intlat_max=(int)Math.Ceiling(lat_pixel);intlon_min=(int)lon_pixel;intlon_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;floatdelta=lat_pixel-lat_min;

//根据外矩形四顶点的经纬度分别插值计算中间线的高程floatwestElevation=ElevationData[lon_min, lat_min]*(1-delta)+ElevationData[lon_min, lat_max]*delta;floateastElevation=ElevationData[lon_max, lat_min]*(1-delta)+ElevationData[lon_max, lat_max]*delta;

//利用插值计算点的高程值

delta=lon_pixel-lon_min;floatinterpolatedElevation=westElevation*(1-delta)+eastElevation*delta;returninterpolatedElevation;

}catch{

}return0;

}

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

读取BIL文件高程数据

//This method initializes the terrain tile add switches to///Initialize floating point/int 16 tiles///publicvoidInitialize()

{if(IsInitialized)return;if(!File.Exists(TerrainTileFilePath))

{//Download elevationif(request==null)

{using( request=newTerrainDownloadRequest(this, m_owner, Row, Col, TargetLevel) )

{

request.SaveFilePath=TerrainTileFilePath;

request.DownloadInForeground();

}

}if(ElevationData==null)

ElevationData=newfloat[SamplesPerTile, SamplesPerTile];if(File.Exists(TerrainTileFilePath))

{//Load elevation filetry{//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=newFileInfo(TerrainTileFilePath);if(tileInfo.Length==0)

{

TimeSpan age=DateTime.Now.Subtract( tileInfo.LastWriteTime );if(age

{//This tile is still flagged badIsInitialized=true;

}else{//remove the empty 'flag' fileFile.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=newBinaryReader(s);if(m_owner.DataType=="Int16"){for(inty=0; y

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)

{thrownewApplicationException(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.thrownewApplicationException(String.Format("Error while trying to read terrain tile {0}", TerrainTileFilePath), ex);

}

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