gis
gis
管理员
管理员
  • 注册日期2003-07-16
  • 发帖数15945
  • QQ554730525
  • 铜币25337枚
  • 威望15352点
  • 贡献值0点
  • 银元0个
  • GIS帝国居民
  • 帝国沙发管家
  • GIS帝国明星
  • GIS帝国铁杆
阅读:5362回复:0

等值线模型的实现

楼主#
更多 发布于:2015-05-26 10:11
等值线通俗易懂、应用广泛,是一种重要的地学模型,在GIS 中有不可忽视的地位。本文以实例作为切入口,用桌面软件、OLE 自动化、COM 二次开发和底层算法编程四种不同的方法生成等值线,并试图以此强化对理论知识的掌握。
还是先了解一些最基本的概念。
1. 数字地形模型(DTM):是地形表面形态属性信息的数字表达,是带有空间位置特征和地形属性特征的数字描述。
2. 数字高程模型(DEM):数字地形模型中地形属性为高程时称为数字高程模型。
3. 规则格网(Grid):DEM 的一种表达方式。将区域空间切分为规则的格网单元,每个格网单元对应一个数值。
4. 不规则三角网(TIN):DEM 的一种表达方式。根据区域有限个点集将区域划分为相连的三角面网络,区域中任意点落在三角面的顶点、边上或三角形内。
5. 等值线(Contour):DEM 的一种表达方式。是数值相同的点依次连接而成的曲线。
规则格网
不规则三角网

从等值线的定义不难看出,要生成等值线有两个步骤:首先找到数值相同的点,然后将它们依次连接。不过在生产实际中,我们通常不是按照数值相同的点来采集数据的。这一点很好理解,举个例子,大家都知道秦岭-淮河是我国800mm 等降水量线,可是我们不可能沿着秦岭的山脊线和淮河的中心线采集数据,又如在野外测量测的一般都是地面特征点,如山脊、山谷和陡坎等,这些点的高程肯定是不相等的。因此,要找到等值点,我们只能根据我们采集到的离散点数据进行推算,这一从已知点推算未知点的过程称为插值。此外,由于通过测量或监测得到的离散点数据通常都是以某种格式的文本文件存储,因此,我们需要把文字记录的点变成有地理位置的图形,这一把文本文件中的记录展绘为图形容器中的点的过程称为展点。另一个问题是,高程相等的点有无数个,怎么可能全找出来呢?因此,为了能够找出有限个等值点,使它们能够表达出该区域的地形,我们必须借助另外两种模型:规则格网(Grid)或不规则三角网(TIN)。
综上所述,生成等值线的过程大致为:读取离散点数据文件->展点->建立规则格网(Grid)或不规则三角网(TIN)->插值找到等值点并跟踪绘制等值线。
本文使用的数据文件Data.txt 中数据格式如下:
1 0.3 6.1 870
2 1.4 6.2 793
……

四列数值分别表示编号、X 坐标、Y 坐标和高程Z。这里需要注意的是,坐标是平面坐标,而不是球面坐标(经纬度)。
一、利用桌面软件:ArcMap
ArcGIS 的桌面软件中提供了3D Analyst 模块,可以建立三维模型并进行分析。第一步我们需要读取文本文件并展点,由于ArcMap 只能从数据库中读取字段信息,所以我们必须先将数据导入数据库,再导入图层。在导入图层时需要指定X 字段和Y 字段,并设置坐标系统(这里用投影坐标而不用地理坐标),这样系统才能将一条条文字记录转化为一个个具有地理坐标的点要素。要素类生成后,只要调出3D Analyst 模块生成TIN,再生成等值线就可以了。具体操作步骤如下:
1.启动Access,新建数据库data.mdb,在空白处右击,导入data.txt,四个字段分别为Code(设为主健)、X、Y、Z,保存为表Point。
2.启动ArcMap,Tools->Add XY Data,在“Choose a table”中添加数据库data.mdb 中的表Point,在“X Field”和“Y Field”中分别选择字段X、Y,Edit->Select,选择/Projected CoordinateSystems/World/Mercator (world).prj,一路OK 后,离散点已被展绘到Point Events 图层中。
3.Tools->Extensions,把3D Analyst 打上勾,Close。右击工具栏,把3D Analyst 打上勾。
4.3D Analyst->Create/Modify TIN->Create TIN From Features,把图层Point Events 打上勾,在“Height Source”中选择字段Z,在“Output TIN”中输入TIN 的保存路径,OK 后生成了TIN。
5.3D Analyst->Surface Analyst->Contour,在“Contour Interval”中输入等高距10,在“OutputTIN”中输入Contour 的保存路径,OK 后便绘制好了等值线。

ArcGIS 素以功能强大著称,但同时也背上了“繁杂”、“不友好”等恶名,不过这里很轻松地就生成了等值线,可见DEM 的建立并不是一件复杂的事。补充一句,在我们的国产软件GeoStar 的GeoTIN 模块中,等值线的生成更加简便快捷,但考虑到有条件接触该软件的学生不多,这里还是选用了ArcMap 演示。
二、利用OLE 自动化:VB+Surfer 8
在桌面软件中绘制等值线非常容易,那么想在自己的开发的系统中实现这一功能是否困难呢?首先介绍一下用Surfer 软件自动化编程。
Surfer 是Golden Soft 生产的一款非常著名的地形制图软件,我们GIS 专业的学生可能知道的不多,测绘、土木、水利等专业的学生应该都接触过。这里用它来演示,主要因为与GIS 软件不同,它是通过格网模型来生成等高线的。Surfer 的桌面软件界面极其亲切,功能强大且操作方便,深受用户喜爱。从Surfer7 开始,该产品使用了OLE 技术,使用户能够方便地在自己的程序中添加它的强大功能。以VB 环境为例,我们只需要在工程->引用中,把“Surfer 8 Type Library”打上勾,便可以在程序中用少量语句调用它的全部功能——注意这里只是调用,Surfer 的界面虽然不出现,但实际上它是在后台运行的。
顺便补充一下,很多软件都支持自动化,使用CreateObject 函数便可创建并返回对一个应用程序对象的引用。如Set xlApp = CreateObject("Excel.Application"),之后就可以在Object 类型的变量xlApp的基础上调用Excel 的功能。
在这里我们适用语句Set objSurfer = CreateObject("Surfer.Application")创建了一个Surfer对象,然后用它的函数对子对象进行了一系列操作。
全部代码如下:
Dim objSurfer, objPlot, objMapFrame As Object
Private Sub cmdRun_Click()
Dim strInFile, strGridFile, strOutFile As String
strInFile = App.path & "\data.txt"
strGridFile = Left(strInFile, Val(InStrRev(strInFile, ".")) - 1) + ".grd" '把扩展名改为.grd
Set objSurfer = CreateObject("Surfer.Application") '创建Surfer 对象
objSurfer.Visible = False 'surfer 软件本身不在前台显示
objSurfer.GridData DataFile:=strInFile, xCol:=2, yCol:=3, zCol:=4, Algorithm:=srfKriging,
OutGrid:=strGridFile '把离散的数据文件转换为格网文件。几个参数分别表示离散点数据文件、X 所在列、Y 所在列、Z 所在列、插值法、输出的格网文件名
Set objPlot = objSurfer.Documents.Add(srfDocPlot) '创建srf 文档
Set objMapFrame = objPlot.Shapes.AddContourMap(strGridFile) '创建等值线图并添加到文档
objPlot.Export FileName:=Left(strInFile, Val(InStrRev(strInFile, ".")) - 1) + ".bmp" '输出bmp
Picture1.Picture = LoadPicture(Left(strInFile, Val(InStrRev(strInFile, ".")) - 1) + ".bmp")'显示bmp
End Sub
不难看出,这段程序用GridData 函数建立了规则格网,然后在这个规则格网的基础上用AddContourMap 函数生成了等值线。其中,在建立规则格网时使用了克里金(Kriging)插值法。

可以看出,用Surfer 绘制出的等值线图明显比ArcMap 的光滑。然而,OLE 技术不能脱离软件本身使用,也就是说,用这种方法制作的系统只能在安装了Surfer 的计算机上使用,这显然使得程序的通用性大大降低。更重要的是,这只是一张栅格图片,不能进行任何操作,而且Surfer 本身也没有任何DEM 分析功能,所以如果仅仅是制图的话,可以考虑使用它,如果要进行GIS 应用,就得另觅他法了。
三、利用COM 二次开发:VB+Supermap Objects 5
使用组件GIS 无疑是最常用的方法。在众多的组件产品中,MapX 和MapObject 都没有DEM 相关接口,而ArcObject 又过于复杂,SuperMap Objects 以其强大的功能和简单的结构赢得了广泛的赞誉。利用它的SuperAnalyst 控件,我们能够非常方便地生成TIN,利用它的Super3D 控件,我们可以十分快捷地绘制等值线。
首先我们需要把原始的离散点数据加工一下,制作成点数据集。尽管这也能编程实现,但是事先用桌面软件SuperMap Deskpro 处理一下显然更好。
1.创建数据源test.mdb,将test.tab 导入,生成点数据集dataP(导入数据集不支持文本文件,可以用工具->从文本文件创建SDB,但是对文本文件的格式有固定要求,这一点不如GeoTIN方便)
2.双击dataP,保存地图,保存工作空间test.smw,退出Deskpro
3.在test.smw 相同目录下创建VB 工程,在窗体上添加SuperWorkspace 控件、SuperMap 控件、SuperAnalyst 控件各一个,将Super3D 控件加载到工具箱中(即在部件中打上勾)。
4.代码如下:
Private Sub Form_Load()
Dim objModel As soModelingOperator '三维表面建模操作对象
Dim objDatasource As soDataSource '数据源对象
Dim objPointDt As soDatasetVector '矢量数据集对象
Dim obj3dAnalyst As New so3DAnalyst '三维分析对象
SuperWorkspace1.Open App.Path & "\test.smw" '打开工作空间
SuperMap1.Connect SuperWorkspace1.Object '连接地图控件
SuperAnalyst1.Connect SuperWorkspace1.Handle '连接空间分析控件
SuperMap1.OpenMap SuperWorkspace1.Maps.Item(1) '打开地图
SuperMap1.ViewBounds = SuperMap1.Bounds '全图显示
Set objDatasource = SuperWorkspace1.Datasources(1) '获取数据源
Set objPointDt = SuperMap1.Layers(1).Dataset '获取点数据集
Set objModel = SuperAnalyst1.SurfaceAnalyst.Modeling '获取建模分析对象
objModel.CreateTin objPointDt, "Z", , objDatasource, "TIN" '离散点数据生成TIN
obj3dAnalyst.TINToContour objDatasource.Datasets("TIN"), objDatasource, "Contour", 10,
10, True '生成等高线
SuperMap1.Layers.AddDataset SuperWorkspace1.Datasources(1).Datasets("Contour"), ture'添加数据集End Sub
与Surfer 不同的是,SuperMap 生成的不再是位图,而是一个线数据集,每条等高线都是一个独立的对象,具有SmLength(长度)和dZvalue(高程)等属性。添加如下代码,可以选中并查询等高线的高程。
Private Sub Command1_Click()
SuperMap1.Action = scaSelect '选择状态
End Sub
Private Sub SuperMap1_DblClick()
Dim objRecordset As soRecordset '记录集对象
Set objRecordset = SuperMap1.selection.ToRecordset(False) '将选择集转化为记录集
MsgBox objRecordset.GetFieldValue("dZvalue") '显示该记录dZvalue 字段的值
End Sub
遗憾的是,在生成的等值线图中,右下角有一处明显的两条等值线相交的错误。为此,我在SuperMap Deskpro 里也试了一下,该错误依然存在。并且,在图的左侧、右侧和上方各有一个离散点在等值线上成为“尖点”。不过,这两个问题在ArcMap 和GeoTIN 中同样
存在,但在Surfer 中则不存在。看来GIS 软件在这方面与专业的制图软件还有一定的差距。

除了SuperMap Objects 和ArcObject,还有一款著名的软件提供了具有强大DEM 分析功能的组件,它就是Matlab 的组件MatrixVB。小巧玲珑的MatrixVB 可以使用Matlab 的很多“令人垂涎”的函数,却能够完全独立于庞大无比的Matlab 桌面软件。由于它本来就以矩阵运算闻名,因此自然是用规则格网生成等值线,用Contour 函数不仅可以生成2D 的等值线,还能生成2.5D 和3D 的,这里限于篇幅,不再赘述。
四、使用算法从底层编程
要深刻理解等值线模型实现的全过程,还是需要研究一下算法,最好能够编程实现。这里主要按照华东师范大学张超教授主编的《地理信息系统实习教程》(高等教育出版社)中的算法,参考了我校林君健、苍桂华老师主编的《摄影测量学》(国防工业出版社)中的一些思路,用VB 编写了“离散点数据->规则格网(Grid)->等值线”的程序,将等值线绘制在picturebox 上。
首先要读取离散点数据,这显然应该用数组来存储,每条记录对应一个数组元素,那么这个数组用什么数据类型呢?考虑到每个点都有Code,X,Y,Z,所以需要自定义类型。Public Type pointdata '点数据类型
code As Integer
x As Single
y As Single
Z As Single
End Type

这里定义了一个pointdata 类型,让我们再定义一个该类型的数组。Public point_data(1000) As pointdata '点数据
下面就可以读取离散点数据了,用Do 循环一行行地,每读一行都将用制表位隔开的四个数值分别存入数组的四个字段,并用一个整形变量total_num 来统计点数。
'读取数据并存入数组point_data(code,x,y,z)
total_num = 0
Do While Not EOF(1)
total_num = total_num + 1
Input #1, point_data(total_num).code, point_data(total_num).x, point_data(total_num).y,
point_data(total_num).Z '点号,X,Y,Z
Loop
好了,离散点数据读取完毕,准备建立规则格网模型,不过在这之前需要考虑一下坐标系统的问题。我们知道,picturebox 采用平面坐标系,而我们的离散点数据的X 和Y 都是地理坐标,所以必须先进行坐标系统的转换,否则将无法绘图。我们知道在二次开发时,地图控件都会提供坐标系统转换的接口,那么picturebox 会不会也有呢?按下F1 键,在MSDN中查“picturebox”,果然发现有个scale 方法,用以定义坐标系统。
根据MSDN 上的说明,我们需要确定坐上角和右下角的坐标。由于我们的离散点数据的坐标原点在左下方,故不难判断坐上角坐标为(xmin, ymax),右下角坐标为(xmax, ymin)。于是问题就成了“找出数组中的最大值和最小值”,这个大家应该都会吧?
xmax = point_data(1).x
xmin = xmax
ymax = point_data(1).y
ymin = ymax
For i = 1 To total_num
If point_data(i).x > xmax Then xmax = point_data(i).x
If point_data(i).x < xmin Then xmin = point_data(i).x
If point_data(i).y > ymax Then ymax = point_data(i).y
If point_data(i).y < ymin Then ymin = point_data(i).y
Next i
'更改坐标系
Picture1.Scale (xmin, ymax)-(xmax, ymin)
下面就要进入主题了。前面说过,规则格网的每个网格点上都有一个高程值,因此,我们必须利用离散点数据对所有格网点进行插值。DEM 内插有多种算法,这里采取最简单的一种——距离加权法。
《地》中的算法描述如下:
“设平面上分布一系列离散点,已知其坐标高程为xi,yi,zi(i=1,2,…,n),P(x,y)为任一格网点,根据周围离散点的高程,通过距离加权插值求P 点高程。周围点与P 点因分布位置的差异,对P(z)影响不同,我们把这种影响称为权函数Wi。权函数主要与距离有关,有时也与方向有关,若是在P 点周围四个方向上均匀取点,那么可不考虑方向因素,这时:

实践证明,Wi = 1 / di2 是较优的选择,di 为离散点至P 点的距离。”
分析如下:
1.步骤:公式中含有Σ,那么肯定需要用一个循环进行累加,循环中的每一轮处理一个离散点,对每一个离散点,需要计算它到待插点的距离的平方,然后求倒数得权(特别注意分子为零的情况要单独讨论),并与它的高程相乘。
2.参数:既然是对待插点处理,参数中肯定要有待插点的坐标,而离散点的数组当然也是必不可少的。
3.变量:将公式中的每个未知数都设为变量,此外为了方便,将分子和分母也分别作为变量。
程序如下:
'对一个待插点进行空间插值。采用距离加权插值法,与所有离散点比较。参数为待插点的横坐标、纵坐标,离散点高程值
Function nds(a As Double, b As Double, point_data() As pointdata) As Single
Dim fenzi As Double, fenmu As Double '分子,分母
Dim disdis(1 To 1000) As Double '距离的平方
Dim W(1 To 1000) As Double '权
Dim i As Integer
fenzi = 0
fenmu = 0
For i = 1 To total_num
disdis(i) = (a - point_data(i).x) * (a - point_data(i).x) + (b - point_data(i).y) * (b -
point_data(i).y) '待插点到某离散点的距离的平方
If disdis(i) = 0 Then '如果格网点与离散点重合
nds = point_data(i).Z
Exit Function
End If
W(i) = 1 / disdis(i)
fenzi = fenzi + W(i) * point_data(i).Z
fenmu = fenmu + W(i)
Next i
nds = fenzi / fenmu
End Function
调用上面的函数,我们就能够根据已知点求出任意未知点了。我们只需要用该函数对每个格网点进行计算,便可建立规则格网模型。
格网的数据结构非常简单,就是一个二维数组,第一维为行,第二维为列,值为格网点高程,其平面坐标可以通过行和列计算出来,不必单独存储。不过这里还有个注意事项,那就是地理坐标的原点在左下方,那么二维数组的原点该设在左下角还是坐上角呢?还有横向的x 该对应i 还是j 呢?还有从0 开始好,还是从1 开始好呢?或许按习惯吧。这里我选择从坐上角开始,i 为行,j 为列,都从0 计,与两本教材上的都不一样。

为了使程序更加灵活,用户应可以自己决定格网的规模,因此最大行数和最大列数分别用变量M 和N 表示。
'建立规则格网。
Public Sub build_grid()
Dim i As Integer, j As Integer
'重定义动态数组
ReDim grid(0 To M, 0 To N)
For i = 0 To M
For j = 0 To N
grid(i, j) = nds(xmin + j * dx, ymin + (M - i) * dy, point_data())
Next j
Next i
End Sub
这样规则格网也建好了,下面就可以做等值线了。前面说过,要生成等值线有两个步骤:首先找到数值相同的点,然后将它们依次连接。此方法就是内插计算全部等高线穿越格网边交点的坐标x,y,然后按等高线的顺序将属于每一条等高线的点找出来,并按等高线的走向将他们顺序连接成线。
在格网边上内插等高线点的坐标,可以用高项多项式法,这里使用精度不高的线性内插。

从图中可以看出,如果值为W 的一条等高线穿过网格,那么A 点两边的格网点高程肯定是一个大于W,一个小于W,B 点也一样。即等值点高程与两边网格点之差异号。只要满足该条件,我们就可以计算出等值点的位置。否则,做上记号表明值为W 的等值线不经过该边。
对于在垂直边上的A 点,其坐标为

对于在水平边上的B 点,其坐标为

由于水平边和垂直边需要区别对待,因此,这里根据《摄》中的算法,再定义两个数组,分别Hor(M,N)和Ver(M,N),分别用来存储水平边和垂直边。这样,很容易就能找出所有的等值点。
Cmax = Int(Hmax / Cstep) * Cstep '如102.1 对应为100
Cmin = Int(Hmin / Cstep + 1) * Cstep '如10.21 对应为20,因为值为10 的等值线不存在
For W = Cmin To Cmax Step Cstep
For j = 0 To N '检查垂直边
For i = 0 To M - 1 '对于垂直边,从左到右有N 个,从上到下只有M-1 个
If (grid(i, j) - W) * (grid(i + 1, j) - W) < 0 Then '存在等值点
Ver(i, j) = (W - grid(i, j)) / (grid(i + 1, j) - grid(i, j)) '求出等值点相
对位置
Else
Ver(i, j) = -1 '无等值点
End If
Next i
Next j
For i = 0 To M '检查水平边
For j = 0 To N - 1 '对于水平边,从上到下有M 个,从左到右只有N-1 个
If (grid(i, j) - W) * (grid(i, j + 1) - W) < 0 Then '存在等值点
Hor(i, j) = (W - grid(i, j)) / (grid(i, j + 1) - grid(i, j)) '求出等值点相对位置
Else
Hor(i, j) = -1 '无等值点
End If
Next j
Next i
Next W
然而,这里还有个“致命”问题,假如等值线刚好从网格点上穿过,那该把它算在哪条边上呢?所以为了避免这种情况的发生,在这段代码之前,需要找出所有与网格点相等的高程点(称为退化点)预处理一下,给他们加上一个非常小的修正值,这个微小改变对等值线的影响完全可以忽略不计,但是却简单漂亮地解决了上面的问题。
'检查退化点并消除,即如果等值线穿过格网点,则加上一个小值,使其穿过边For W = Cmin To Cmax Step Cstep
For i = 0 To M
For j = 0 To N
If grid(i, j) = W Then grid(i, j) = grid(i, j) + 0.00001
Next j
Next i
Next W
现在已经找到了全部等值点,下面要完成最后也是最艰难的任务——等值点的追踪。该过程非常复杂,由于时间紧张,我尚未完成正确代码,绘制出的等值线图也与上面用软件生成的图有较大差异,因此这里不便继续,有兴趣的读者可到GIS 大学论坛上下载源码自行研习。

小结:
很多GIS 初学者都对3D 很感兴趣,那么不妨自己好好研究一下DEM,这里只是对DEM中等值线模型的建立做了一个粗浅的研究性学习,希望能够起到抛砖引玉的效果。
喜欢0 评分0
游客

返回顶部