DEM

数字高程模型（DEM）内插程序

数字高程模型（DEM）内插程序报告

姓名：

学号： 18103617

专业：测绘工程

日期： 2012年12月

编写数字高程模型（DEM ）内插程序

一. 实验目的

掌握反距离加权（IDW ）数字高程模型内插原理及其内插子程序的设计方法，了解其它逐点高程内插方法的基本原理。 二.实验 内容

根据测量学实习所测的南师大仙林校区北区高程采样点（X i ,Y i ,Z i ）,将区域划分为5m ×5m 的规则格网，以每个待求格网中心顶点为中心，按最多不超过8个点、最少不小于5个点确定搜索圆半径，采用反距离加权（IDW ）法依据落入搜索圆内的各高程采样点的高程求取各待求网格中心点的高程。最后输出该区域的规则格网DEM 。 三.数据准备

将该区域的高程采样点坐标按以下格式准备（由于数据量大，只显示部分数据）：

点号 x y z

1 397144.429 3555459.261 38.310

2 397178.172 3555494.435 38.171

3 397157.91

4 3555472.37

5 38.70

6 4 397151.716 3555465.459 38.706 5 397149.153 3555463.114 38.609 6 397149.059 3555462.875 38.531

7 397146.69

8 3555460.996 38.454 8 397166.911 3555482.147 38.430 …… ……

……

(1270)

397312.564 3555189.716

22.236

四. 操作步骤 （1）读入已知点的坐标，根据所采集的采样点计算插值区域; （2）按5m ×5m 分辨率将以上插值区域划分为规则格网； （3）按行、列顺利依次对各格网顶点确定搜索圆（保证落入圆内的采样点数目为5-8

个）；

（4）依据落入搜索圆内的采样点高程Zi ，按反距离加权内插法计算格网顶点的高程，

并予以记录待定点Z ： ∑∑==?

=

m

i

m

i i i d

d z z 1

i 21

2

1

1

；

（5）返回（3），直到所有格网单元均已内插。

五. 程序运行界面以及结果

读入数据的文件格式：

输出结果

六. 程序代码

#include "stdafx.h"

#include"stdio.h"

#include"stdlib.h"

#define Zis 10 //数据点个数

#define us 9 //使用数据点个数const int n=6;

double D[Zis][4];

double M[us][6];//系数阵

double M1[6][us];//系数阵转置

double Zi[us][1];//高程值

double P[us][us];//权值矩阵

double M1P[6][us];

double M1PM[us][6];

double M1PM_1M1[6][us];

double M1PM_1M1P[6][us];

double X[6][1];

double z;

void zhaodian(double x,double y)

{

double temp[4];

//计算出数据点与定点距离平方

for(int i=0;i

D[i][3]=(D[i][0]-x)*(D[i][0]-x)+(D[i][1]-y)*(D[i][1]-y);

//按距离平方大小对数据进行从小到大排序

for(int i=0;i

{

for(int j=i;j

{

if(D[i][3]>D[j][3])

{

for(int k=0;k<4;++k)

{

temp[k]=D[i][k];

D[i][k]=D[j][k];

D[j][k]=temp[k];

}

}

}

}

for(int i=0;i

{

for(int j=0;j

{

if(i==j)

P[i][j]=D[i][3];

else

P[i][j]=0;

}

}

//使用距离定点最近的8个数据值计算系数阵

for(int i=0;i

{

M[i][0]=(D[i][0]-x)*(D[i][0]-x);

M[i][1]=(D[i][0]-x)*(D[i][1]-y);

M[i][2]=(D[i][1]-y)*(D[i][1]-y);

M[i][3]=D[i][0]-x;

M[i][4]=D[i][1]-y;

M[i][5]=1;

Zi[i][0]=D[i][2];

}

}

template void ReadD(T *Matrix,int m,int n)

{

ifstream readin("D.txt");

while(!readin.eof )

{

for(int i=0;i

for(int j=0;j

readin>>Matrix[i][j];

}

}

template void xiangcheng(T1 *M1,T2 *M2,T2 *Re,int a,int b,int c)//矩阵相乘

{

int i,j,k;

for(i=0;i

{for(j=0;j

{

Re[i][j]=0;

for(k=0;k

Re[i][j]+=M1[i][k]*M2[k][j];

}

}

return;

}

templatevoid zhihuan(T1 *M1,T2 *M2,int m,int n)//矩阵置换

{ int i,j;

for(i=0;i

for(j=0;j

M2[j][i]=M1[i][j];

return;

}

void inverse(double c[n][n])//矩阵求逆

{ int i,j,h,k;

double p;

double q[n][12];

for(i=0;i

for(j=0;j

q[i][j]=c[i][j];

for(i=0;i

for(j=n;j<12;j++)

{if(i+6==j)

q[i][j]=1;

else

q[i][j]=0;}

for(h=k=0;k

{if(q[i][h]==0)

continue;

p=q[k][h]/q[i][h];

for(j=0;j<12;j++)

{ q[i][j]*=p;

q[i][j]-=q[k][j];

}

}

for(h=k=n-1;k>0;k--,h--) // 消去对角线以上的数据for(i=k-1;i>=0;i--)

{if(q[i][h]==0)

continue;

p=q[k][h]/q[i][h];

for(j=0;j<12;j++)

{q[i][j]*=p;

q[i][j]-=q[k][j];}}

for(i=0;i

{ p=1.0/q[i][i];

for(j=0;j<12;j++)

q[i][j]*=p;}

for(i=0;i

for(j=0;j

c[i][j]=q[i][j+6];

}

int _tmain(int argc, _TCHAR* argv[])

{

FILE *fp;

if((fp=fopen("d.txt","r"))==NULL)

printf("cannot find the file\n");

for(int i=0;i

for(int j=0;j<3;j++)

fscanf(fp,"%lf",&D[i][j]);

fclose(fp);

zhaodian(x,y);

zhihuan(M,M1,8,6);

xiangcheng(M1,P,M1P,6,8,8);

xiangcheng(M1P,M,M1PM,6,8,6);

inverse(M1PM);

xiangcheng(M1PM,M1,M1PM_1M1,6,6,8);

xiangcheng(M1PM_1M1,P,M1PM_1M1P,6,8,8);

xiangcheng(M1PM_1M1P,Zi,X,6,8,1);

z=X[5][0];

printf("（%lf,%lf)高程值为：%lf\n",x,y,z);

printf("系数阵M为\n");

for(int i=0;i

{

for(int j=0;j<6;++j)

printf("%6.0lf\t",M[i][j]);

printf("\n");

}

printf("权阵为\n");

for(int i=0;i

{

for(int j=0;j

printf("%4.0lf\t",P[i][j]);

printf("\n");

}

printf("\n求得X矩阵为\n");

for(int i=0;i<6;++i)

printf("%8.4lf\n",X[i][0]);

return 0;

}