鲁棒高效的矢量地图叠加分析算法_朱效民

1007 4619 2012 03 448 19Journal of Remote Sensing遥感学报 A robust and eff icient method for vector map overlay ZHU Xiaomin1 2 ZHAO Hongchao1 FANG Jinyun1 1 Institute of Computing Technology Chinese Academy of Sciences Beijing 100190 China 2 Graduate University of Chinese Academy of Sciences Beijing 100049 China Abstract Robust and effi cient internal memory algorithm for vector map overlay is proposed The algorithm employs the plane sweep idea improved from the traditional plane sweep algorithm to calculate the intersection points by which all the special cas es are handled properly Then rings of the results are constructed by the intersection points and the information and original rings with no intersection point are ignored or added to the result as outer rings contour or holes All the generated rings have ID in formation which simplifi es the following two processes fi nding the matching outer ring for each hole and attributes propagation With this algorithm we can determine all the intersection points for any overlay operation immediately not by a one to one loop We implemented the algorithm and the comparisons with the one by one method using spatial access method demonstrated its effi ciency Besides we implemented the whole function including the algorithm and the necessary processes of reading data from and writing data to the disk and compared it with ESRI s ArcGIS by which correctness and effi ciency of our approach are demonstrated Key words vector map overlay plane sweep algorithm polygon overlay intersection union difference CLC number TP301 6 Document code A Citation format Zhu X M Zhao H C and Fang J Y 2012 A robust and efficient method for vector map overlay Journal of Remote Sensing 16 3 448 466 1INTRODUCTION In Geographical Information Systems GIS spatial data is often represented as layers of thematic maps e g land use region alism Each layer describes a certain aspect of the modeled real world Storing map layers separately makes it impossible to di rectly solve topological queries that related to features belonged to different layers Therefore a map overlay strategy is needed which combines two or more maps of different topics into a single new map The goals are to derive new maps to fi nd correlation between the information encoded in different maps and to process complex queries It is very important in GIS it is a basic building block for analyzing geographical data and it is also one of the key operations in a GIS as it allows the user to integrate spatial data As shown in Fig 1 an overlay operation can generate four areas with different County A County B a County BC AR AC Rain Cloud b Weather attributes Area of county A with Rain AR Area of county B with Rain BR Area of county A with clouds AC and Area of county B with clouds BC including two disjoint polygons Map overlay operation is widely used in GIS Generally there are two kinds of map overlay methods raster based or vector based Raster based method is simple but not precise and it may need many spaces to store all rasters In contrast vector based map overlay is precise but diffi cult and most of research on map over BR BC c Overlay d Result Fig 1 A map overlay operation Received 2010 04 13 Accepted 2011 08 26 Foundation National High Technology Research and Development Program of China 863 program No 2009AA12Z226 2011AA120302 First author biography ZHU Xiaomin 1982 male Ph D candidate he majors in spatial analysis spatial computation algorithms and Parallel algorithms E mail zhuxiaomin Corresponding author biography FANG Jinyun 1967 male associate professor his research interests are GIS spatial analysis and Massive multi dimensional data processing E mail fangjy ZHU Xiaomin et al A robust and eff icient method for vector map overlay2 lay focused on the vector based method Although there are many overlay operations the kernel of over lay on two polygon sets In this paper we call them R and B is to divide the area into three parts the area which belongs to both red and blue the intersection area R B the area which only belongs to red the difference area R B and the area which only belongs to blue the difference area B R as shown in Fig 2 The two operations of intersection and difference are the basis Besides the geometric computation we should compute the information of the sources for each area and attach the attributes to the newly pro duced area In summary the input of an overlay operation is two polygon sets in any of which each polygon has a specific ID to identify itself The output is a polygon set and each polygon carries one for difference operation or two IDs for intersection opera tion to indicate the input polygon s which generate it DifferenceIntersectionDifference Fig 2 Essences of overlay 2RELATED WORK As vector map overlay is important and plays an important role in a GIS previous work on vector map overlay techniques has been proposed Some of those techniques are based on topological structures Finke van Roessel 1991 and oth ers Dong et al 2004 Oosterom 1994 Nievergel Kriegel et al 1992 Chan Becker et al 1999 Rodrigues et al 2007 are all based on geometric objects Since the most popular organizational form of spatial data is based on geometric object for example the Shape fi le format and extra cost is needed for the transformation between object based data and Topological based data the techniques based on objects are state of the art For the object based methods there are generally two ideas ei ther computing by a one to one loop adopting some spatial assess methods Dong et al 2004 to avoid the unnecessary computation between disjoint objects or computing all the objects as a whole Oosterom 1994 Nievergelt Kriegel et al 1992 Chan Becker et al 1999 usually adopting a plane sweep methods Shamos Bentley Chan 1994 Andrews et al 1994 The one by one loop methods just compute the overlay between two objects once and they are not effi cient enough because any object may participate in the intersection computing process many times For example the algorithm of Dong et al 2004 used the quad tree to acquire the relational polygons of the current one and conducted polygon clipping one by one using entry exit property In Rodrigues et al 2007 the Terralib also adopted the one by one method using plane sweep method to determine intersections between two poly gons During the implementation of Terralib version 3 2 0 they used the brute force method to get the intersection points consider ing the plane sweep method as a potential improvement Oosterom s algorithm 1994 mainly focuses on the intersection computing method and generates new line segments at an intersec tion point using r tree to store the segments Then it constructs all the segments into rings and fi nds the contour for each hole with the point in ring test making it inefficient Other methods often use the plane sweep idea Nievergelt and reparata 1982 extended Bentley and Ottoman s 1979 plane sweep algorithm and first proposed the region fi nding algorithm by defi ning different actions on different kinds of intersections bend start end or intersec tion Whereas it is just appropriate for two polygons it does not consider special cases such as collinear segments Then Kriegel et al 1992 extended the idea and gave treatment of all kinds of event points which is an extension of Nievergelt and Preparata s different intersections categories However it did not consider the collinear but not coincide segments of the same polygon set or different polygon sets and does not give a proper order for those linear segments Chan Finke Andrews et al 1994 Its kernel idea is to use a virtual vertical line to sweep all the line segments always keeping an above below order of all the active segments and to com pute intersections only between adjacent segments In this paper we just adopt Bentley for ArcGIS function we just list the total time cost We compare the overlay results of the test and all of them are the same In an overall view of time cost our method works better than ArcGIS Table 3 Comparisons on real world data sets Total time the algorithm By this method we get the result as a polygon set each of OperationResult No cost s GeometricMethod which has two IDs to identify its source features ID Then we combine the geometric result with features attributes construct a new layer and then save it This is the whole function of map overlay including geometric computation and geographical computation 5EXPERIMENTAL EVALUATION The most time consuming process is the intersection computing process and then its complexity is O N K logN in which N is the number of polygon edges and K is the number of intersection points The space complexity is O N K The comparison is on real world data sets the Land use dataset and the County dataset China formatted Shape fi le and the feature Intersect D1 D2 Difference D1 D2 Difference D2 D1 Symmetric difference D1 D2 D2 D1 Union D1 D2 1574547563Our method 157454144 ArcGIS 4116361Our method 411132 ArcGIS 27035956Our method 2703134 ArcGIS 31146461Our method 3114133 ArcGIS 1605688675Our method 160568144 ArcGIS 456Journal of Remote Sensing遥感学报2012 16 3 6CONCLUSION In this paper we introduce a new approach for vector map over lay The method uses the vector structure of STL to store the data and uses the improved plane sweep algorithm to compute intersec tions A new method of computing the entry exit property is also proposed for both the general case and the non general cases Also we record the polygon ID information with the intersection the ring and the result polygon and this simplifi es the process of fi nd ing contour for holes and attributes propagation We implement the algorithm with some necessary auxiliary processes and do a com parison with ArcGIS and the comparison indicates its correctness and effi ciency REFERENCES Andrews D S Snoeyink J Boritz J Chan T Denham G Harrison J and Zhu C 1994 Further comparison of algorithms for geometric intersection problems Proceedings of the 6th International Sym posium on Spatial Data Handling 709 724 Becker L Giesen A Hinrichs K H and Vahrenhold J 1999 Algorithms for performing polygonal map overlay and spatial join on massive data sets Proceedings of the International Symposium on Advances in Spatial Databases Hong Kong China 270 285 Bentley J L and Ottmann T A 1979 Algorithms for reporting and counting geometric intersections IEEE Transactions on Comput ers 28 9 643 647 Chan E P F and Ng J N H 1997 A general and efficient implementa tion of geometric operators and predicates Proceedings of the 5th International Symposium on Advances in Spatial Databases Berlin Germany 69 93 Chan T M 1994 A simple trapezoid sweep algorithm for reporting red blue segment intersections Proceedings of 6th Canadian Con ference on Computational Geometry Saskatoon Saskatchewan Canada 263 268 Dong P Li J P Bai Y Q Qian Z G and Yang C J 2004 Vector map overlay algorithm based on improved quadtree indexing Journal of Computer Aided Design and Computer Graphics 16 4 530 534 Feito F Torres J C and Ure a A 1995 Orientation simplicity and in clusion test for planar polygons Computers and Graphics 19 4 595 600 Finke U and Hinrichs K 1993 A spatial data model and a topological sweep algorithm for map overlay Lecture Notes in Computer Science 692 162 177 Kriegel H Brinkhoff T and Schneider R 1992 An efficient map over lay algorithm based on spatial access methods and computational geometry Geographic Database Management Systems Berlin Springer Verlag 194 211 Liu Y K Gao Y and Huang Y Q 2003 An efficient algorithm for poly gon clipping Journal of Software 14 4 845 856 Nievergelt J and Preparata F P 1982 Plane sweep algorithms for inter secting geometric figures Communications of the ACM 25 10 739 747 Oosterom P 1994 An r tree based map overlay algorithm Proceed ings of EGIS 94 Paris EGIS Foundation 318 327 Rodrigues V L Andrade M V A de Queiroz G R and de Magalh es M A 2007 A more efficient method for map overlay in terralib David Junior C A Monteiro A M V eds Advances in Geoinformatics Berlin Springer 37 46 Shamos M I and Hoey D 1976 Geometric intersection problems Proceeding of 17th IEEE Annual Symposium on Foundation of Computer Science Houston Texas 208 215 Van Roessel J W 1991 A new approach to plane sweep overlay topo logical structuring and line segment classification Cartography and Geographic Information Science 18 1 49 67 Wang Z Q Xiao L J and Hong J Z 1998 Simplicity orientation and inclusion test algorithms for polygons Chinese Journal of Com puters 21 2 183 187 朱效民 等 鲁棒高效的矢量地图叠加分析算法457 鲁棒高效的矢量地图叠加分析算法 朱效民1 2 赵红超1 方金云1 1 中国科学院计算技术研究所 北京 100190 2 中国科学院 研究生院 北京 100049 摘 要 提出了一个鲁棒高效的内存矢量地图叠加分析算法 采用改进的平面扫描算法计算交点 解决了重叠边 交 点位于端点等所有特殊情形 利用交点及其携带的信息来构造结果环 并且将没有产生交点的输入环忽略 或者增加 到结果的外环 或内环 集合中去 所有结果环都带有标识码 增加该标识码信息可以简化后续的两个过程 内外环的匹 配以及属性的继承 与一一循环方法相比 本文方法对任何叠加操作可以一次计算得到所有的交点 此外还实现了叠 加分析操作 并且用一组真实地理数据的不同操作与ESRI的ArcGIS的叠加分析操作进行了比较 计算结果的要素数完 全一致 计算时间耗费约为ArcGIS时间耗费的50 60 关键词 矢量地图叠加 平面扫描算法 多边形叠加 多边形交并差 中图分类号 TP301 6文献标志码 A 引用格式 朱效民 赵红超 方金云 2012 鲁棒高效的矢量地图叠加分析算法 遥感学报 16 3 448 466 Zhu X M Zhao H C and Fang J Y 2012 A robust and effi cient method for vector map overlay Journal of Remote Sensing 16 3 448 466 1引言 在地理信息系统 GIS 中 空间数据经常用带有主题 的地图图层来表示 如土地利用 行政区划等 每个图 层描述了一个建模的真实世界的特定主题 如果将地图 图层分开存储 则无法直接处理涉及不同图层要素的 拓扑查询 这就需要一个地图叠加操作 地图叠加操 作是将多个不同主题的地图合并生成一个新的地图 其目标是通过产生新的地图 找到地图内要素的相互 关系 并且进行复杂的查询 由于地图叠加允许用户 集成空间数据 因此它在GIS中非常重要 是分析地理 数据的基础组成模块 也是GIS中最重要的操作之一 如图1所示 一个叠加操作能够产生4块具有不同属性 的区域AR BR AC和BC 包含两个相离的多边形 其 中A B代表县域 R C分别代表降雨 阴天 地图叠加在GIS中应用广泛 有栅格和矢量两种 类型的叠加分析方法 栅格方法简单但是不够准确 并且它可能需要很多的内存空间来存储所有的栅格数 据 基于矢量的地图叠加方法精确但难于实现 因此 大多地图叠加的研究都是基于矢量方法的 尽管有许多叠加操作 交 并 差 对称差等 但是两个多边形集合 本文中称为R和B 的叠加的核 心即将两组多边形覆盖的区域区分为3部分 仅属于 R的部分 R B之差R B 既属于R又属于B的部分 R B的公共交部分 R B 仅属于B的部分 B R 之差B R 如图2所示 因此交 差操作是基础 其 他操作都是交 差操作的组合 叠 加 分 析 除 几 何 计 算 外 即 叠 加 的 图 形 显 示 结 果 还计算结果多边形的来源 并且将其源信息附 加到新生成的区域 叠加操作的输入是两个多边形集 合 每个集合中的多边形有其唯一的标识码来标记自 身 其输出为一个多边形的集合 结果集合中的每个 多边形有一个 对于求差操作 或者两个标识码 对于求 交操作 标记该结果多边形的源多边形 收稿日期 2010 04 13 修订日期 2011 08 26 基金项目 国家高技术研究发展计划 863 计划 编号 2009AA12Z226 2011AA120302 第一作者简介 朱效民 1982 男 博士研究生 主要研究方向为空间计算 算法及并行算法设计 E mail zhuxiaomin 通信作者简介 方金云 1967 男 副研究员 博士生导师 主要研究方向为地理信息系统 空间计算和海量空间数据处理 E mail Fangjy 458Journal of Remote Sensing遥感学报2012 16 3 A B a 县 BC AR AC 降雨 阴天 b 天气 几何体对 然后依次计算 另一种采用平面扫描方法 Shamos和Hoey 1976 Bentley和Ottmann 1979 Chan 1994 Andrews 等 1994 一次整体性进 行计算 Oosterom 1994 Nievergelt和Preparata 1982 Kriegel 等 1992 Chan和Jimmy 1997 Becker 等 1999 一一循环的方法每次只在两个几 何体之间计算 每个几何体可能参与多次的交点计 算过程 因此不够高效 比如 董鹏等人的算法 董 鹏 等 2004 使用四叉树得到可能与当前多边形有交 点的多边形 然后基于出入属性 一一进行多边形 的裁剪操作 在 Rodrigues 等 2007 中 Terralib 也 BR BC c 叠加 d 结果 图1 一个地图叠加操作 差交差 图2 叠加的本质 2相关研究 众多矢量地图叠加分析方法中 Finke和Hinrichs 1993 Van Roessel 1991 的方法是基于拓扑结构 的 董鹏 等 2004 Oosterom 1994 Nievergelt 和 P r e p a r a t a 1 9 8 2 K r i e g e l 等 1 9 9 2 C h a n 和 Jimmy 1997 Becker 等 1999 Rodrigues 等 2007 的方法是基于几何对象的 而最流行的空间数据组 织和存储方式是基于几何对象的 比如Shape File格 式 因此最流行的叠加分析方法是基于几何对象 的 该类方法无需在两种组织和存储形式下对数据 进行转化 基于几何体的方法通常有两种 一种采用如四叉 树 董鹏 等 2004 空间获取方法 得到需要计算的 采用了一一计算的方法 不过利用平面扫描方法来 计算两个多边形之间的交点 然而在其具体实现中 Terralib version 3 2 0 仍然采用了暴力求解法计算 交点 将平面扫描方法视为一个潜在的优化方法 Oosterom 1994 的算法主要关注交点的计算方法 在 每个交点处产生新的线串 采用R树存储线串 然后 将所有的线串组合生成环 并且采用点在环内测试 方法进行内外环的匹配 其他的方法通常采用平面 扫描思路 Nievergelt和Preparata 1982 拓展了Bentley 和Ottmann 1979 的平面扫描算法并且首次提出了通过 在不同的交点处定义不同的操作的区域查找算法 之后Kriegel 等 1992 扩展了这个思路并且给出了 各种事件点的处理方法 以上两个平面扫描思路的 算法都有其局限性 Nievergelt和Preparata的算法仅 适合两个多边形的叠加 Kriegel等人的方法对重叠 但不重合的线段等特殊情形无法正确处理 Chan和 Ng 1997 Becker 等 1999 通过增加一些针对特殊情 况的处理方法完善了Nievergelt和Preparata的方法 这些仍然没有覆盖所有的特殊情形 此外 一一循 环的方法的交点计算与扫描线算法相比效率较低 但内外环匹配 属性继承等过程无需计算 而这两 个过程对于整体进行计算的方法则是必需的且需要 额外的计算 一般采用点在环内测试 点在多边形 内的包含性测试方法 如 Oosterom 1994 Finke和 Hinrichs 1993 本文提出了一个新的矢量地图叠加分析算法 通 过一次扫描所有线段计算得到交点 并且改进Bent ley和Ottmann算法来处理所有的特殊情形 增强了算 法的鲁棒性 此外 还记录标识码信息以完成内外环 匹配及属性继承 避免了复杂的点在环内或点在多边 形内的包含性测试 朱效民 等 鲁棒高效的矢量地图叠加分析算法459 3基础的定义和计算方法 3 1 平面扫描算法 平面扫描算法是一个高效的计算交点的方法 该方法首先由Shamos和Hoey 1976 提出 Bentley和 Ottmann 1979 改进了这个算法 之后 该算法得到 了进一步的改进 Chan 1994 Andrews 等 1994 其核心思路是用一条虚拟的竖直的直线扫描所有的线 段 始终保持所有活动的线段按照由上而下的顺序 交点计算仅发生在相邻的线段之间 3 2 环 多边形 的方向 环的方向有顺时针和逆时针两种 对于一个多边 形而言 外环经常与内环的方向相反 如此多边形的 内部区域总是位于任何多边形边的同一侧 不管该边 是多边形的外环还是内环的边 如图3所示 多边形 的外环是顺时针 内环是逆时针 因此内部区域位于 每条边的右侧 传统的判断环的方向的方法为带符号的面积法 Feito 等 1995 本文采用的是王志强等人 1998 提 出的新方法 即计算具有最大 最小 x y值的坐标及 其前一个 后一个坐标构成的三角形的方向 如图3 所示 通过计算三角形 3个端点为 D0 D D1 的 方向得到整个环的方向 边形的交组成的边仍然保持输入多边形边的方向 差 的结果 R B 中属于R的边与输入的方向相同 而属 于B的边则与输入时的方向相反 如图4中的B1 任何两个合法的多边形集合 合法多边形集合是 指 集合中的每个多边形都没有自相交区域 集合 内任意两个多边形没有重叠区域 的交有两个源多边 形 一个属于R 另外一个属于B 任何一个R与B的 差结果多边形只属于一个R 但是可能被许多B多边 形相减得到 因此对于差结果而言 本文只考虑其源 标识码而忽略作为减数的若干裁剪多边形 如果采用 扫描而非依次计算的方法计算R B之差 需要记录 B多边形集合自身的 交点 不是交点而仅仅是接触 点 如图5 分别图示了R B求差操作的3种情形 只利用R B之间的交点就可以求得差 需要记录多 R2 R1 B1 B2 RB 图4 利用内外边构造生成结果 B A C D1D0 D 图 3 多边形的方向及其内部区域示意图 a 利用R B之间的交点来求差 3 3 基本规则 两个多边形集合R 如图4所示的实线多边形 和 B 虚线多边形 交的区域是由属于R且位于B内部的 边 R1 以及属于B且位于R内部的边 B1 两部分组成 R B之差结果是由属于R且位于B外部的边 R2 与属 b 记录B自身的接触点来 构造R B之差 c 记录B自身的接触点来构 造R B之差 B中包含不与 R相交的多边形 于B且位于R内部的边 B1 组成的 如图4所示 如果两个输入多边形的方向相同 那么这两个多 RB 图5 不同的差操作 460Journal of Remote Sensing遥感学报2012 16 3 边形集合B自身的接触点来构造R B之差 且生成差 结果的B多边形都与R相交 需要记录多边形集合B 自身的接触点来构造R B之差 且B中不包含与R相 交的多边形 需要利用与R的真正交点以及B自身的 伪交点 将B多边形的边拆分成线串 4 内存算法 4 1 出入属性计算 任何一条边都是多边形的边界 是多边形内外区 域的分割 任何一条与多边形边相交的线段必然从交 点处 进入 或者 走出 多边形 因此任何一个交 点都有出入属性 该属性与线段的走向相关 如图6 所示 线段的走向如箭头所示 交点P是一个入点 如果线段如虚线方向所示 那么P就是一个出点 本 文的出入属性的含义为 R多边形在交点处出入B多 边形 多边形的内部 实际上 在交点的位置 沿着线段 前进的方向前进足够小的距离的点P必然位于多边 形的内部 以上描述主要针对通常情形 对于特殊情形 交 点恰恰是端点 如图8所示 其出入属性的计算遵循 下面方法 从计算角度而言 多边形的两条边及其延 长线将平面区域分为4部分 LL LR RR和RL L 即为左侧 R即为右侧 如果两条边的方向与多边形 E V B a 交点出入属性计算 V P P b 交点是一个入点 但末端点 并没有位于B内部 RB 图6 出入属性 RB 图7 新的交点出入属性计算方法 LR 如刘勇奎等人 2003 所述 传统的交点出入属性 计算需要一次点的包含性测试 本文给出了一个新的 RR LL 仅根据产生交点的两线段计算出入属性的方法 这需 要一个前提即多边形都是顺时针 因此其内部区域都RL 位于任何一条边的右侧 如图7 a 所示 实线段 虚 线段分别属于R B两个多边形集合 通过计算末端 点V位于B多边形的线段的左侧还是右侧即可得到交 点的属性 如果V位于B的线段的右侧 那么交点是 一个入点 否则即为出点 上述 左右 的计算可以通过计算B线段的始点 B 终点 E 以及R线段的终点V构成的三角形的方向 来判断 不过 位于B线段的右侧不意味着位于B线 段所在的多边形的内部 如图7 b 所示 交点是一 个入点 但是线段的末端点 V 并不位于B线段所在 a 多边形的边将平面分为4个区域 R2 R1 b 顺时针 c 逆时针 RB 图8 特殊情形的交点出入属性计算 朱效民 等 鲁棒高效的矢量地图叠加分析算法461 的方向一致 图8 b 是顺时针的 那么RR意味着为内 部