CN108985306B

CN108985306B - Method for extracting intersected polygons based on improved boundary algebra

Info

Publication number: CN108985306B
Application number: CN201810731268.0A
Authority: CN
Inventors: 周琛; 李满春; 陈振杰; 姜朋辉; 程亮
Original assignee: Nanjing University
Current assignee: Nanjing University
Priority date: 2018-07-05
Filing date: 2018-07-05
Publication date: 2020-03-31
Anticipated expiration: 2038-07-05
Also published as: CN108985306A

Abstract

The invention relates to an extraction method of an intersected polygon based on an improved boundary algebra, which comprises the following steps: numbering the polygons in all layers in sequence; calculating MBR, array containing all layershDstDS、pIDArrayAndRLEGrouprespectively storing the attribute value, the polygon ID and the run length of the grid unit; rasterizing all polygons in sequence by using a boundary algebra algorithm, wherein the attribute value of each polygon is 1 in the rasterizing process; in an arrayhDstDSAcquiring a grid unit contained in the current polygonal MBR, reading and acquiring attribute values of the grid unit line by line, and performing corresponding processing according to different attribute values; slave arrayRLEGroupExtracting corresponding sets of intersecting polygons from stored runs, i.e. arrays in each runpGroupI.e. to an intersecting set of polygons. The method is low in calculation complexity, and is particularly suitable for extracting the intersected polygons of the large-scale polygon data set.

Description

Method for extracting intersected polygons based on improved boundary algebra

Technical Field

The invention relates to an extraction method of an intersected polygon based on an improved boundary algebra, belonging to the technical field of image processing.

Background

The polygon intersection calculation is an important calculation type in a Geographic Information System (GIS), is widely applied to GIS space analysis types such as polygon topology inspection, buffer generation, overlay analysis and the like, and has typical calculation complex and dense algorithm characteristics. The general steps of this type of spatial computation include: (1) traversing all polygon data sets, searching and determining intersected polygons; (2) extracting the intersected polygons and completing the solution of the results of the intersected parts. At present, a great number of scholars develop research aiming at the processing efficiency and precision of the algorithm, and the algorithm can be divided into a vector-based method and a grid-based method according to different using methods.

The vector-based method is mainly described in a space-oriented manner, and the judgment of the space relationship among the polygonal objects needs to be performed by positioning, analyzing and retrieving based on a computational geometry algorithm. The basic processing process comprises the following steps: the method comprises the steps of representing the space shape of a polygon by using a minimum Bounding Rectangle (MBR for short) of the polygon, filtering out candidate polygons which are possibly intersected, refining the candidate polygons to determine an intersected polygon group, and cutting the intersected polygon group. Although the computing algorithm principle of different polygon spaces is different, the method has the same algorithm characteristic that the essential operations are to filter out polygons with an intersection relation and cut a polygon point set. Thus, a core problem of computationally intensive polygon space computation is clipping between polygons and determination of intersecting polygons. The existing mature polygon clipping algorithm is mainly performed on two polygons; aiming at massive polygon data, the general method is to adopt a brute force algorithm to traverse the source data one by one to judge the spatial relationship between the current polygon and other polygons, and repeatedly call a cutting algorithm to cut intersecting polygons, so that not only is the computational complexity high and the efficiency low, but also the method is difficult to be applied to the intersection calculation of a large-scale polygon data set. However, the number of judgments of subsequent polygons can only be reduced by using the spatial index, and the formed polygon group still contains a large number of disjoint polygons, so that the method still wastes a large amount of time on the computation of the disjoint polygons, thereby reducing the computation efficiency.

The raster data model is a space object representation method that is equally important as the vector data model. Compared with the vector geographic data format, the grid data is high in distribution rule and independence. Therefore, the space calculation efficiency can be greatly improved by adopting the grid-based calculation method. The grid-based computing method is to implement spatial analysis in a simple grid image through simple algebraic operation of attribute values between grid cells. The general calculation process of the method comprises three steps: (1) rasterizing and filling the vector polygon layers respectively to obtain rasterization results of the two processing layers; (2) performing raster calculation on the rasterized images of the two image layers according to the attribute values of the rasterized images to obtain a calculation result of an appointed space; (3) and carrying out grid vectorization on the calculation result, extracting the polygon boundary and the island hole to which the polygon boundary belongs, and constructing a polygon topological relation so as to finish the calculation. At present, more scholars adopt the method to realize the spatial calculation of the polygon data, but the existing method can only process the grid calculation such as the intersection, the combination and the like of the polygons in two or more layers, and can not process the intersection calculation under the condition that the polygons in the same layer are mutually overlapped.

Disclosure of Invention

The invention aims to solve the technical problems that: the defects of the technology are overcome, and the method for extracting the intersected polygons based on the improved boundary algebra is low in calculation complexity and suitable for large-scale polygon data sets.

In order to solve the technical problems, the technical scheme provided by the invention is as follows: an extraction method of an intersected polygon based on an improved boundary algebra comprises the following steps:

step one, numbering polygons in all layers in sequence, wherein each polygon corresponds to a unique polygon ID;

step two, calculating MBR (Membrane bioreactor) containing all layers, generating an array hDstDS with an initial value of 0 to store the attribute value of a grid unit according to the grid size, wherein the size of the array hDstDS is nYSize row multiplied by nXSize column, and simultaneously generating an array pIDAlay with the same size as the array hDstDS and an initial value of-1 to store the polygon ID;

generating a new array RLEGroup to store the run length; wherein, the run length with attribute value pNum of grid cell in a grid line is defined as RLE { StartX, EndX, LocateY, (pID)₀，pID₁，…，pID_pNum-1) pNum }; wherein StartX and EndX are the abscissa of the starting grid cell and the ending grid cell in the grid row with property value pNum, LocateY is the ordinate of the grid row, (pID)₀，pID₁，…，pID_pNum-1) For storing arrays of intersecting polygon IDs pGroup

Thirdly, rasterizing all the polygons in the first step in sequence by using a boundary algebra algorithm, endowing each polygon with an attribute value of 1 in the rasterizing process, and updating the value of the array hDstDS as the attribute value of a grid unit;

in the process of rasterizing all polygons in sequence, the process when rasterizing one polygon is as follows: acquiring a grid unit contained in the current polygon MBR in an array hDstDS, reading the grid unit line by line, and calculating the storage position location of the grid unit in hDstDS as locateY × nXSize + StartX and acquiring the attribute value of hDstDS [ location ] if the grid unit is located inside the polygon for any grid unit (StartX, locateY) in a grid line with the ordinate as locateY;

according to the attribute value of the grid cell, the processing is respectively carried out in three cases:

1) if the value of the property of the grid cell is <2, the grid cell is only located inside one polygon and is not processed;

2) if the attribute value of the grid cell is 2, the grid cell is located inside two intersecting polygons; continuing to scan the subsequent grid cell in the grid line with the current grid cell as a starting point to form a grid sequence from the grid cell (StartX, LocateY) to the grid cell (EndX, LocateY), wherein StartX is less than or equal to EndX; any grid unit (LocateX) in the grid sequence_i，LocateY)，StartX≤LocateX_iEndX, must be located inside the polygon, and satisfies the following conditions:

thus, all grid cells in the grid sequence are located inside two intersecting polygons; meanwhile, the ID of the two intersected polygons are pIDARray [ locate ] and pCurrentID respectively; creating a new run to store the trellis sequence, which may be denoted as { StartX, EndX, LocateY, (pIDARrray [ locate ], pCurrentID),2}, and storing the run into the array RLEGroup;

3) if the attribute value of the grid cell is value>2, the grid cell is located at three or threeInside the above intersecting polygons; continuing to scan the subsequent grid cell in the grid line with the current grid cell as a starting point to form a grid sequence from the grid cell (StartX, LocateY) to the grid cell (EndX, LocateY), wherein StartX is less than or equal to EndX; any grid unit (LocateX) in the grid sequence_iLocateY), wherein StartX is not more than LocateX_iEqual to or less than EndX, must be positioned in the polygon, and satisfies the following conditions:

thus, all grid cells in the grid sequence are positioned inside value intersecting polygons, wherein one polygon ID is pCurrentID, and value-1 polygon IDs are stored in the existing run;

slave grid cell (max (StartX)_jStartX), LocateY) to grid cell (min (EndX)_jEndX), locateY) to search for a run RLE containing value-1 polygon IDs_j{StartX_j,EndX_j,LocateY,(pID₀,pID₁,…,pID_value-2) Value-1, which is first broken down into two new runs, respectively { StartX }_j,max(StartX_j,StartX)-1,LocateY,(pID₀,pID₁,…,pID_value-2) Value-1 and { min (EndX)_j,EndX)+1,EndX_j,LocateY,(pID₀,pID₁,…,pID_value-2) Value-1, where min (EndX)_j,EndX)≤EndX_j，StartX_j≤min(StartX_jStartX) and places two new runs into the array rlegorup, and then RLE the original run_jDelete from array RLEGroup; while generating a new run max (StartX)_j,StartX),min(EndX_j,EndX),LocateY,(RLE_jpGroup, pCurrentID), value }, and put it into an array RLEGroup;

in the process of rasterizing all polygons in sequence, after rasterization of one polygon is completed, the array pIDARay is updated, and the current polygon ID value is assigned to the grid unit positioned in the current polygon in the array pIDARay;

and step four, extracting corresponding intersected polygon groups from the runs stored in the array RLEGroup, namely that the array pGroup in each run corresponds to an intersected polygon group.

Compared with vector polygon data, the grid data has a regular data structure, the grid data is more suitable for GIS space analysis, and more polygon data have higher calculation efficiency. Therefore, converting vector polygons into grid data to quickly search for intersecting polygons may improve GIS spatial analysis efficiency.

The boundary algebra is a polygon rasterization method based on an integral idea, and polygon attribute values are assigned to grid units in the polygon and on the boundary through simple addition and subtraction algebra operation, so that polygon rasterization is realized. The vector polygon is generally composed of an outer ring and a plurality of inner rings, and the processing procedure of the boundary algebra on any polygon includes the following steps, as shown in fig. 1: (1) integrally rasterizing a surface area covering the polygon, and initializing a grid unit value of a grid array to be zero; (2) starting from any point on the boundary, searching the boundary of the polygon outer ring in the anticlockwise direction, and searching the boundary of the polygon inner ring in the clockwise direction; (3) when the direction of the boundary line segment is an uplink, adding an attribute value a to all grids with the same row coordinate on the left side of the line segment; when the direction of the boundary line segment is downward, all grids with the same line coordinate on the left side of the line segment are subtracted by the attribute value a; when the boundary line segment is parallel to the grid line, no operation is performed; (4) and (4) repeating the step (3) until all the boundaries of the polygon are processed. After the above calculation is completed, the grid cell with the attribute value a is the calculation result of the current polygon rasterization. Compared with other polygon rasterization algorithms, the boundary algebra has the advantage of higher calculation efficiency; more importantly, the algorithm can quickly identify the grid cells located at the intersection of the polygons.

The invention utilizes the thought of the boundary algebra, the attribute value given to the grid unit is fixed as 1, and in the processing process of the boundary algebra, if the attribute value of the grid unit is 1, the grid unit is only positioned in one polygon; if the grid cell attribute value is n (n >1), then the grid cell is located within the intersection of n polygons, as shown in FIG. 2. The invention utilizes the characteristic to realize the quick positioning of the grid unit at the intersection part of the polygon.

Drawings

The invention will be further explained with reference to the drawings.

FIG. 1 is a schematic diagram of a prior art boundary algebra.

Fig. 2 is a schematic diagram of the principle of the present invention.

FIG. 3 is a schematic flow chart of an embodiment of the present invention.

Detailed Description

Examples

As shown in fig. 1, the method for extracting an intersected polygon based on an improved boundary algebra of this embodiment includes the following steps:

step one, numbering polygons in all layers in sequence, wherein each polygon corresponds to a unique polygon ID; fig. 1 shows a total of three polygons, the IDs of which are 0, 1 and 2, respectively.

Step two, calculating MBR (minimum bounding rectangle) containing all layers, generating an array hDstDS with an initial value of 0 to store the attribute value of a grid unit according to the grid size, wherein the size of the array hDstDS is nYSize row multiplied by nXSize column, and simultaneously generating an array pIDARray with the same size as the array hDstDS and an initial value of-1 to store the polygon ID. The size of the array hDstDS and the array pIDArray in this example is 6 × 6.

Generating a new array RLEGroup to store the run length; wherein, the run length with attribute value pNum of grid cell in a grid line is defined as RLE { StartX, EndX, LocateY, (pID)₀，pID₁，…，pID_pNum-1) pNum }; wherein StartX and EndX are the abscissa of the starting grid cell and the ending grid cell in the grid row with property value pNum, LocateY is the ordinate of the grid row, (pID)₀，pID₁，…，pID_pNum-1) To store an array pGroup of intersecting polygon IDs.

in the process of rasterizing all polygons in sequence, the process when rasterizing one polygon is as follows:

during rasterization of a polygon, the assigned grid cell in its MBR is located at the polygonThe process of determining the interior of the polygon, i.e., the grid cell is located within the polygon, is performed simultaneously with the rasterization process, and no additional calculation process is required. 3) If the attribute value of the grid cell is value>2, the grid unit is positioned inside three or more crossed polygons; continuing to scan the subsequent grid cell in the grid line with the current grid cell as a starting point to form a grid sequence from the grid cell (StartX, LocateY) to the grid cell (EndX, LocateY), wherein StartX is less than or equal to EndX; any grid unit (LocateX) in the grid sequence_iLocateY), wherein StartX is not more than LocateX_iEqual to or less than EndX, must be positioned in the polygon, and satisfies the following conditions:

FIG. 3(a) shows rasterized polygon 0, which has no intersecting polygons. FIG. 3(b) shows rasterized polygon 1, where there is a raster sequence of two raster lines in the MBR for polygon 1 where polygon 0 and polygon 1 intersect, by a run RLE₀And RLE₁And recording is carried out. FIG. 3(c) shows rasterized polygon 2, where a raster sequence of intersecting polygons exists in the MBR for polygon 2, as by a run RLE₂、sub-RLE₁、RLE₃、RLE₄And RLE₅And recording is carried out.

Considering that different runs may store repeated polygon groups, the polygon group of the type needs to be deleted after the extraction of the intersected polygon group is completed, so as to avoid subsequent repeated calculation. When two polygon groups store polygons with the same number and ID, one polygon group should be deleted; when a polygon set stores polygons that are a subset of another polygon set, the polygon set should also be deleted. In the formed polygon groups, if one polygon group includes two polygons, the two polygons intersect; if a polygon group includes more than two polygons, every two polygons intersect.

In the process of rasterizing all polygons in sequence, after rasterization of one polygon is completed, the array pIDARay is updated, and the current polygon ID value is assigned to the grid unit positioned inside the current polygon in the array pIDARay.

Step four, extracting corresponding intersecting polygon groups from the runs stored in the array RLEGroup, that is, the array pGroup in each run corresponds to an intersecting polygon group, as shown in FIG. 3 (d). After the mutually independent polygon groups are extracted, the calculation of the intersection result can be realized for each polygon group to obtain the intersection result polygon.

In order to verify the extraction method of the embodiment, the inventor uses 2 pieces of real land use data and 1 piece of simulation data to extract the intersected multi-shape edge group, then performs intersection calculation, and compares the intersection calculation with the intersection calculation performed by the traditional spatial index method. Wherein, the data 1 is 2009 land utilization status data of Shanghai city, the data volume is 1.04GB, and the number of polygons is 1,371,765; data 2 is Shanghai city county-level administrative division data, and its data size is 376KB and the number of polygons is 17. Neither data 1 nor data 2 has intersecting polygons within the graph layer. Data 3 is analog data, primarily by moving, rotating, etc. a number of polygons in data 1, thereby forming intersecting polygons in data 3. The data size of data 3 is 1.04GB, the number of polygons is 1,371,765, and 163,934 polygon groups are included.

And the verification adopts a high-performance IBM computing node as a computing tool, and the hardware configuration of the IBM computing node is as follows: CPU 2, whose specification is Intel (R) Xeon (R) CPU E5-2620 (dominant frequency 2.00GHz, six-core twelve-thread); the memory is 16GB (4 memory banks with 4GB, and the specification is DDR3RDIMM1600 MHz); the hard disk is 2TB (2 1TB hard disks, the specification is 6Gbps 2.5' 7.2Krpm NLSAS), and the network is an integrated double-port gigabit Ethernet. Software configuration: the operating system is Centos Linux 6.3, and the file system is a Lustre system. The reading and writing of the vector and the raster data are realized by adopting an open source geographic space data conversion library GDAL (geospatial data extraction library); the polygon group space geometry calculation is realized by open source algorithm library GEOS (GeometryEngineeropen Source).

The result shows that the running time of the method of the embodiment is reduced by about 50% compared with the traditional spatial index method, namely, the method of the embodiment has higher calculation efficiency compared with the traditional method.

The present invention is not limited to the specific technical solutions described in the above embodiments, and other embodiments may be made in the present invention in addition to the above embodiments. It will be understood by those skilled in the art that various changes, substitutions of equivalents, and alterations can be made without departing from the spirit and scope of the invention.

Claims

1. An extraction method of an intersected polygon based on an improved boundary algebra comprises the following steps:

generating a new array RLEGroup to store the run length; wherein, the run length with attribute value pNum of grid cell in a grid line is defined as RLE { StartX, EndX, LocateY, (pID)₀，pID₁，…，pID_pNum-1) pNum }; wherein StartX and EndX are the abscissa of the starting grid cell and the ending grid cell in the grid row with property value pNum, LocateY is the ordinate of the grid row, (pID)₀，pID₁，…，pID_pNum-1) For storing an array of intersecting polygon IDs pGroup;

2) if the attribute value of the grid cell is 2, the grid cell is located inside two intersecting polygons; continuing to scan the subsequent grid cell in the grid line with the current grid cell as the starting point to form a grid sequence from the grid cell (StartX, LocateY) to the grid cell (EndX, LocateY), wherein StartX is less than or equal to EndX, respectively; any grid unit (LocateX) in the grid sequence_i，LocateY)，StartX≤LocateX_iEndX, must be located inside the polygon, and satisfies the following conditions:

3) if the attribute value of the grid cell is value>2, the grid unit is positioned inside three or more crossed polygons; continuing to scan the subsequent grid cell in the grid line with the current grid cell as a starting point to form a grid sequence from the grid cell (StartX, LocateY) to the grid cell (EndX, LocateY), wherein StartX is less than or equal to EndX; any grid unit (LocateX) in the grid sequence_iLocateY), wherein StartX is not more than LocateX_iEqual to or less than EndX, must be positioned in the polygon, and satisfies the following conditions:

2. The method of claim 1, wherein the method comprises: traversing all the polygon groups in the run, and deleting one of the two polygon groups when the polygons with the same storage number and the same polygon ID exist in the two polygon groups; when a polygon stored in one polygon group is a subset of another polygon group, the polygon group as the subset is deleted.