读取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;
}