爱华网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(DrawArgsdrawArgs)
{
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,"Couldnotloadeffect"+m_effectPath+":"+errs);
Log.Write(Log.Levels.Warning,"Effecthasbeendisabled.");
m_effectPath=null;
m_effect=null;
}
}
if(ImageStores[0].LevelZeroTileSizeDegrees<180)
{
//Checkforlayeroutsideview
doublevrd=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(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(QuadTileqtinm_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<tileSpread;i++)
{
for(doublej=middleCenterLat-i*ImageStores[0].LevelZeroTileSizeDegrees;j<middleCenterLat+i*ImageStores[0].LevelZeroTileSizeDegrees;j+=ImageStores[0].LevelZeroTileSizeDegrees)
{
for(doublek=middleCenterLon-i*ImageStores[0].LevelZeroTileSizeDegrees;k<middleCenterLon+i*ImageStores[0].LevelZeroTileSizeDegrees;k+=ImageStores[0].LevelZeroTileSizeDegrees)
{
//根据经纬度和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，则更新QuadTile
QuadTileqt=(QuadTile)m_topmostTiles[key];
if(qt!=null)
{
qt.Update(drawArgs);
continue;
}
//Checkfortileoutsidelayerboundaries，获取外边框四点经度或纬度坐标
doublewest=-180.0f+curCol*ImageStores[0].LevelZeroTileSizeDegrees;
if(west>m_east)
continue;
doubleeast=west+ImageStores[0].LevelZeroTileSizeDegrees;
if(east<m_west)
continue;
doublesouth=-90.0f+curRow*ImageStores[0].LevelZeroTileSizeDegrees;
if(south>m_north)
continue;
doublenorth=south+ImageStores[0].LevelZeroTileSizeDegrees;
if(north<m_south)
continue;
//结合中不存在，创建新的QuadTile
qt=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(Exceptioncaught)
{
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(DrawArgsdrawArgs)
{
if(m_isResetingCache)
return;
try
{
doubletileSize=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) publicvirtualvoidComputeChildren(DrawArgsdrawArgs)
{
//判断是否是最高级别
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) ///<summary>
///ReturnstheQuadTileforspecifiedlocationifavailable.
///Triestoqueueadownloadifnotavailable.
///</summary>
///<returns>InitializedQuadTileifavailablelocally,elsenull.</returns>
privateQuadTileComputeChild(doublechildSouth,doublechildNorth,doublechildWest,doublechildEast)
{
QuadTilechild=newQuadTile(
childSouth,
childNorth,
childWest,
childEast,
this.Level+1,
QuadTileSet);
returnchild;
}
QuadTile.cs 中的CreateTileMesh()方法用来创建瓦片格网的，分别在Initialize() 、Update()方法中调用。
410行 这里调用的CreateElevatedMesh();是添加DEM数据创建高程格网的。
///<summary>
///Buildsflatorterrainmeshforcurrenttile
///</summary>
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()
创建具有高程值的格网 ///<summary>
///Buildtheelevatedterrainmesh
///</summary>
protectedvirtualvoidCreateElevatedMesh()
{
isDownloadingTerrain=true;
//vertexCountElevated值为40；向四周分别扩充一个样本点
//Getheightdatawithoneextrasamplearoundthetile
doubledegreePerSample=LatitudeSpan/vertexCountElevated;
//获取具有高程值的TerrainTile对象（这是最关键部分，深入分析）
TerrainTiletile=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;
//Calculatemeshbaseradius(bottomvertices)
//Findminimumelevationtoaccountforpossiblebathymetry
floatminimumElevation=float.MaxValue;
floatmaximumElevation=float.MinValue;
foreach(floatheightinheightData)
{
if(height<minimumElevation)
minimumElevation=height;
if(height>maximumElevation)
maximumElevation=height;
}
minimumElevation*=verticalExaggeration;
maximumElevation*=verticalExaggeration;
if(minimumElevation>maximumElevation)
{
//Compensatefornegativeverticalexaggeration
minimumElevation=maximumElevation;
maximumElevation=minimumElevation;
}
doubleoverlap=500*verticalExaggeration;//500mhightiles
//Radiusofmeshbottomgrid
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=newBoundingBox((float)South,(float)North,(float)West,(float)East,
(float)layerRadius,(float)layerRadius+10000*this.verticalExaggeration);
QuadTileSet.IsDownloadingElevation=false;
//Buildcommonsetofindexesforthe4childmeshes
intvertexCountElevatedPlus2=vertexCountElevated/2+2;
vertexIndexes=newshort[2*vertexCountElevatedPlus2*vertexCountElevatedPlus2*3];
intelevated_idx=0;
for(inti=0;i<vertexCountElevatedPlus2;i++)
{
for(intj=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(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行代码
publicvirtualTerrainTileGetElevationArray(doublenorth,doublesouth,doublewest,doubleeast,intsamples)
{
TerrainTileres=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<samples;x++)
{
for(inty=0;y<samples;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);
360docimg_10_
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行
360docimg_11_ getTerrainAccessorsFromXPathNodeIterator(XPathNodeIterator iter, string cacheDirectory)代码 privatestaticTerrainAccessor[]getTerrainAccessorsFromXPathNodeIterator(XPathNodeIteratoriter,stringcacheDirectory)
{
System.Collections.ArrayListterrainAccessorList=newSystem.Collections.ArrayList();
//下面是读取DEM配置XML，并根据配置信息创建TerrainTileService对象和TerrainAccessor对象
while(iter.MoveNext())
{
stringterrainAccessorName=iter.Current.GetAttribute("Name","");
if(terrainAccessorName==null)
{
//TODO:Throwexception?
continue;
}
XPathNodeIteratorlatLonBoxIter=iter.Current.Select("LatLonBoundingBox");
if(latLonBoxIter.Count!=1)
{
//TODO:Throwexception?
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));
XPathNodeIteratortileServiceIter=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对象
TerrainTileServicetts=newTerrainTileService(
serverUrl,
dataSetName,
levelZeroTileSizeDegrees,
(int)samplesPerTile,
fileExtension,
(int)numberLevels,
Path.Combine(cacheDirectory,terrainAccessorName),
World.Settings.TerrainTileRetryInterval,
dataFormat);
TerrainAccessornewTerrainAccessor=newNltTerrainAccessor(
terrainAccessorName,
west,
south,
east,
north,
tts,
higerResolutionSubsets);
terrainAccessorList.Add(newTerrainAccessor);
}
//TODO:AddFloatingpointterrainAccessorcode
//TODO:AddWMSAccessorcodeandmakeitworkinTerrainAccessor(whichitcurrentlydoesn't)
}
if(terrainAccessorList.Count>0)
{
return(TerrainAccessor[])terrainAccessorList.ToArray(typeof(TerrainAccessor));
}
else
{
returnnull;
}
}
再来看看DEM的配置在哪里和XML内容吧
路径：
C:Program FilesNASAWorld Wind 1.4ConfigEarth.xml
配置内容：
?xmlversion="1.0"encoding="UTF-8"?>
<WorldName="Earth"EquatorialRadius="6378137.0"LayerDirectory="Earth"xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"xsi:noNamespaceSchemaLocation="WorldXmlDescriptor.xsd">
<TerrainAccessorName="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)方法。
360docimg_13_ 获取特定点的高程值 ///<summary>
///Getterrainelevationatspecifiedlocation.
///</summary>
///<paramname="latitude">Latitudeindecimaldegrees.</param>
///<paramname="longitude">Longitudeindecimaldegrees.</param>
///<paramname="targetSamplesPerDegree"></param>
///<returns>Returns0ifthetileisnotavailableondisk.</returns>
publicoverridefloatGetElevationAt(doublelatitude,doublelongitude,doubletargetSamplesPerDegree)
{
try
{
if(m_terrainTileService==null||targetSamplesPerDegree<World.Settings.MinSamplesPerDegree)
return0;
if(m_higherResolutionSubsets!=null)
{
foreach(TerrainAccessorhigherResSubinm_higherResolutionSubsets)
{
if(latitude>higherResSub.South&&latitude<higherResSub.North&&
longitude>higherResSub.West&&longitude<higherResSub.East)
{
returnhigherResSub.GetElevationAt(latitude,longitude,targetSamplesPerDegree);
}
}
}
//自己可以看看如何完成TerrainTile的初始化构建的
TerrainTilett=m_terrainTileService.GetTerrainTile(latitude,longitude,targetSamplesPerDegree);
TerrainTileCacheEntryttce=(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，这里应该减2
doubledf2=(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）文件来获取的。
360docimg_15_ 读取BIL文件高程数据 ///<summary>
///Thismethodinitializestheterraintileaddswitchesto
///Initializefloatingpoint/int16tiles
///</summary>
publicvoidInitialize()
{
if(IsInitialized)
return;
if(!File.Exists(TerrainTileFilePath))
{
//Downloadelevation
if(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))
{
//Loadelevationfile
try
{
//TerrainDownloadRequest'sFlagBadTile()createsemptyfiles
//asawaytoflag"bad"terraintiles.
//Removetheempty'flag'filesafterpresettime.
try
{
FileInfotileInfo=newFileInfo(TerrainTileFilePath);
if(tileInfo.Length==0)
{
TimeSpanage=DateTime.Now.Subtract(tileInfo.LastWriteTime);
if(age<m_owner.TerrainTileRetryInterval)
{
//Thistileisstillflaggedbad
IsInitialized=true;
}
else
{
//removetheempty'flag'file
File.Delete(TerrainTileFilePath);
}
return;
}
}
catch
{
//Ignoreanyerrorsintheaboveblock,andcontinue.
//Forexample,ifsomeonehadtheempty'flag'file
//open,thedeletewouldfail.
}
//读取BIL文件数据的关键代码，可以被我们借鉴使用
using(Streams=File.OpenRead(TerrainTileFilePath))
{
BinaryReaderreader=newBinaryReader(s);
if(m_owner.DataType=="Int16")
{
for(inty=0;y<SamplesPerTile;y++)
for(intx=0;x<SamplesPerTile;x++)
ElevationData[x,y]=reader.ReadInt16();
}
if(m_owner.DataType=="Float32")
{
for(inty=0;y<SamplesPerTile;y++)
for(intx=0;x<SamplesPerTile;x++)
{
ElevationData[x,y]=reader.ReadSingle();
}
}
IsInitialized=true;
IsValid=true;
}
return;
}
catch(IOException)
{
//IfthereisanIOexceptionwhenreadingtheterraintile,
//theneithersomethingiswrongwiththefile,orwith
//accesstothefile,sotryandremoveit.
try
{
File.Delete(TerrainTileFilePath);
}
catch(Exceptionex)
{
thrownewApplicationException(String.Format("Errorwhiletryingtodeletecorruptterraintile{0}",TerrainTileFilePath),ex);
}
}
catch(Exceptionex)
{
//Someothertypeoferrorwhenreadingtheterraintile.
thrownewApplicationException(String.Format("Errorwhiletryingtoreadterraintile{0}",TerrainTileFilePath),ex);
}
}
}
另注：SRTM的高程数据存放路径如下图：（DEM跟影像也是分级存放的，存放方式一致）
360docimg_17_
至此，如何加载DEM数据创建网格的过程已经分析完了。
接下来，继续分析QuadTile.cs中CreateElevatedMesh()和Render方法，内容主要是DirectX编程，稍后添加……
Tag标签: World Wind
无痕客
关注 - 6
粉丝 - 14
关注博主
2
0
0 (请您对文章做出评价)
« 上一篇：WorldWind学习系列十三：地形数据（DEM）加载和应用（入门篇）
» 下一篇：WorldWind学习系列十五：如何切割影像和DEM数据及其在WW中的应用配置
