由DEM格网数据生成等高线

楼主^#

更多发布于：2008-01-15 23:25

现阶段DEM数据的来源不再仅仅局限于原来保存的纸质地形图的跟踪数字化，可能是采用激光测距仪配合地面控制点生成，由此产生的DEM格网数据可能需要生成等高线图，以作为基础地理数据或底图。此工艺流程与原来跟踪等高线来生产DEM的目标刚好相反，是在拥有DEM数据后反而希望生成等高线。下面介绍具体的流程和参考代码：
1) 针对格网中每个基本单元，创建一个参考变量，即YuanSu[j](0<=i<=Row; 0<=j<=Column)同时分配内存。此参考变量赋值0或1，其中1代表此位置处的邻域高程数值将发生变化，那么此位置将绘制等高线。
col = (char **)calloc(m_iRows, sizeof(char *)); //指向指针的指针初始化
 row = (char **)calloc(m_iRows, sizeof(char *)); //指向指针的指针初始化
 for(y=0; y<m_iRows; y++) //m_iRows:格网的行数
 {
 col[y] = (char *)calloc(m_iColumns, sizeof(char)); //每个单元分配内存，m_iColumns : 格网的列数
 row[y] = (char *)calloc(m_iColumns, sizeof(char)); //每个单元分配内存
 }
2）在以高程数值从最小值开始依次递增的最外部一层循开始，接着对整个格网数据进行循环，在给参考变量YuanSu[j]赋数值0或1之后，然后以YuanSu[j] == 1为判断条件，寻找并记录等高线。
 for(zValue=zMin, ID=0; zValue<=zMax; zValue+=zStep) //从最小Z数值计算起，知道最大Z数值
 {
 //-------------------------------------------------
 // Step1: Find Border Cells //找到边界单元，边界：高度趋势发生改变的地方。
 for(y=0; y<m_iRows-1; y++)
 for(x=0; x<m_iColumns-1; x++)
 if( m_p3dDEMPoints[y*m_iColumns+x][2] >= zValue )
 {
 row[y][x] = m_p3dDEMPoints[y*m_iColumns+(x+1)][2] < zValue ? 1 : 0; //‘’就代表边界
 col[y][x] = m_p3dDEMPoints[(y+1)*m_iColumns+x][2] < zValue ? 1 : 0;
 }
 else
 {
 row[y][x] = m_p3dDEMPoints[y*m_iColumns+(x+1)][2] >= zValue ? 1 : 0; //‘’就代表边界
 col[y][x] = m_p3dDEMPoints[(y+1)*m_iColumns+x][2] >= zValue ? 1 : 0;
 }
 //-------------------------------------------------
 // Step2: Interpolation + Delineation
 for(y=0; y<m_iRows-1; y++)
 for(x=0; x<m_iColumns-1; x++)
 {
 if(row[y][x]) //如果数值为
 { 
 for(i=0; i<2; i++) //两次调用函数，ID值增加
 {
 if(i==0) m_iCountourfind = 0;
 else m_iCountourfind = 1;
 ContourFind(x, y, zValue, true, ID++);
 }
 row[y][x] = 0; //代表已经处理过
 
 }
 if(col[y][x]) //如果数值为
 { 
 for(i=0; i<2; i++) //两次调用函数，ID值继续增加
 {
 if(i==0) m_iCountourfind = 0;
 else m_iCountourfind = 1;
 ContourFind(x, y, zValue, false, ID++);
 }
 col[y][x] = 0; //代表已经处理过
 }
 }
 }
3）寻找等高线的过程，是在本位置的高程数值，本条等高线的高程ZValue和邻近位置的高程数值做插值处理后，得到X、Y位置并做记录。
 do
 {
 //-------------------------------------------------
 // Interpolation...
 d = m_p3dDEMPoints[y*m_iColumns+x][2]; //取本位置的高程数值
 d = (d - z) / (d - m_p3dDEMPoints[zy*m_iColumns+zx][2]); //以本位置的高程数值，ZValue和邻近位置的高程数值做插值处理
 xPos = m_dXMin + m_iXCellSize * (x + d * (zx - x)); //通过高度差的比例，线性插值求得X、Y位置
 yPos = m_dYMin + m_iYCellSize * (y + d * (zy - y)); //
 
 if(m_iCountourfind == 0)
 {
 if (m_iPointNum_one >= m_iBufferSize_one) 
 {
 m_iBufferSize_one += 1024;
 m_point3dtmp_one = (POINT3d *)realloc(m_point3dtmp_one,m_iBufferSize_one * sizeof(POINT3d)); //为顶点增加新的内存分配 
 }
 if(m_point3dtmp_one!=NULL)
 {
 m_point3dtmp_one[m_iPointNum_one][0] = xPos; //记录下X位置
 m_point3dtmp_one[m_iPointNum_one][1] = yPos; //记录下Y位置
 m_point3dtmp_one[m_iPointNum_one][2] = z+0.2; //此处高度为ZValue,增加.2个单位的目的是为了显示时抬高一些，便于观察。
 m_iPointNum_one++;
 }
 }
 else
 {
 if (m_iPointNum_two >= m_iBufferSize_two) 
 {
 m_iBufferSize_two += 1024;
 m_point3dtmp_two = (POINT3d *)realloc(m_point3dtmp_two,m_iBufferSize_two * sizeof(POINT3d)); //为顶点增加新的内存分配 
 }
 if(m_point3dtmp_two!=NULL)
 {
 m_point3dtmp_two[m_iPointNum_two][0] = xPos; //记录下X位置
 m_point3dtmp_two[m_iPointNum_two][1] = yPos; //记录下Y位置
 m_point3dtmp_two[m_iPointNum_two][2] = z+0.2; //此处高度为ZValue,增加.2个单位的目的是为了显示时抬高一些，便于观察。
 m_iPointNum_two++;
 }
 }
 
 //-------------------------------------------------
 // Naechstes (x/y) (col/row) finden...
 if( !ContourFindNext(Dir, x ,y ,doRow) ) //x、y将发生增或减的变化，doRow也改变数值
 doContinue = ContourFindNext(Dir, x, y, doRow);
 Dir = (Dir + 5) % 8;
 //-------------------------------------------------
 // Loeschen und initialisieren...
 if(doRow)
 {
 row[y][x] = 0; //表示已经处理过
 zx = x + 1; //重新计算数值
 zy = y;
 }
 else
 {
 col[y][x] = 0; //表示已经处理过
 zx = x; //重新计算数值
 zy = y + 1;
 }
 }
 while(doContinue); 
4）寻找等高线过程结束的条件：
bool doContinue;
 if(doRow)
 {
 switch(Dir) //邻域上个方向的判断
 {
 case 0:// Norden
 if(row[y+1][x ])
 { y++; Dir=0; doContinue=true; break; }
 case 1:// Nord-Ost
 if(col[y ][x+1])
 { x++; doRow=false;Dir=1; doContinue=true; break; }
 case 2:// Osten ist nicht...
 case 3:// Sued-Ost
 if(y-1>=0) if(col[y-1][x+1])
 { x++; y--; doRow=false;Dir=3; doContinue=true; break; }
 case 4:// Sueden
 if(y-1>=0) if(row[y-1][x ])
 { y--; Dir=4; doContinue=true; break; }
 case 5:// Sued-West
 if(y-1>=0) if(col[y-1][x ])
 { y--; doRow=false;Dir=5; doContinue=true; break; }
 case 6:// Westen ist nicht...
 case 7:// Nord-West
 if(col[y ][x ])
 { doRow=false;Dir=7; doContinue=true; break; }
 default:
 Dir=0; doContinue=false;
 }
 }
 else
 {
 switch(Dir) //邻域上个方向的判断
 {
 case 0:// Norden ist nicht...
 case 1:// Nord-Ost
 if(row[y+1][x ])
 { y++; doRow=true; Dir=1; doContinue=true; break; }
 case 2:// Osten
 if(col[y ][x+1])
 { x++; Dir=2; doContinue=true; break; }
 case 3:// Sued-Ost
 if(row[y ][x ])
 { doRow=true; Dir=3; doContinue=true; break; }
 case 4:// Sueden ist nicht...
 case 5:// Sued-West
 if(x-1>=0) if(row[y ][x-1])
 { x--; doRow=true; Dir=5; doContinue=true; break; }
 case 6:// Westen
 if(x-1>=0) if(col[y ][x-1])
 { x--; Dir=6; doContinue=true; break; }
 case 7:// Nord-West
 if(x-1>=0) if(row[y+1][x-1])
 { x--; y++; doRow=true; Dir=7; doContinue=true; break; }
 default:
 Dir=0; doContinue=false;
 }
 }
 return(doContinue); 原始DEM显示图: <IMG src="http://www.cnblogs.com/images/cnblogs_com/wuhanhoutao/120009/r_DEM_Original.jpg" border=0> 生成的等高线图: <IMG src="http://www.cnblogs.com/images/cnblogs_com/wuhanhoutao/120009/r_To_Countour.jpg" border=0>

喜欢0 评分0

A friend is never known till a man has need. ...CL

举报回复

发帖回复

« 返回列表

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

返回顶部

由DEM格网数据生成等高线

最新喜欢：