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

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明的一种用于地理信息系统中的空间叠加分析方法和系统进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明的一种用于地理信息系统(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=\frac{(y-P_2.y)\times (P_2.x-P_1.x)}{(P_2.y-P_1.y)}+P_2.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分类描述表

①前2种分别表示当前追踪点为追踪几何体(面)的内部或外部点，而非边界点。这种情形是不需要追踪的情形；

②接下来的8种组合分为1组，均表示追踪边界的转折点(又称边界点)，这些点是追踪过程得到的边界点；

③以下4种情况分为一组，所描述的是追踪过程中的中间点(它们是一条直线段的中间点，虽然他们不直接作为追踪边界点，但是得到第2组情况的必经之点)；

④最后这两种情形是不出现的情况。如果追踪过程中碰到这样的情况可直接跳过。

图9是本发明中基于窗口的栅格追踪的走向图，如图9所示，当前点为B，一点为A，则下一点的走向由当前窗口的走向决定，在图中显然下一点为C。对于其他的走向问题可参照上图的窗口决定。

图10是本发明中一个能覆盖所有组合的基于窗口的追踪实例图，作为一种可实施方式，如图10所示，下面以一个能覆盖所有组合形式的基于2×2的栅格窗口的追踪情况为例，对边界矢量化进行说明，追踪情况见表2。

表2追踪情况表

如图10所示，从行号最小最左边的游程的左边界开始追踪，图中的Ga点为起始点，Gb点为追踪的当前点，则在以Gb为中心的四个格子对应的是第5种情况，下一步要向左走一个格子；然后便是以Fb为当前点，其周围的四个格子对应的是第7种情况，下一步要向上走一个格子；依次类推，一步左，再一步上，直到走到Ah点，此时的情况对应的是第13种情况，下一步要向上走；同样的一直走到Ap点，此时对应的是第8种情况，向右走一个格子后，以Cp为当前点，对应的是第2种情况，要向上走一个格子；同样的右一步上一步交替的走到Gu点，此时情况变成第10种；一直向右走到Qu点，对应的是第9种情况，下一步要向下走；以Qt为当前点对应第3种情况，向右走一个格子；下一步和右一步交替走到Vo点，对应第11种情况，一直向下走到Vg点；以Vg点为当前点对应第6种情况，向左走一个格子，然后对应第4种情况，向下走一个格子；同样依次类推，直到走到Oa点，对应第12种情况，一直向左走到起始点Ga点。当发现走到起始点，认为该多边形矢量化完成了。

上述是对一个独立的多边形的矢量化过程，而对于真实数据，即许多多边形相邻接的，例如地图，用游程插入栅格场的做法就有个弊端：就是相邻多边形都合起来了。每个多边形都应该有自己的属性信息，如果一味的合并，就会使某些结果没有意义。

较佳地，本发明提出了一种解决多边形相邻接问题的方法，提出一种改进的判断栅格单元值的方法，用GridValue函数实现。

由于判断栅格是否有值是根据栅格的左下角的点是否在相应的游程上来决定的，用GridValue这个函数实现的，所以只需要在这个函数上进行修改即可，本发明加入了判断多边形ID信息的限制条件。图11是本发明中邻接多边形的边界追踪示意图，如图11所示(为了更看清楚，把栅格的格子放大画出来)，两个相邻的矩形0号和1号，当矢量化右侧1号矩形时，虽然P这个格子是有值的(因为其在0号多边形内部，所以显然在0号多边形当前行的游程链表上)，但是对1号多边形这个值其实是没有意义的(否则矢量化后两个多边形就合并成一个了)，所以加入了判断多边形ID的限制条件。当走到Fe点的时候，判断Q这个格子是否有值的时候，要看Ed点是否在1号多边形的游程链表上。由于这一点正是两个多边形的相邻边上的点(即公共点)，所以用上述做法，发现这个点在0号多边形的游程上，而误认为没有值。然而这个点也在1号多边形的游程上，所以对点在相邻边上并且判断是在非当前多边形的游程链表上时，再多判断一步这个链表的下一个链表是否在当前多边形上，如果这个条件满足则也认为是有值的。

所述改进的判断栅格单元值的方法，可以通过下述伪代码实现：

bool GridValue(参数)

//判断在实游程链表格式记录的栅格数据场上，L行、C列的值是true/false？

{

  用二分查找法找到L行C列所在的游程链表的号，赋给segno；

    如果没有查到，直接返回false；

  If(segno的游程的所属多边形id！＝当前矢量化的多边形的id)

  {

      //C是否在两个多边形的公共边上

    bool horizonal_sep＝(segno+1＜L行的游程数量&&

                      C＝＝(int)L行第segno个游程的右值&&

                      C＝＝(int)L行第segno+1个游程的左值)；

    if(！hori zonal_sep)

       return false；

    if(L行第segno+1个游程的所属多边形的id＝＝当前矢量化的多边形的

id)  即两个多边形有重合并且此点在重合边上

       return true；

       return false；

  }else

     return true；

 }

使用上述方法，在判断左右边界时，不用再用二分法查找L行C列在哪一段，再判断是否是此段的左右边界，而是直接判断游程的id号是否等于当前多边形的那段游程的左右值，因为判断跟别的多边形比没有意义，这也是加了多边形ID信息的一个优势。

步骤314.重复步骤312-213，直至不存在任何lsign(或rsign)为TRUE的游程，即没有新的追踪初始点时终止循环；

步骤320.面的生成；

面的生成采用现有的技术，即我们将上述追踪到的弧设置一个边界属性(是外或内边界)、归属面的序号(每次产生一个新面，则面序号递增)，在所有弧段追踪完毕后，根据上述属性，建立弧段与面之间的关联关系。

图12是本发明中两个输入图层的示意图，图13是对图12中的两个图层求交、差、并的结果示意图，作为一种可实施方式，如图12和13所示，利用本发明的方法将输入的图层进行叠加分析，假设两个输入图层分别是：中国县界图(大小20.2M，3707个要素集)和3个大圆(大小2M，3个要素集)。这两个图层均是ESRI的shape文件格式。

相应于本发明的用于GIS中的空间叠加分析方法，还提供一种用于GIS中的空间叠加分析系统，图14是本发明的空间叠加分析系统的结构示意图，如图14所示，所述空间叠加分析系统，包括：

矢转栅模块1，用于将输入GIS中的图层的矢量数据转换成栅格数据，并采用游程编码表示；

叠加模块2，用于对该采用游程编码表示的栅格数据执行叠加操作；

叠加图层输出模块3，用于将叠加后的栅格数据再转换成矢量数据，得到经过叠加的图层。

所述矢转栅模块1，包括：

扫面模块11，用于扫描输入GIS的一个图层中的每一个多边形，对每一个多边形再循环读入该多边形的每一个环(外环和内环)，对于每个环，计算y值扫描线与每个弧段的断点的横坐标，以构造断点坐标加入断点栅格场；

游程编码模块12，用于将所述多边形的所有环的断点加入到栅格场，建立经过所述y值扫描线扫描得到的行的游程集合并组织成链表形式；

叠加模块2，包括：

合并插入模块21，用于对两个有重叠的游程单元间合并，否则把当前游程单元插入到适当的位置，并且在插入的过程中把游程单元的属性信息保存下来；

删除模块22，用于在当前行游程链表内删除与待删除游程单元区域重叠的部分；

所述叠加图层输出模块3，包括：

边界追踪模块31，用于追踪经过叠加后的用游程编码表示的栅格数据，得到图层弧段；

拓扑重构模块32，用于根据所述图层弧段，生成图层。

现有技术中，从输出效果上看，多边形的叠加操作类似于多边形裁剪。但是多边形裁剪主要是针对两个(或者数量较少的几个)多边形进行的，在规模化的多边形计算面前无能为力。由于多边形的叠加操作，本质上是多边形所包围的点集的集合操作，因此可以采用“暴力”算法，通过反复调用两个多边形裁剪的算法来完成。但是，这样的算法实现其时间复杂度是O(n²)。经过实际测试，这样的算法在实际应用中实用性很差。因此，如同简单的两个线段求交点的算法并不能推广到平面中n条线段的求交点问题一样，给定两个多边形图层，求叠加结果，必须寻找复杂度更低的算法。

而本发明中，对于输入的图层，设其栅格化扫描间隔为dxy，则栅格场行数一般为某个给定的常数$M=\frac{\mathrm{Height}\left(\mathrm{LayerBox}\right)}{\mathrm{dxy}}$(所以一旦扫描间隔确定后，M就是个常数了)。将一个断点插入已排序的断点链表中需要执行折半查找，其复杂度为logP(P为链表中元素个数)。假设输入的图层(线段数较多的图层)包括了n条线段，则在最坏的情况下，将图层数据进行栅格化的复杂度是M×n×logn，即O(nlogn)。图15是本发明中一图层数据示意图，如图15所示，多边形ABCDE包含5条线段，其中4条在栅格化时需要扫描近似M行。

表3是基于栅格的方法(本发明的方法)和ArcGIS的方法的测试性能对比。根据对表3数据的分析和实验，取栅格步长为0.01的时候，结果要素数量与真实数据的结果的一致性能达到97％以上，此时栅格数量为4794*6234。如果要求数据更精确的话，取栅格步长为0.0075，结果要素数量与真实数据的结果是基本吻合的，此时栅格数量为6369*8289。所以下表中精度1代表步长0.01的测试，精度2代表步长0.0075的测试。

表3基于栅格的方法和ArcGIS的方法的测试性能对比表

图16是本发明的方法与ArcGIS的性能测试对比图，其中，表示本发明的运行时间1，表示本发明的运行时间2，表示ArcGIS运行时间，表示本发明内存占用1，表示本发明内存占用2，表示ArcGIS内存占用，左侧纵坐标表示运行时间(s)，右侧纵坐标表示内存占用(MB)，如图16所示，取上表中某组数据做性能测试对比图。输入图层A和B分别是：中国县界图(大小4.1M，2525个要素)和土地利用图(大小100M，122552个要素)。这两个图层均是ESRI的shape文件格式。下图是基于栅格的方法(本方法运行时间1和内存占用1代表使用精度1的结果，本方法运行时间2和内存占用2代表使用精度2的结果)和ArcGIS的方法的测试性能对比图。

由以图16和表3可以看出：

(1)栅格方法在时间和内存上都有较大的优势，随着栅格步长的精确，结果要素也就越接近矢量方法。就性能而言，在大数据量下，本发明的方法与ArcGIS系统中对应的模块相比，本发明的方法的效率比ArcGIS提高40％-50％左右，而且内存耗费大约是其40％-50％。经分析，主要原因在于本算法采用游程链表的数据结构，节省了叠加分析计算的复杂性和耗时。

(2)此算法的时间和内存会随着栅格步长的减小而增加(因为栅格场行数$M=\frac{\mathrm{Height}\left(\mathrm{LayerBox}\right)}{\mathrm{dxy}}$跟栅格步长成反比)。由于栅格数量对算法执行效率影响较大，所以应在满足精度要求条件下选择尽可能少的栅格数。

图17是本发明方法三个关键步骤运行时间比例图，如图17所示，算法运行时间主要花费在初始化输入图层的栅格场和栅格叠加结果矢量化方面，而真正在做叠加操作的时候，由于只是游程链表的简单操作，所以速度比较快，所用时间不到10％。这个发现为网络地图服务中的空间分析功能实现提供了新的方向与潜力。目前主流的网络地图提供方式(比如google map等电子地图与ArcIMS等网络地图服务器)皆提供了分级预生成的方式，栅格图像通过预生成或第一次请求时动态生成的方式缓存在服务器端。这种情形下，栅格图像不仅是矢量数据的可视化表达，更可视为是矢量数据的栅格化近似。在这种网络地图服务模式下提供叠加分析功能以进行辅助决策是一项非常有意义的工作。本文算法契合了这种应用模式。从测试数据可以看出，在预生成栅格化表达的前提下，本算法的叠加效率很高，这为网络地图服务模式下的空间分析功能性能提升提供了很好的探索方向。

本发明的有益效果在于：

1.本发明提出了一种用于GIS中的空间叠加分析算法，此方法利用栅格算法的优势，能避免计算几何算法的缺点，跟业界领先软件-ArcGIS相比，此算法有较大的优势和好的性价比；

2.本发明提出的基于栅格的叠加算法也同时可以通过根据预定义的地图比例尺事先生成的地图瓦片文件的方法来为网络地图服务器上提供空间分析功能；

3.本发明提出的基于栅格的叠加算法，在预生成栅格化表达的前提下，本算法的叠加效率很高，这为网络地图服务模式下的空间分析功能性能提升提供了很好的探索方向。

通过结合附图对本发明具体实施例的描述，本发明的其它方面及特征对本领域的技术人员而言是显而易见的。

以上对本发明的具体实施例进行了描述和说明，这些实施例应被认为其只是示例性的，并不用于对本发明进行限制，本发明应根据所附的权利要求进行解释。