读取USGS DEM数据显示三维分层设色地图

转自 http://www.cnblogs.com/wuhanhoutao/archive/2007/10/31/944835.html
美国USGS的DEM数据读取的程序代码已经提供，关键是对文件头的处理要注意每个字段的含义，到数据部分需要注意XY轴的最大和最小数值，Z轴的最大最小数值提取的时候，要防止NodataValue数据的干扰，因为部分数据是NodataValue，会非常小或大。

读取美国USGS的DEM数据程序

string quadName,processInfo;
double seGeoCornerX,seGeoCornerY;
long processCode;
string sectionIndicator, mapCenterCode;
long levelCode,elevPattern,groundRefSysCode,groundRefSysZone,groundRefSysUnits;
long elevUnits, numPolySides;
DEMPointVector demCorners;
double counterclockAngle;
long elevAccuracyCode;

double spatialResX,spatialResY,spatialResZ;
long profileRows,profileColumns;
long largeContInt,maxSourceUnits,smallContInt,minSourceUnits;
long sourceDate, inspRevDate;
string inspFlag;
long valFlag;
long suspectVoidFlg;
long vertDatum,horizDatum;
long dataEdition;
long perctVoid;
long westEdgeFlag,northEdgeFlag,eastEdgeFlag,southEdgeFlag;
double vertDatumShift; //通过文件流得到DEM文件头信息

char bufstr[1024];
char temp[1024];
long i;

DEMUtil::getRecord(dem,bufstr); //从文件流中获得信息

strncpy(temp, bufstr, 40); //从获得信息中拷贝到临时变量
temp[40] = '\0';
quadName = temp;

strncpy(temp,bufstr+40,40); //从获得信息中拷贝到临时变量
temp[40] = '\0';
processInfo = temp;

DEMUtil::getDouble(bufstr, 109, 13, seGeoCornerX);
DEMUtil::getDouble(bufstr, 122, 13, seGeoCornerY);
processCode = DEMUtil::getLong(bufstr, 135, 1); //135 = 122 + 13

strncpy(temp,bufstr+137,3);
temp[3] = '\0';
sectionIndicator = temp;

strncpy(temp,bufstr+140,4);
temp[4] = '\0';
mapCenterCode = temp;

levelCode = DEMUtil::getLong(bufstr, 144, 6); //level编码
elevPattern = DEMUtil::getLong(bufstr, 150, 6);
groundRefSysCode = DEMUtil::getLong(bufstr, 156, 6);
groundRefSysZone = DEMUtil::getLong(bufstr, 162, 6);
groundRefSysUnits = DEMUtil::getLong(bufstr, 528, 6); //得到平面上使用的度量单位
elevUnits = DEMUtil::getLong(bufstr, 534, 6); //得到Z轴上的度量单位
numPolySides = DEMUtil::getLong(bufstr, 540, 6);

for (i = 0; i < 4; i++)
{
double x,y;
long pos = 546 + (i * 48);
DEMUtil::getDouble(bufstr, pos, 24, x);
DEMUtil::getDouble(bufstr, pos + 24, 24, y);
demCorners.push_back(DEMPoint(x,y)); //各角部顶点数据
}

DEMUtil::getDouble(bufstr, 738, 24, m_dminElevation); //得到最小高程数值
DEMUtil::getDouble(bufstr, 762, 24, m_dmaxElevation); //得到最大高程数值
DEMUtil::getDouble(bufstr, 786, 24, counterclockAngle );
elevAccuracyCode = DEMUtil::getLong(bufstr, 810, 6);
DEMUtil::getDouble(bufstr, 816, 12, spatialResX); //X轴分辨率
DEMUtil::getDouble(bufstr, 828, 12, spatialResY); //Y轴分辨率
DEMUtil::getDouble(bufstr, 840, 12, spatialResZ); //Z轴分辨率
profileRows = DEMUtil::getLong(bufstr, 852, 6); //得到行数
profileColumns = DEMUtil::getLong(bufstr, 858, 6); //得到列数
largeContInt = DEMUtil::getLong(bufstr, 864, 5);
maxSourceUnits = DEMUtil::getLong(bufstr, 869, 1);
smallContInt = DEMUtil::getLong(bufstr, 870, 5);
minSourceUnits = DEMUtil::getLong(bufstr, 875, 1);
sourceDate = DEMUtil::getLong(bufstr, 876, 4);
inspRevDate = DEMUtil::getLong(bufstr, 880, 4);

strncpy(temp, bufstr+884,1);
temp[1]='\0';
inspFlag = temp;

valFlag = DEMUtil::getLong(bufstr, 885, 1);
suspectVoidFlg = DEMUtil::getLong(bufstr, 886, 2);
vertDatum = DEMUtil::getLong(bufstr, 888, 2); //垂直
horizDatum = DEMUtil::getLong(bufstr, 890, 2); //水平
dataEdition = DEMUtil::getLong(bufstr, 892, 4);
perctVoid = DEMUtil::getLong(bufstr, 896, 4);
westEdgeFlag = DEMUtil::getLong(bufstr, 900, 2); //西部
northEdgeFlag = DEMUtil::getLong(bufstr, 902, 2); //北部
eastEdgeFlag = DEMUtil::getLong(bufstr, 904, 2); //东部
southEdgeFlag = DEMUtil::getLong(bufstr, 906, 2); //南部
DEMUtil::getDouble(bufstr, 908, 7, vertDatumShift);

m_p3dDEMGrid->SetXMin( demCorners[0].getX() );
m_p3dDEMGrid->SetYMin( demCorners[0].getY() );
m_p3dDEMGrid->SetXCellSize( spatialResX );
m_p3dDEMGrid->SetYCellSize( spatialResY );
m_p3dDEMGrid->SetColumns( profileColumns );

long retval;

retval = FillGeographic(dem,incrementalRead);

long ImportDEMDataUSGS::FillGeographic(istream& dem,bool incrementalRead)
{
if (m_bReadingFileHeadOnlyOnce)
{
m_lCurProfile = 0;
}

while (m_lCurProfile < m_p3dDEMGrid->GetColumns()) //断面数据还没有读完
{
m_vecProfiles.push_back(DEMProfile());
dem >> m_vecProfiles.back(); //存储断面数据
m_lCurProfile++;
if (incrementalRead) //倘若采用增量读法，则马上返回
return m_p3dDEMGrid->GetColumns() - m_lCurProfile + 1;
}

int width = m_p3dDEMGrid->GetColumns();
int height = m_vecProfiles[0].getNumberOfElevations(); //假设所有的断面和第一个断面有一样的一连串高程数值
m_p3dDEMGrid->SetRows( height );

m_ppt3dPoint = new POINT3d[width * height]; //added

double dXMin = m_p3dDEMGrid->GetXMin();
double dYMin = m_p3dDEMGrid->GetYMin();
double dXCellSize = m_p3dDEMGrid->GetXCellSize();
double dYCellSize = m_p3dDEMGrid->GetYCellSize();
int i,j;
long Index;
for (i = 1; i < width; i++)
{
DEMElevationVector const& elev = m_vecProfiles[i].getElevations(); //得到这个断面的一连串高程数值
for (j = 0; j < height; j++)
{
Index = width * j + i;
m_ppt3dPoint[Index][0]=dXMin + i*dXCellSize;
m_ppt3dPoint[Index][1]=dYMin + j*dYCellSize;
m_ppt3dPoint[Index][2]=elev[j]/m_dPlaneHeightRadio; //保证三轴单位统一
}
}

i=0; //为了填充数组第一列
DEMElevationVector const& elev_1 = m_vecProfiles[1].getElevations();
for(j = 0; j < elev_1.size();j++)
{
Index = width * j + i;
m_ppt3dPoint[Index][0]=dXMin + i*dXCellSize;
m_ppt3dPoint[Index][1]=dYMin + j*dYCellSize;
m_ppt3dPoint[Index][2]=elev_1[j]/m_dPlaneHeightRadio;

}
m_vecProfiles.erase(m_vecProfiles.begin(), m_vecProfiles.end());
return 0;
}