[原创]常见DEM格式分析与三维分层设色显示

楼主^#

更多发布于：2007-11-12 11:51

在城市地质调查项目中需要读取CNSDTF DEM数据,这种格式与美国USGS DEM数据格式相比较,有些简单,不包含椭球体,投影方式等元数据,实际上就是投影后的平面格网数据,所以读取入系统不困难.以后将介绍ERDAS IMG格式的读取,这个就显得有些复杂了. std::locale::global(std::locale(""));
 ifstream _demfile; //DEM文件流 _demfile.open(lpszPathName); //读取文件 string sDataName,sCNDataName; string sVersion;float iVersion; string sAlpha;float fAlpha; string sUnit; string sUnitName; string sCompress; float dCompress; string sX0; double dX0; string sY0; double dY0; string sDX; double dDX; string sDY; double dDY; string sRow; int iRows; string sColumn; int iColumns; string sValueTypeName; string sValueType; string sHzoom; int iHzoom; string sMinV; string sMaxV;
 if(_demfile) { _demfile>>sDataName>>sCNDataName; _demfile>>sVersion>>iVersion; _demfile>>sAlpha>>fAlpha; _demfile>>sUnit>>sUnitName; _demfile>>sCompress>>dCompress; _demfile>>sX0>>dX0; _demfile>>sY0>>dY0; _demfile>>sDX>>dDX; _demfile>>sDY>>dDY; _demfile>>sRow>>iRows; _demfile>>sColumn>>iColumns; _demfile>>sValueTypeName>>sValueType; _demfile>>sHzoom>>iHzoom; _demfile>>sMinV; _demfile>>sMaxV;
 } else { AfxMessageBox("无法读取文件数据头！"); return false; } dY0 = dY0 - dDY*iRows; //从文件中读取的dYO记录的是最大数值 long index = 0,PointsTotalNumber = 0; int i,j; double tmpHeight = 0; POINT3d* p3dDEMGridPoint = new POINT3d[iRows*iColumns]; //准备读文件中数据 double dNodata = -99999; double dZMin = 9999,dZMax = -9999; for(j=iRows-1;j>=0;j--) //由于此文件记录Y数值的方法与其他不同，需要翻转数据 { for(i=0;i<iColumns;i++) { index = j*iColumns+i; _demfile>>tmpHeight; if(iHzoom != 0) tmpHeight = tmpHeight/iHzoom; //原始的缩放因子 p3dDEMGridPoint[index][2] = tmpHeight; if( tmpHeight > -9999) //也就是不等于 dNodata { if(dZMin > tmpHeight) dZMin = tmpHeight; //记录高程上的最小数值 if(dZMax < tmpHeight) dZMax = tmpHeight; //记录高程上的最大数值 } } }

 
<IMG src="http://www.cnblogs.com/images/cnblogs_com/wuhanhoutao/111912/r_layercoloring1.jpg" border=0>
美国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.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; } 
<IMG src="http://hiphotos.baidu.com/wuhanhoutao/pic/item/22cdc0ee53e9c0222cf534f8.jpg" border=0>
 SDTS(Spatial Data Transfer Standard ) DEM 读取也是颇费力的一件事情,首先这个体系并不是一个单一的数据文件,而是由一系列的相同后缀(.ddf)构成文件的组合.其中,*CATD.DDF是首要文件,相当于一个索引,告知其他文件的功能.剩下的文件中,*LDEF.DDF可以获得格网的长宽数值。*IREF.DDF文件可以得到XY轴的分辨率。*CELL.DDF是真正带有数据的文件。可以提取高程数值。另外还涉及到读取8211文件的问题,所以说,是件比较费力的事情. 在<a href="http://wuhanhoutao.cnblogs.com/" target="_blank" >http://wuhanhoutao.cnblogs.com/</A>提供程序框架和结果图供大家参考。愿与大家多多交流.

喜欢0 评分0

举报回复

发帖回复

« 返回列表

您需要登录后才可以回帖，登录或者注册

返回顶部

[原创]常见DEM格式分析与三维分层设色显示

最新喜欢：