一种用于地理信息系统中的空间叠加分析方法和系统.pdf

本发明公开了一种用于GIS中的空间叠加分析方法和系统。该方法，包括下列步骤：将输入GIS中的图层的矢量数据转换成栅格数据，并采用游程编码表示；对该采用游程编码表示的栅格数据执行叠加操作；将叠加后的栅格数据再转换成矢量数据，得到经过叠加的图层。其利用了栅格算法的优势，能够避免计算几何算法的缺点，提高叠加效率。

1.  用于地理信息系统中的空间叠加分析方法，其特征在于，包括下列步骤：
步骤100.将输入地理信息系统中的图层的矢量数据转换成栅格数据，并采用游程编码表示；
步骤200.对该采用游程编码表示的栅格数据执行叠加操作；
步骤300.将叠加后的栅格数据再转换成矢量数据，得到经过叠加的图层。

2.  根据权利要求1所述的用于地理信息系统中的空间叠加分析方法，其特征在于，所述步骤100，包括下列步骤：
步骤110.扫描输入地理信息系统的一个图层中的每一个多边形，对每一个多边形再循环读入该多边形的每一个环，对于每个环，计算y值扫描线与每个弧段P₁P₂，设P₁.y＜P₂.y的断点I的横坐标I.x=(y-P2.y)×(P2.x-P1.x)(P2.y-P1.y)+P2.x,]]>以构造断点坐标加入断点栅格场；
步骤120.将所述多边形的所有环的断点加入到栅格场，建立经过所述y值扫描线扫描得到的行的游程集合并组织成链表形式。

3.  根据权利要求1所述的用于地理信息系统中的空间叠加分析方法，其特征在于，所述步骤200中，所述叠加操作，是对采用游程编码表示的栅格数据进行交、并、差布尔运算。

4.  根据权利要求1所述的用于地理信息系统中的空间叠加分析方法，其特征在于，所述步骤300，包括下列步骤：
步骤310.游程栅格边界矢量化，得到图层弧段；
步骤320.根据所述图层弧段，生成图层。

5.  根据权利要求4所述的用于地理信息系统中的空间叠加分析方法，其特征在于，所述步骤310，包括下列步骤：
步骤311.对每个游程单元，用两个数组表示各个游程单元的左边界、右边界是否被使用过，其初始值都是没有使用过；
步骤312.从游程链表中选择符合预设条件的游程单元，将其左边界或右边界作为追踪的初始点；
步骤313.从初始点开始进行矢量追踪，在追踪过程中，将经历过的游程左边界或右边界的相应标志值设置为使用过；
步骤314.重复步骤312-213，直至不存在任何相应标志值设置为使用过的游程，没有新的追踪初始点时终止循环。

6.  用于地理信息系统中的空间叠加分析系统，其特征在于，所述系统，包括：
矢转栅模块，用于将输入地理信息系统中的图层的矢量数据转换成栅格数据，并采用游程编码表示；
叠加模块，用于对该采用游程编码表示的栅格数据执行叠加操作；
叠加图层输出模块，用于将叠加后的栅格数据再转换成矢量数据，得到经过叠加的图层。

7.  根据权利要求6所述的用于地理信息系统中的空间叠加分析系统，其特征在于，所述矢转栅模块，包括：
扫面模块，用于扫描输入地理信息系统的一个图层中的每一个多边形，对每一个多边形再循环读入该多边形的每一个环，对于每个环，计算y值扫描线与每个弧段的断点的横坐标，以构造断点坐标加入断点栅格场；
游程编码模块，用于将所述多边形的所有环的断点加入到栅格场，建立经过所述y值扫描线扫描得到的行的游程集合并组织成链表形式。

8.  根据权利要求6所述的用于地理信息系统中的空间叠加分析系统，其特征在于，所述叠加模块，包括：
合并插入模块，用于对两个有重叠的游程单元间合并，否则把当前游程单元插入到适当的位置，并且在插入的过程中把游程单元的属性信息保存下来；
删除模块，用于在当前行游程链表内删除与待删除游程单元区域重叠的部分。

9.  根据权利要求6所述的用于地理信息系统中的空间叠加分析系统，其特征在于，所述叠加图层输出模块，包括：
边界追踪模块，用于追踪经过叠加后的用游程编码表示的栅格数据，得到图层弧段；
拓扑重构模块，用于根据所述图层弧段，生成图层。

10.  根据权利要求9所述的用于地理信息系统中的空间叠加分析系统，其特征在于，所述边界追踪模块，采用二值栅格的窗口化追踪方法。

一种用于地理信息系统中的空间叠加分析方法和系统
技术领域
本发明涉及地理信息系统(Geographic Information System，GIS)领域，特别是涉及一种用于地理信息系统中的空间叠加分析方法和系统。
背景技术
地理信息系统(Geographic Information System或Geo-Informationsystem，GIS)有时又称为“地学信息系统”或“资源与环境信息系统”。它是一种特定的十分重要的空间信息系统。它是在计算机硬、软件系统支持下，对整个或部分地球表层(包括大气层)空间中的多种地理空间实体数据及其关系，包括空间定位数据、图形数据、遥感图像数据、属性数据等，进行采集、储存、管理、运算、分析、显示和描述的技术系统。用于分析和处理在一定地理区域内分布的各种现象和过程，解决复杂的规划、决策和管理问题。叠加分析是地理信息系统空间分析中的核心部分，在地理信息系统空间分析中具有很重要的地位。
叠加分析是指在统一空间参考系统下，通过对两个数据进行的一系列集合运算，产生新数据的过程。这里提到的数据可以是图层对应的数据集，也可以是地物对象。叠加分析的目的就是通过对空间数据的加工或分析，提取用户需要的新的空间几何信息。在叠加分析中至少涉及到三个数据，其中一个数据的类型可以是点、线、面等，被称作输入数据；另一个数据是面数据被称作叠加数据；还有一个数据就是叠加结果数据，包含叠加后数据的几何信息和属性信息。
就目前公开发表的地理信息系统中有关叠加分析的文献而言，叠加分析算法普遍为基于矢量的方法。然而，基于矢量方法是面向物体的描述，一般地说物体间的关系尤其是几何关系的分析比较困难，相关的定位检索比较耗时，算法也较为复杂。例如，利用矢量方法建立拓扑关系，往往要进行频繁的点位匹配、坐标排序、角度计算等操作，对算法的效率有影响；对一些特殊情况，如建立弧-多边形拓扑关系时，可能出现含“岛”或具有环状结构的多边形，还可能有一些死胡同或孤立弧，利用矢量方法处理这类情况较为复杂。GIS中与矢量方法同等重要的是栅格方法。多数情况下栅格数据结构要比矢量数据结构简单得多，基于空间位置的分析也相对容易一些。因此，利用栅格方法进行GIS中的空间叠加分析，应该是GIS研究与实践中一项有意义的尝试。
发明内容
本发明的目的在于提供一种用于地理信息系统中的空间叠加分析方法和系统。是一种新的地图/多边形叠加算法，其利用了栅格算法的优势，能够避免计算几何算法的缺点，提高叠加效率。
为实现本发明的目的而提供的一种用于地理信息系统中的空间叠加分析方法，包括下列步骤：
步骤100.将输入地理信息系统中的图层的矢量数据转换成栅格数据，并采用游程编码表示；
步骤200.对该采用游程编码表示的栅格数据执行叠加操作；
步骤300.将叠加后的栅格数据再转换成矢量数据，得到经过叠加的图层。
所述步骤100，包括下列步骤：
步骤110.扫描输入地理信息系统的一个图层中的每一个多边形，对每一个多边形再循环读入该多边形的每一个环，对于每个环，计算y值扫描线与每个弧段P₁P₂，设P₁.y＜P₂.y的断点I的横坐标I.x=(y-P2.y)×(P2.x-P1.x)(P2.y-P1.y)+P2.x,]]>以构造断点坐标加入断点栅格场；
步骤120.将所述多边形的所有环的断点加入到栅格场，建立经过所述y值扫描线扫描得到的行的游程集合并组织成链表形式。
所述步骤200中，所述叠加操作，是对采用游程编码表示的栅格数据进行交、并、差布尔运算。
所述步骤300，包括下列步骤：
步骤310.游程栅格边界矢量化，得到图层弧段；
步骤320.根据所述图层弧段，生成图层。
所述步骤310，包括下列步骤：
步骤311.对每个游程单元，用两个数组表示各个游程单元的左边界、右边界是否被使用过，其初始值都是没有使用过；
步骤312.从游程链表中选择符合预设条件的游程单元，将其左边界或右边界作为追踪的初始点；
步骤313.从初始点开始进行矢量追踪，在追踪过程中，将经历过的游程左边界或右边界的相应标志值设置为使用过；
步骤314.重复步骤312-213，直至不存在任何相应标志值设置为使用过的游程，没有新的追踪初始点时终止循环。
为实现本发明的目的，还提供一种用于地理信息系统中的空间叠加分析系统，包括：
矢转栅模块，用于将输入地理信息系统中的图层的矢量数据转换成栅格数据，并采用游程编码表示；
叠加模块，用于对该采用游程编码表示的栅格数据执行叠加操作；
叠加图层输出模块，用于将叠加后的栅格数据再转换成矢量数据，得到经过叠加的图层。
所述矢转栅模块，包括：
扫面模块，用于扫描输入地理信息系统的一个图层中的每一个多边形，对每一个多边形再循环读入该多边形的每一个环，对于每个环，计算y值扫描线与每个弧段的断点的横坐标，以构造断点坐标加入断点栅格场；
游程编码模块，用于将所述多边形的所有环的断点加入到栅格场，建立经过所述y值扫描线扫描得到的行的游程集合并组织成链表形式。
所述叠加模块，包括：
合并插入模块，用于对两个有重叠的游程单元间合并，否则把当前游程单元插入到适当的位置，并且在插入的过程中把游程单元的属性信息保存下来；
删除模块，用于在当前行游程链表内删除与待删除游程单元区域重叠的部分。
所述叠加图层输出模块，包括：
边界追踪模块，用于追踪经过叠加后的用游程编码表示的栅格数据，得到图层弧段；
拓扑重构模块，用于根据所述图层弧段，生成图层。
所述边界追踪模块，采用二值栅格的窗口化追踪方法。
本发明的有益效果是：
1.本发明提出了一种用于GIS中的空间叠加分析算法，此方法利用栅格算法的优势，能避免计算几何算法的缺点，跟业界领先软件-ArcGIS相比，此算法有较大的优势和好的性价比；
2.本发明提出的基于栅格的叠加算法也同时可以通过根据预定义的地图比例尺事先生成的地图瓦片文件的方法来为网络地图服务器上提供空间分析功能；
3.本发明提出的基于栅格的叠加算法，在预生成栅格化表达的前提下，本算法的叠加效率很高，这为网络地图服务模式下的空间分析功能性能提升提供了很好的探索方向。
附图说明
图1是本发明空间叠加分析方法的流程图；
图2是利用本发明对输入的两个图层进行叠加分析的实施例的流程示意图；
图3是本发明中矢量数据转换成栅格数据的方法流程图；
图4是本发明中游程链表的叠加操作的结果示意图；
图5是本发明中栅格数据再转换成矢量数据得到经过叠加的图层的方法的流程图；
图6是本发明中边界矢量化的方法流程图；
图7是本发明中基于游程编码的游程栅格边界追踪图；
图8是本发明中基于栅格追踪的几种情况的示意图
图9是本发明中基于窗口的栅格追踪的走向图；
图10是本发明中一个能覆盖所有组合的基于窗口的追踪实例图；
图11是本发明中邻接多边形的边界追踪示意图；
图12是本发明中两个输入图层的示意图；
图13是对图12中的两个图层求交、差、并的结果示意图；
图14是本发明的空间叠加分析系统的结构示意图；
图15是本发明中一图层数据示意图；
图16是本发明的方法与ArcGIS的性能测试对比图；
图17是本发明方法三个关键步骤运行时间比例图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明的一种用于地理信息系统中的空间叠加分析方法和系统进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本发明的一种用于地理信息系统(Geographic Information System，GIS)中的空间叠加分析方法和系统，是通过将矢量数据转换成栅格数据，并对栅格数据进行叠加操作，再将经过叠加操作的栅格数据转换为矢量数据的方法实现GIS中基于栅格的空间叠加分析。
本发明要解决的技术问题包括：
1.基于扫描线和游程编码的矢量数据栅格化；
2.基于游程链表的图层叠加操作；
3.基于游程链表的栅格数据矢量化。
下面结合上述目标详细介绍本发明用于GIS中的空间叠加分析方法，图1是本发明空间叠加分析方法的流程图，图2是利用本发明对输入的两个图层进行叠加分析的实施例的流程示意图，如图1和图2所示，所述空间叠加分析方法，包括下列步骤：
步骤100.将输入GIS中的图层的矢量数据转换成栅格数据，并采用游程编码表示；
矢量数据和栅格数据是GIS中最重要的两种空间数据，它们之间的相互转换是GIS中的基础问题并得到广泛研究，矢量转栅格是其中的一个方面。在点、线、面三种基本类型要素中，面的处理相对复杂一些。目前已提出的针对面要素的矢栅转换方法有内部点扩散法、射线算法、扫描算法、边界代数算法等，它们采用规则格网系统，并以直接栅格编码方式(行列式)组织数据。直接栅格编码结构简单，操作方便，但存在数据量大、精度有限等缺点，采用压缩栅格编码是克服这些缺点的有效途径。
GIS中较常用的压缩栅格编码有链编码、游程编码、哈夫曼编码、四叉树编码、LZW编码等，它们通常是针对不同特征的栅格数据而设计的。对于GIS中的面状要素，一种简单而有效的压缩方法是游程编码，它可达到较高的压缩效率，且数据结构简单，易于操作。在游程编码环境中，很多空间计算过程以游程为操作单元，与以网格为单元的运算相比，这有助于提高相关算法的效率。因此，本发明中采用基于行扫描的游程编码表示栅格数据，并采用有序单向链表结构方便游程的基本操作并提高效率。
图3是本发明中矢量数据转换成栅格数据的方法流程图，如图3所示，所述步骤100，包括下列步骤：
步骤110.扫描输入GIS的一个图层中的每一个多边形，对每一个多边形再循环读入该多边形的每一个环(外环和内环)，对于每个环，计算y值扫描线与每个弧段P₁P₂(设P₁.y＜P₂.y)的交点I(断点)的横坐标(I.x=(y-P2.y)×(P2.x-P1.x)(P2.y-P1.y)+P2.x),]]>以构造断点坐标加入断点栅格场IP_GRID(每行皆为断点链表，断点结构包含坐标值val＝I.x和指向下一个断点的指针next)；
步骤120.将所述多边形的所有环的断点加入到栅格场，建立经过所述y值扫描线扫描得到的行的游程集合并组织成链表形式。
当把多边形的所有环的断点都加入到栅格场IP_GRID后(由于加入过程同时进行了排序，所以每行的断点是按序排列的)，依据点-面包含关系判定的射线法原理得证，每个游程单元的边界点由一对相邻断点组成，其中以奇数序号断点为起始点的游程属于面，所以可根据断点序号的奇偶性，快速建立该行的游程集合并组织成链表形式。每段游程RL_Linktable包含：游程的左值leftval和右值rightval、所属多边形编号rlid、指向下一段游程的指针next。
较佳地，若某游程rl满足条件：rl.rightval-rl.leftval≤dxy，即表示游程单元的边界点的两个相邻断点的值分别取整后在一个格子里，这种情况认为小于栅格精度，所以则不将此游程加入到链表中，否则会使矢量化追踪走向错误从而导致追踪标记错误，产生死循环或者多余内环等。这种做法的额外好处是能简化数据，在精度要求不是很高但速度要求比较严格的应用场合，此方法可提供数据简化功能。
较佳地，本发明中采用实数值记录游程边界点的方法，可以更有效提高游程栅格化的位置精度。这点也正是用直接栅格编码所不能实现的，因为在直接栅格编码中，栅格的位置是通过数组下标隐式表示的，不能通过改变数据类型来提高精度；而游程编码完全可以用更高精度的实数型数据进行记录。
较佳地，本发明中还为输入图层中的多边形设置了要素ID信息，并将其保存在链表结构，这样做的好处之一是矢量化追踪的时候可以根据多边形的不同属性进行区分，之二是在输出结果图层中也能保存其属性信息，以便以后做叠加分析的Web版本。
所述步骤100的操作，可以通过下述伪代码实现：
For(每个多边形)
 {
   For(每个环)
   {
      For(每个弧段)
      {
        令该弧段两个端点为P1，P2，且P1.y＜＝P2.y
        For(y＝P1.y；y＜P2.y)；y++)
        {
          计算Y与P1P2的交点横坐标，以此作为断点
          在栅格场中插入断点(X)
          }
   }
}//一个多边形的各个环都处理完毕
 //根据各行的断点，插入总场
 for(i；i＜栅格场总行数；i++)
 {
   if(该行断点个数＞＝2)
   {
       每次选取2个断点，形成游程链表区域
       if(执行的是裁剪操作)
          在栅格场中移除游程(XLEFT，XRIGHT)
      else
      在栅格场中插入游程(XLEFT，XRIGHT)
   }
   }
}//这个图层的所有多边形都做完了
步骤200.对该采用游程编码表示的栅格数据执行叠加操作；
由于步骤100中已经将输入的图层在栅格场中表示，并采用游程编码表示，因此本步骤中对栅格数据执行叠加操作，即转变为对游程链表中的游程单元的叠加操作。游程链表的基本操作是基于某一行的游程而言的。
图4是本发明中游程链表的叠加操作的结果示意图，如图4所示，对游程链表中的游程单元的叠加操作，分为以下情况：
合并以及插入。两个游程单元间有重叠就合并，否则把当前游程单元插入到适当的位置，并且在插入的过程中把游程单元的属性信息保存下来。
在已有行游程基础上加入一个新的游程，合并游程重叠的部分，可以看作是一个普通有序链表的插入过程的演化。不同之处在于现在的游程的值有2个域值(相当于数轴上的一个区间)，所以插入过程就是进行区间的叠置：2个游程间有重叠就合并，否则把当前游程插入到适当的位置。并且在插入的过程中把游程的属性信息也保留下来，如游程所属的多边形的id号等，以便矢量化的使用。
游程单元的删除操作。在当前行游程链表内删除与待删除游程单元区域重叠的部分。
在一个有序的游程链表中，删除待删除游程所指示的游程空间。这个过程正好和插入过程相反。
较佳地，基于步骤100中采用游程编码表示栅格数据，本发明中提供一种采用基于游程链表的交、并、差等布尔操作算法实现叠加操作，如图4所示。在基于游程链表表示的图层栅格数据进行叠加时，其本质是对两个已排序的区间(游程)链表进行集合操作，其复杂度为2n。如果栅格场内共有M行游程链表，每行链表最多有接近n段游程，本算法此步实现的复杂性可表示为2×M×(n-3)，即复杂度为O(n)。
作为一种可实施方式，假设输入A、B两个图层，某行的游程链表如图4所示，标号为游程属于的多边形ID号。图4中列出了求交、并、差的结果。求交时，把两段游程单元的相交的部分作为新的游程单元加入结果集的过程中同时把多边形的id号也当作属性信息(id0，id1)添加进去，这样在结果中就能知道是这两个图层的哪两个多边形的交集了。例如第一段结果集是由A的第0号多边形和B的第0号多边形求交得出的，所以其id0＝0，id1＝0。求差时比较简单，把减数的游程删除，结果集只需要保留被减图层的ID号信息即可。求并时，是这里面比较复杂的，不能一味的合并，而要把各部分有区别的加入到结果集，比如属于A不属于B的部分，在加入结果集的同时，需要把这段游程的id1置为-1，这样的好处是求并的结果集和求交的结果集数据结构都一样，可以共享矢量化过程。
步骤300.将叠加后的栅格数据再转换成矢量数据，得到经过叠加的图层。
将表达面域的栅格数据向矢量数据转换，是提取具有相同属性编码的栅格集合的矢量边界及边界与边界之间拓扑关系的过程，一般包含有多边形边界提取、边界线搜索、拓扑关系生成、去除冗余点并进行曲线圆滑等步骤，核心问题是边界追踪和拓扑重构。文献“基于二值栅格数据的边界矢量化和面信息生成的一种实现方法”采用了基于窗口匹配技术的边界追踪算法和提出了基于闭合边界曲线嵌套关系判定的方法，但无法区分相邻多边形的情况。本发明在边界追踪和拓扑重构两个核心过程中提出了基于要素属性信息的矢量化过程，能够巧妙的区分邻接多边形的矢量化。
本步骤中，将所述叠加后的栅格数据再转换成矢量数据的方法，包括边界矢量化(即边界追踪)和面的生成(即拓扑重构)两个部分。
图5是本发明中栅格数据再转换成矢量数据得到经过叠加的图层的方法的流程图，如图5所示，所述步骤300，包括下列步骤：
步骤310.游程栅格边界矢量化；
本发明中游程栅格和普通矩阵式二值栅格的区别主要在于栅格图形的表现形式有所不同。显然，二值栅格在表达栅格图形的过程中比游程链表直观，形象。然而，游程栅格的矢量化过程与二值栅格的矢量化方法大体相同，不同的仅仅是本处所处理的栅格场结构不同而已；另一方面，游程栅格边界矢量追踪过程中只需要处理每一个游程端点(游程起点和终点)，而大量栅格图中的内部点不需要判别和追踪。因此，矢量化过程中仍然采用二值栅格的窗口化追踪方法，只是具体到每个游程的端点时需要做一个转化：将游程端点处的值转化为二值栅格中的行列值，以便直接进行追踪。
图6是本发明中边界矢量化的方法流程图，如图6所示，所述步骤310，包括下列步骤：
步骤311.对每个游程单元，用lsign，rsign两个数组表示各个游程单元的左边界、右边界是否被使用过(追踪过)，其初始值都是FALSE，即没有使用过；
步骤312.从游程链表中选择符合条件的游程单元(左边界或右边界没有被使用过)，将其左边界或右边界作为追踪的初始点；
图7是本发明中基于游程编码的游程栅格边界追踪图，如图7所示，如果是从左边界开始追踪，表明发现了一个新面(从L1点即第一个游程的左边界开始追踪，其结果是一个面的外边界)；如果从右边界开始追踪，如图4中从R1点追踪，其结果将是一个面的内边界。
在实际追踪过程中，总是从行号最小最左边的游程的左边界开始追踪，那么追踪出的第一个外边界一定是一个面的外边界；当我们从某个游程的右边界追踪出一个多边形时，它一定是一个内边界(因为该游程的左边界已经被追踪过)，该内边界与该游程左边界追踪出的多边形属于同一个多边形，这就将所有的内边界与外边界建立了关联关系。
步骤313.从初始点开始进行矢量追踪。即根据当前点关联的4个栅格的数值来判断走向，获得下一追踪点。在追踪过程中，将经历过的游程左边界或右边界的相应标志值lsign(或rsign)设置为TRUE；
在追踪过程中，判断栅格单元是否有值是根据栅格单元的左下角的点是否在相应的游程上来决定的，本发明采用2×2的栅格窗口，在每个窗口内的四个栅格数据可以唯一地确定下一弧段的搜索方向和该弧段的拓扑关系。这一方法可以加快搜索速度，拓扑关系也易于建立。图8是本发明中基于栅格追踪的几种情况的示意图，如图8所示，对于一个2×2的栅格窗口，可能出现图8所示16种组合方式，分别用0到15序号代表每种情况。
对图8中的16种组合方式进行分组，可大致分4类，见表1：
表1分类描述表