读取美国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();
    doub