CN104915991A - Spatial subdivision method for intersecting linear geographic elements - Google Patents

Abstract

The invention discloses a spatial subdivision method for intersecting linear geographic elements. The method includes the following steps that: intersecting linear elements are dispersed so as to form point sets based on the principle of a Voronoi distribution model; point addition processing is performed on intersecting positions of every two intersecting linear elements; Delaunay triangulation and Voronoi diagram generation are performed based on the point sets; and a generated Voronoi diagram is clipped, and sub regions generated by point sets which belong to the same linear element are merged, so that sub regions of different linear elements can be obtained, and therefore, spatial subdivision of the intersecting linear geographic elements can be realized.

Description

A kind of space partitioning method for crossing wire geographic element

Technical field

The invention belongs to geographical information technology and application, be specifically related to a kind of method that crossing wire geographic element based on intersection auxiliary point carries out geographical space subdivision.

Background technology

In geographical information technology and application, the more solution relating to allocation optimum problem.And in allocation optimum, first it is envisaged that the determination of each key element coverage separately.How to accomplish the coverage rationally determining each key element? at present, using maximum is space partitioning method, and correlative study mainly concentrates on the uniformly subdivision of a key element, and less to the space partitioning method research based on Linear element.

Wang new life etc. proposes a kind of Approximation Algorithm for Constructing Voronoi Diagrams with General Generators (" Approximation Algorithm for Constructing Voronoi Diagrams with General Generators based on GIS ", Wang new life, Liu Jiyuan, the village are generous, Jiang Youhua, geographical science is in progress, 2004,23 (4)), for the subdivision of line feature, mainly by line feature is separated into point set at equal intervals, then generate the first Voronoi figure of generation arbitrarily and realize.But the situation being difficult in the method avoid intersection point to be considered to belong to multiple generation unit and Voronoi figure cannot be generated, and need realize by ArcGIS platform, automaticity is lower.Zhang Hui etc. devise a kind of algorithm (" a kind of practical algorithm building the first Voronoi figure of generation arbitrarily " building the first Voronoi figure of generation arbitrarily, Zhang Hui, Hu Wei, Pu Yingxia, king tie minister, geography and Geographical Information Sciences, 2011,27 (4)).The method and Wang Xin generation method similar, and for cross spider key element, devise and a kind ofly cross spider to be interrupted and antinode does the disposal route of trace skew in point of intersection, though avoid intersection point to be considered to belong to and multiplely unit to occur and the situation of Voronoi figure cannot be generated, but due to the migration processing of intersection point, change intersection point and belong to multiple this fact that unit occurs, bring certain division error.In addition, these two kinds of methods, owing to being adopt method line feature being separated into point set at equal intervals, can affect, and efficiency are lower to original factor shape.

Summary of the invention

For this reason, the present invention, mainly for crossing wire geographic element, proposes a kind of crossing wire geographic element space partitioning method based on intersection auxiliary point.On the one hand, directly adopt the point of original line key element to generate Voronoi figure as discrete point, the shape of original factor can not be changed, and execution efficiency is higher; On the other hand, when occurring that line feature intersects situation, by the cross-point locations both sides respectively in two line features respectively equal proportion increase by two trace separately but vision is difficult to the auxiliary point distinguished, when not changing intersection point and belonging to this fact of multiple generation unit, avoid the situation that cannot generate Voronoi figure, and there is good subdivision effect.

Main thought of the present invention is: according to the principle of Voronoi apportion model, and will intersect Linear element discrete is point set; And the crossover location to two line features intersected, carry out the process of increasing point; Carry out Delaunay triangulation based on point set again and generate Voronoi figure; Cutting is carried out to the Voronoi figure generated, finally the subregion that the point set belonging to same line feature generates is carried out merger, obtain different line feature subregion separately, thus realize the uniformly subdivision intersecting wire geographic element.

Concrete steps of the present invention are as follows:

Step one: load GIS line feature set L={l _i| i=1,2 ..., n}, n are the number of line feature.

Step 2: ask for different line feature l in line feature set L _i, l _jbetween point of crossing C _ij(x, y), and carry out interpolation process, detailed process is:

(a) first, for line feature l _i, read line feature l _ion point set, search line feature l _itwo points that overcrossing point is contiguous, are designated as C respectively _a(x _a, y _a), C _b(x _b, y _b);

B (), according to following formula (1), calculates and is inserted into an I _c(x _c, y _c), I _d(x _d, y _d) coordinate, wherein (x, y) is point of crossing C _ijcoordinate, (x _a, y _a) and (x _b, y _b) be the C of two near point of crossing _a, C _bcoordinate, (x _c, y _c) and (x _d, y _d) for being inserted into an I _c, I _dcoordinate, s is step factor, s desirable 10,100,1000 ... 10 ⁿ.

$\left\{\begin{array}{c}\mathrm{\Δd}=\frac{\sqrt{{(y_b-y_a)}^2+{(x_b-x_a)}^2}}{s}\\ k=\frac{y_b-y_a}{x_b-x_a}\\ x_c=x+\mathrm{\Δd}\\ y_c=k\×x_c+y_a-k\×x_a\\ x_d=x-\mathrm{\Δd}\\ y_d=k\×x_d+y_a-k\×x_a\end{array}\right.---\left(1\right)$

C () will put I _a, I _bbe inserted into line feature l _iin some C _a, C _bbetween;

(d) secondly, for another line feature l at this place, point of crossing _j, repeated execution of steps (a)-(c), is inserted into line feature l by the corresponding insertion point generated _jin;

Step 3: by discrete for line feature set L be point set P={p _i(x _i, y _i) | i=1,2 ..., k}, k are the number (comprising intersection auxiliary point) of a key element, and use the line feature that " line feature " attribute of point set P belongs to record this point;

Step 4: based on point set P, generate Delaunay triangulation network D, the detailed process of algorithm is as follows:

A () introduces 3 auxiliary point p ₁₁(x ₁₁, y ₁₁), p ₂₂(x ₂₂, y ₂₂), p ₃₃(x ₃₃, y ₃₃), build an initial convex hull S, auxiliary point coordinate can be tried to achieve by formula (2), and wherein MaxX/MaxY is maximum x/y coordinate figure in point set P, AvgX/AvgY is the mean value of x/y coordinate in point set P, and MinX/MinY is minimum x/y coordinate figure in point set P.This initial convex hull S must comprise in point set P institute a little;

$\left\{\begin{array}{c}{.x}_{11}=\mathrm{MaxX}+\mathrm{MaxX}-\mathrm{AvgX}\\ y_{11}=\mathrm{AvgY}\\ x_{22}=\mathrm{MinX}\\ y_{22}=\mathrm{MaxY}+\mathrm{MaxY}-\mathrm{MinY}\\ x_{33}=\mathrm{MinX}-\frac{\mathrm{MaxX}-\mathrm{AvgX}}{2}\\ y_{33}=\mathrm{MinY}-\mathrm{MaxX}+\mathrm{MinY}\end{array}\right.---\left(2\right)$

B () finds out range points p from point set P ₁2 lines are become the initial baseline of Delaunay triangulation network, follow LOP principle, then find out thirdly by nearest point, build first triangle.

C (), based on first leg-of-mutton another two baseline, is followed local optimization procedure principle, is found their the 3rd points respectively.The all newly-generated limits of circular treatment, follow local optimization procedure principle successively, find thirdly, until a little all become the point of Delaunay triangulation network;

Step 5: based on Delaunay triangulation network D, generates Voronoi and schemes V;

Step 6: by " line feature " property value of each point in point set P, is assigned to the corresponding Thicssen polygon t comprising this point _i" former line feature " attribute;

Step 7: Voronoi schemes the cutting process of V.By solving given rectangular extent R (the minimum rectangle scope as figure layer data) and each Thicssen polygon t in figure V _icommon factor, generate the Thicssen polygon set t after cutting _i'.Each Thicssen polygon t _ithe concrete steps of cutting are:

A () obtains rectangular extent R and Thicssen polygon t _iall intersection points, and deposit in point set D;

B () obtains the Thicssen polygon t be included in rectangular extent R respectively _isummit and be included in Thicssen polygon t _ithe summit of interior rectangular extent R, continues to add in point set D;

(c) by point set D first put into point set D again, form end to end point set D;

D (), by point set D by the sequence of clockwise or counterclockwise order, the polygon be made up of point set D is the Thicssen polygon t after cutting _i';

E the above step of () cycling, until complete the polygonal cutting process of all Thicssen;

Step 8: according to " former line feature " attribute, will have the Thicssen polygon t of same " former line feature " property value successively in Voronoi figure _i' merge into a polygon, each polygonal region is the coverage of a line feature.

Compared to prior art, outstanding advantages of the present invention is mainly manifested in three aspects: one is that the present invention directly adopts the point of original line key element to generate Voronoi figure as discrete point, can not change the shape of original line key element; Two is when line feature intersects, the present invention by the cross-point locations both sides in two line features respectively equal proportion increase by two trace separately but vision is difficult to the auxiliary point distinguished, do not change intersection point and belong to this fact of multiple generation unit, avoid the situation that cannot generate Voronoi figure; Three is that automaticity of the present invention is high, and prior art still needs to realize by ArcGIS platform.

Accompanying drawing explanation

The process flow diagram of Fig. 1 the inventive method;

Fig. 2 experiment line figure layer data;

Line chart layer data after Fig. 3 point set;

Fig. 4 increases a process schematic diagram;

Fig. 5 increases the point set line chart layer data after a process;

Fig. 6 is the Delaunay triangulation network D generated;

Fig. 7 is that the Voronoi generated schemes V;

Fig. 8 Voronoi is schemed V carry out cutting after figure V';

The design sketch of Fig. 9 line chart layer data uniformly subdivision.

Embodiment

Below in conjunction with accompanying drawing also by description example for the uniformly subdivision of wire geographic element, further illustrate effect of the present invention.This example selects the line chart layer of a shp form to be experimental data.

Specific implementation process is as follows:

1, open the line chart layer file of layer " road " by name, as shown in Figure 2.Stored in line feature set L, gather L in this example and comprise 4 line features, be respectively l ₁, l ₂, l ₃, l ₄.Its ID attribute is respectively 1,2,3,4.

2, ask for the point of crossing C in L between different line feature, and carry out interpolation process, the point set of 4 line features wherein in this example as shown in Figure 3.There are 4 point of crossing in example, are respectively a C ₁(-96.07,0.33), some C ₂(-24.39,5.74), some C ₃(-75.71 ,-51.78) and some C ₄(-0.39 ,-42.68).For each point of crossing, carry out interpolation process, obtain vicinity two points of their place line features according to formula (1) respectively, in the present embodiment, get s=100, try to achieve 16 interpolated points.With point of crossing C ₁for example, based on point of crossing C1, interpolation goes out l respectively ₁on two points (-93.36851,0.336116), (-93.56850,0.32959) and l ₃on two points (-93.36851,0.11350), (-93.56850,0.55220), interpolation results is as shown in Figure 4.All the other C ₂around four interpolated points are (-24.29142,5.74238), (-24.49142,5.74008), (-24.29142,5.57623), (-24.49142,5.90623); C ₃four interpolated points are around (-75.60735 ,-51.76809), (-75.80734 ,-51.79939), (-75.60735 ,-52.28374), (-75.80734 ,-51.28374); C ₄four interpolated points are around (-0.29368 ,-42.68370), (-0.49368 ,-42.68370), (-0.29368 ,-42.97601), (-0.49368,-42.39139) layer after the increasing point of, last, generation as shown in Figure 5.

3, by discrete for line feature set L be point set P, in this example, point set P={ (-156.29 ,-8.18), (-128.03,-2.23), (-94.98,0.28) ..., (9.38,-77.19) } have 50 points, comprising 16 interpolated points, and the ID attribute of 4 line features is assigned to " line feature " attribute of its corresponding discrete point.

4, based on an elements combination P, generate Delaunay triangulation network D.

1) 3 auxiliary point: p are added ₁₁, p ₂₂, p ₃₃, form initial convex hull S.In the present embodiment, can calculate three point coordinate by formula (2), be p respectively ₁₁(185.48 ,-23.86), p ₂₂(-156.29,195.26), p ₃₃(-214.44 ,-249.10).

2) from point set P, range points p is found out ₁2 lines are become the initial baseline of Delaunay triangulation network by nearest point.In the present embodiment, some p ₂(-128.03 ,-2.23) range points p ₁recently.Follow local optimization procedure principle, then find out thirdly, this example mid point p ₁₀(-19.47,5.80) are the 3rd point, some p ₁, some p ₂with a p ₁₀build first triangle.

3) based on first leg-of-mutton another two baseline, follow local optimization procedure principle, find their the 3rd points respectively.The all newly-generated limits of circular treatment, follow local optimization procedure principle successively, find thirdly, until a little all become the point of Delaunay triangulation network; In this example, generate 101 triangles altogether based on point set P, and leave in container listTriangle.Wherein 11 triangles generated based on auxiliary point do not show, and Fig. 6 shows 90 triangles not comprising auxiliary point.

5, based on Delaunay triangulation network D, generate Voronoi and scheme V.

After generating Voronoi figure, in the present embodiment, by " line feature " attribute of every bit in point set P, be assigned to the corresponding Thicssen polygon t comprising this point _i" former line feature " attribute.As with a p ₁for example, its " line feature " property value is 1, and the Thicssen polygon of its correspondence is t ₁, then t ₁" former line feature " property value be 1.

Generate 53 Thicssen polygons in the present embodiment altogether, according to these 53 groups of data, generate the line chart layer file of layer " V figure " by name, as shown in Figure 7.

6, Voronoi schemes the cutting process of V.

Using the minimum rectangle scope of mileage chart layer data as R scope in the present embodiment, R={ (-156.29,47.14), (69.17,47.14), (69.17 ,-100.98), (-156.29,-100.98), (-156.29,47.14) }, the dot cycle of 50 in point set P is proceeded as follows successively (with a p ₁for example):

1) p is read ₁thicssen polygon t ₁, analyze and obtain t ₁be overlapping relation with the spatial relation of R.

2) rectangular domain R and Thicssen polygon t is obtained ₁intersection point, totally 2 points, be respectively point (-156.29,37.89) and point (-156.29 ,-30.95) and deposit in point set C ₁in.

3) t be included in rectangular domain R is obtained ₁summit and be included in t ₁in the summit of rectangular domain R, continue to obtain point (-156.29 ,-30.95) and point (-137.09 ,-29.31) and continue to deposit in point set C ₁in.

4) by C ₁in first point (-156.29,37.89) deposit in C1, form the equal point set of ending.

5) after point set C being pressed sequence counterclockwise, obtain point set C={ (-156.29 ,-30.95), (-156.29,37.89), (-150.12,32.62), (-137.09,-29.31), (-156.29 ,-30.95) } totally 5 points, form ending and to connect new region t ₁'.

The above step of cycling, the new polygon t of each generation _i' (represent p _ithe new region that point is corresponding).V figure after process, as shown in Figure 8.

7, according to " former line feature " attribute of the Thicssen polygon set T' after cutting, merge Thicssen polygon t' in Voronoi figure.

Carry out the Voronoi after cutting process and scheme V, be made up of 50 Thicssen polygon t'.According to " former line feature " attribute of 50 Thicssen polygon t', the Thicssen polygon t' deriving from same line feature is merged.In the present embodiment, rectangular domain R is just divided into 4 sub regions, and in corresponding the present embodiment, ID attribute is 1,2 respectively, the coverage of the line feature of 3,4.As shown in Figure 9.

Claims (4)

1., for a space partitioning method for crossing wire geographic element, it is characterized in that, concrete steps are as follows:

Step one: load GIS line feature set L={l _i| i=1,2 ..., n}, n are the number of line feature;

Step 2: ask for different line feature l in line feature set L _iand l _jbetween point of crossing C _ij(x, y), and carry out interpolation process;

Step 3: by discrete for line feature set L be point set P={p _i(x _i, y _i) | i=1,2 ..., k}, k are the number of a key element, wherein comprise the number of intersection auxiliary point, and use the line feature that " line feature " attribute of point set P belongs to record this point;

Step 4: based on point set P, generates Delaunay triangulation network D;

Step 5: based on Delaunay triangulation network D, generates Voronoi and schemes V;

Step 6: by " line feature " property value of each point in point set P, is assigned to the corresponding Thicssen polygon t comprising this point _ithe property value of " former line feature ";

Step 7: scheme V to Voronoi and carry out cutting process, by solving given rectangular extent R and each Thicssen polygon t in figure V _icommon factor, generate the Thicssen polygon set T after cutting _i';

Step 8: according to the Thicssen polygon set T after cutting _i' " former line feature " attribute, successively Voronoi is schemed the Thicssen polygon t in V with identical " former line feature " property value _i' merge into a polygon, each polygonal region is the coverage of a line feature, thus realizes the uniformly subdivision intersecting geographic element.

2. a kind of space partitioning method for crossing wire geographic element according to claim 1, it is characterized in that, the detailed process of described step 2 is:

(a) first, for line feature l _i, read line feature l _ion point set, search line feature l _itwo points that overcrossing point is contiguous, are designated as C respectively _a(x _a, y _a), C _b(x _b, y _b);

B (), according to following formula (1), calculates and is inserted into an I _c(x _c, y _c), I _d(x _d, y _d) coordinate, wherein (x, y) is point of crossing C _ijcoordinate, (x _a, y _a) and (x _b, y _b) be the C of two near point of crossing _a, C _bcoordinate, (x _c, y _c) and (x _d, y _d) for being inserted into an I _c, I _dcoordinate, s is step factor, s desirable 10,100,1000 ... 10 ⁿ,

$\left\{\begin{array}{c}\mathrm{\Δd}=\frac{\sqrt{{(y_b-y_a)}^2+{(x_b-x_a)}^2}}{s}\\ k=\frac{y_b-y_a}{x_b-x_a}\\ x_c=x+\mathrm{\Δd}\\ y_c=k\×x_c+y_a-k\×x_a\\ x_d=x-\mathrm{\Δd}\\ y_d=k\×x_d+y_a-k\×x_a\end{array}\right.---\left(1\right)$

C () will put I _aand I _bbe inserted into line feature l _iin some C _a, C _bbetween;

(d) secondly, for another line feature l at this place, point of crossing _j, repeated execution of steps (a)-(c), is inserted into line feature l by the corresponding insertion point generated _jin.

3. a kind of space partitioning method for crossing wire geographic element according to claim 1, is characterized in that, in described step 4, the detailed process generating Delaunay triangulation network D is as follows:

A () introduces 3 auxiliary point p ₁₁(x ₁₁, y ₁₁), p ₂₂(x ₂₂, y ₂₂), p ₃₃(x ₃₃, y ₃₃), build an initial convex hull S, auxiliary point coordinate can be tried to achieve by formula (2), wherein MaxX/MaxY is maximum x/y coordinate figure in point set P, AvgX/AvgY is the mean value of x/y coordinate in point set P, MinX/MinY is minimum x/y coordinate figure in point set P, this initial convex hull S must comprise in point set P institute a little;

$\left\{\begin{array}{c}{x}_{11}=\mathrm{MaxX}+\mathrm{MaxX}-\mathrm{AvgX}\\ y_{11}=\mathrm{AvgY}\\ x_{22}=\mathrm{MinX}\\ y_{22}=\mathrm{MaxY}+\mathrm{MaxY}-\mathrm{MinY}\\ x_{33}=\mathrm{MinX}-\frac{\mathrm{MaxX}-\mathrm{AvgX}}{2}\\ y_{33}=\mathrm{MinY}-\mathrm{MaxX}+\mathrm{MinY}\end{array}\right.---\left(2\right)$

B () finds out range points p from point set P ₁2 lines are become the initial baseline of Delaunay triangulation network, follow the principle of local optimization procedure, then find out thirdly by nearest point, build first triangle;

C (), based on first leg-of-mutton another two baseline, is followed local optimization procedure principle, is found their the 3rd points respectively, the all newly-generated limits of circular treatment successively, follow local optimization procedure principle, find thirdly, until a little all become the point of Delaunay triangulation network.

4. a kind of space partitioning method for crossing wire geographic element according to claim 1, is characterized in that, in described step 6, to each Thicssen polygon t _ithe concrete steps of cutting are:

A () obtains rectangular extent R and Thicssen polygon t _iall intersection points, and deposit in point set D;

B () obtains the Thicssen polygon t be included in rectangular extent R respectively _isummit and be included in Thicssen polygon t _ithe summit of interior rectangular extent R, continues to add in point set D;

(c) by point set D first put into point set D again, form end to end point set D;

D (), by point set D by the sequence of clockwise or counterclockwise order, the polygon be made up of point set D is the Thicssen polygon t after cutting _i';

E the above step of () cycling, until complete the polygonal cutting process of all Thicssen.