CN110334311A - The area computation method of Subgraph spot on a kind of earth ellipsoid face - Google Patents
The area computation method of Subgraph spot on a kind of earth ellipsoid face Download PDFInfo
- Publication number
- CN110334311A CN110334311A CN201910587514.4A CN201910587514A CN110334311A CN 110334311 A CN110334311 A CN 110334311A CN 201910587514 A CN201910587514 A CN 201910587514A CN 110334311 A CN110334311 A CN 110334311A
- Authority
- CN
- China
- Prior art keywords
- area
- spot
- ellipsoid
- trapezoidal
- calculation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 47
- 238000004364 calculation method Methods 0.000 claims abstract description 63
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 13
- 230000006870 function Effects 0.000 description 22
- 230000008569 process Effects 0.000 description 7
- 230000001186 cumulative effect Effects 0.000 description 5
- 238000011835 investigation Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 230000006872 improvement Effects 0.000 description 4
- 239000000203 mixture Substances 0.000 description 3
- PCTMTFRHKVHKIS-BMFZQQSSSA-N (1s,3r,4e,6e,8e,10e,12e,14e,16e,18s,19r,20r,21s,25r,27r,30r,31r,33s,35r,37s,38r)-3-[(2r,3s,4s,5s,6r)-4-amino-3,5-dihydroxy-6-methyloxan-2-yl]oxy-19,25,27,30,31,33,35,37-octahydroxy-18,20,21-trimethyl-23-oxo-22,39-dioxabicyclo[33.3.1]nonatriaconta-4,6,8,10 Chemical compound C1C=C2C[C@@H](OS(O)(=O)=O)CC[C@]2(C)[C@@H]2[C@@H]1[C@@H]1CC[C@H]([C@H](C)CCCC(C)C)[C@@]1(C)CC2.O[C@H]1[C@@H](N)[C@H](O)[C@@H](C)O[C@H]1O[C@H]1/C=C/C=C/C=C/C=C/C=C/C=C/C=C/[C@H](C)[C@@H](O)[C@@H](C)[C@H](C)OC(=O)C[C@H](O)C[C@H](O)CC[C@@H](O)[C@H](O)C[C@H](O)C[C@](O)(C[C@H](O)[C@H]2C(O)=O)O[C@H]2C1 PCTMTFRHKVHKIS-BMFZQQSSSA-N 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 239000000284 extract Substances 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 241001339245 Callirhoe digitata Species 0.000 description 1
- 235000002259 Callirhoe involucrata Nutrition 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000001983 electron spin resonance imaging Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000005580 one pot reaction Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
- G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
- G06F16/29—Geographical information databases
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F8/00—Arrangements for software engineering
- G06F8/30—Creation or generation of source code
- G06F8/31—Programming languages or programming paradigms
- G06F8/315—Object-oriented languages
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
- G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Databases & Information Systems (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Algebra (AREA)
- Pure & Applied Mathematics (AREA)
- Remote Sensing (AREA)
- Computing Systems (AREA)
- Geometry (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The area computation method of Subgraph spot on a kind of earth ellipsoid face of the present invention, algorithm thinking divides top layer, middle layer and bottom: bottom is substantially the trapezoidal areal calculation of ellipsoid, this method is by decimal number variable and function, break ellipsoid trapezoid area precise calculation computational accuracy limitation, the area value for being up to 23 effective digitals is obtained, ensure that the areal calculation precision of any figure spot；It is that each ellipsoid is trapezoidal greatly that top layer, which decomposes any figure spot, is added up by the algebraical sum of major trapezoidal area and obtains any figure spot area；Middle layer then passes through the area that recursive call underlying algorithm obtains big trapezoidal segment.Method presents the entire flows of any figure spot areal calculation, and are verified by example, and polygon clockwise is made to obtain positive area, and polygon counterclockwise obtains negative product.This method calculates thinking and is illustrated, and figure spot area precision controlling is simple and reliable, is easy to software realization and popularization.
Description
Technical field
The invention belongs to earth ellipsoid face any figure spot areal calculation fields, specifically, in order to avoid Gauss-Ke Lvge
The influence of distortion of projection, in the more demanding application of area precision, when needing to calculate the earth ellipsoid face area of soil patch,
Use the inventive method.
Background technique
Subgraph spot areal calculation document is broadly divided into two classes on earth ellipsoid face: one kind is directly calculated based on coordinate,
Such as, area (Shi Yimin, Zhu Ziyang the geodetic that (2006) such as a people calculate convex polygon on ellipsoid using geodetic coordinates is applied
Coordinate calculates algorithm [J] Tongji University's journal (natural science edition) of the convex area of a polygon of ellipsoid, 2006,34 (04):
504-507), woods green equal (2007) utilizes the method for curve surface integral to the meter of region area on ellipsoid by rectangular space coordinate
Calculation method studied (woods is green, on Ma Jingsong earth ellipsoid face region area algorithm research [J] survey and draw notification, 2007,
(06): 8-10)；The trapezoidal indirect calculating (state affairs on another kind of to be based on two meridians on ellipsoid and two parallel circles be boundary
Second national land investigation leading group, institute office map sheet Theoretical Area and figure spot ellipsoid area calculation formula and requirement
[Z] .2008-3-28), which can substantially be divided into three levels: top layer, middle layer and bottom.Top layer provides entirely
The calculating thinking of figure spot polygon, it may be assumed that each side of polygon ABCD and any given warp L0Surround a big ladder diagram
(see Fig. 1, by taking the side AB as an example, A, B point two o'clock are along weft direction in warp L for block0On subpoint be respectively A1、B1, then quadrangle
ABB1A1The as side AB and L0The big trapezoidal segment surrounded), to all big trapezoidal its algebraical sums of segment areal calculation surrounded, just
Obtain the area of polygon ABCD.Middle layer gives the calculation method of single big trapezoidal segment: single big trapezoidal segment
ABB1A1It is cut by weft, splits into many small trapezoidal segment AE1F1A1, calculate its area S1, S1It is cumulative just to obtain ladder diagram
Block ABB1A1Area.Bottom gives the calculation method of small trapezoidal segment: being converted to two meridians and two parallel circles is boundary
The trapezoidal calculating of ellipsoid, calculate its area using approximate formula.Two class area computation methods compare, and the latter is easier to
Understand and applies.However, the second class method, does not provide the big trapezoidal details specifically split, should specifically split into how many
It is a small trapezoidal? or it is small trapezoidal high in how many suitable after splitting? history keeps just equal (2018) and has carried out exploration to the problem
Research, studies have shown that can be substituted using the thinking of Approximate calculation of definite integral as the height for splitting small trapezoidal segment is smaller and smaller
Approximate formula carries out the areal calculation of big trapezoidal segment, and computational accuracy gradually increases, and obtains trapezoidal high-precision of sample
Spend reliable value 2661732.9601161m2(the improvement of figure spot area calculation method on Shi Shouzheng, Shi Yishao, Zhao Wei ellipsoid
[J] Wuhan University Journal information science version, 2018,43 (5): 779-785).
Technical problem: calculating any ellipsoid figure spot area and need to calculate big trapezoidal segment area, and big trapezoidal segment face
Product computational accuracy is split dependent on small trapezoidal near infinite and small trapezoidal areal calculation.The small trapezoidal area meter of ellipsoid
Calculation is substantially trapezoidal come what is completed using ellipsoid, and ellipsoid trapezoid area formula has precise calculation and approximate calculation public
Point of formula, under the premise of using independent variable and return value being the system function of double-precision floating points, using approximate calculation public affairs
Formula computational accuracy is higher.However, double-precision floating points only have 15-16 effective digitals, for single ellipsoid trapezoid area meter
Precision is fine for calculation, near infinite split it is small it is trapezoidal for, the precision may influence its it is cumulative after big trapezoidal faces
Product precision.So far, the areal calculation of Subgraph spot is limited only to the theory study stage on ellipsoid, domestic mainstream
GIS software SuperMap and international well-known GIS software ArcGIS are without providing the calculating function of available figure spot ellipsoid area
Energy.
For this purpose, I integrates the knowledge that existing ellipsoid figure spot area calculates indirectly, a kind of earth ellipsoid face is invented and has taken up an official post
It is intended to the area computation method of spot, and by the case verification reliability of this method, it, can if GIS software uses this method
To be precisely calculated the earth ellipsoid face area of any figure spot.
Summary of the invention
The invention calculates any figure spot area using ellipsoid is trapezoidal indirectly, and calculation method can substantially be divided into three layers
It is secondary: top layer, middle layer and bottom.
Top layer thinking is consistent with document, and (office, second national land investigation leading group, State Council map sheet is theoretical
Area and figure spot ellipsoid area calculation formula and require [Z] .2008-3-28), for any figure spot polygon, split into and its
The equal several big trapezoidal segments of number of edges, calculate the area algebraical sum of all big trapezoidal segments to get the area for arriving Subgraph spot.
Middle layer calculates big trapezoidal segment using the method in the specification of number of patent application 201910046051.0
Area.Without loss of generality, with trapezoidal segment ABB big in Fig. 11A1For, illustrate this method.The midpoint for taking the side AB is E ((B1+
B2)/2, (L1+L2)/2), E point is along warp direction in weft B1、B2On subpoint be respectively EA、EB(see Fig. 2), passes through warp
EAEBFormer big trapezoidal segment is cut, with the hatched triangle Δ EBE on the right of the warp cut awayBTo fill up the warp left side
Open triangle Δ EAEA, then ellipsoid segment trapezoidal greatly is transformed to the trapezoidal E of ellipsoidAEBB1A1, the ellipsoid trapezoid area
For the approximate evaluation value of big trapezoidal segment.It changes for a thinking, big trapezoidal segment area true value can be using accurate expression as ellipsoid
Trapezoidal area (approximate evaluation value) subtracts Open triangle area, adds hatched triangle area, it may be assumed that
S segment ABB trapezoidal greatly1A1The trapezoidal E of=SAEBB1A1-SΔEAEA+SΔEBEB
View discovery is reviewed to hatched triangle and Open triangle, the two can be considered as the big trapezoidal special case of ellipsoid,
Its area computation method is big trapezoidal completely the same with conventional ellipsoid.One warp side of hatched triangle and Open triangle can
Be considered as it is big it is trapezoidal in L0Side, the weft of the two in a corresponding big trapezoidal weft, the two bevel edge another
Vertex then corresponds to big another trapezoidal weft side, and it is 0 that only the edge lengths, which are degenerated, in this way, cutting the shade three of bowl spares generation
Angular and Open triangle can be considered as short bottom edge to degenerate being the big trapezoidal of 0 length, accordingly it is also possible to by bevel edge
Benefit is cut, it is trapezoidal (see Fig. 2) to be converted into ellipsoid, is calculated using big trapezoidal segment areal calculation true value expression formula.The algorithm
It is realised in that recurrence, two bottoms, (hatched triangle and Open triangle are as big for the high big trapezoidal area of ellipsoid reduced
Trapezoidal special case, it is trapezoidal greatly to be classified as ellipsoid) the segment areal calculation expression formula progress trapezoidal greatly of recursive call original ellipsoid
It calculates.
Bottom provides the trapezoidal high-resolution method of ellipsoid.There are four the trapezoidal area formula of ellipsoid is main, one
A precise calculation and three approximate formulas of equal value (figure spot area meter on Shi Shouzheng, Shi Yishao, Zhao Wei ellipsoid
Improvement [J] Wuhan University Journal information science version of calculation method, 2018,43 (5): 779-785).All four formula come
Derived from following integral expression formula:
In formula, T is ellipsoid trapezoid area, L1、L2And B1、B2Respectively surround two trapezoidal warp and weft, a
For semimajor axis of ellipsoid (unit: m), b is semiminor axis of ellipsoid (unit: m), and e is the first eccentricity, e2=(a2-b2)/a2。
The integral is directly solved to get precise calculation is arrived:
And power series expansion is first carried out to the integral formula, then carry out subitem integral, it arranges, can be obtained three theories etc.
The approximate formula of valence.The approximate calculation wherein recommended again with office, second national land investigation leading group, State Council
Formula is the most commonly used, specific formula is as follows:
Wherein A, B, C, D, E are constant, are calculated as follows
A=1+ (3/6) e2+(30/80)e4+(35/112)e6+(630/2304)e8
B=(1/6) e2+(15/80)e4+(21/112)e6+(420/2304)e8
C=(3/80) e4+(7/112)e6+(180/2304)e8
D=(1/112) e6+(45/2304)e8
E=(5/2304) e8
In formula, a is semimajor axis of ellipsoid (unit: m), and b is semiminor axis of ellipsoid (unit: m), and e is the first eccentricity, e2=
(a2-b2)/a2, Δ L is trapezoidal segment through poor (unit: radian), (B2-B1) it is the trapezoidal segment meridional difference (unit: radian), Bm=
(B1+B2)/2。
Compared with approximate formula, precise calculation needs additionally using letters such as logarithm operation, extraction of square root operations
Number.Since the variable and return value of the system function of computer default are double-precision floating points, only 15-16 significant figure
Word causes the computational accuracy of precise calculation not as good as approximate formula.
This method looks for another way, and comprehensively using the decimal system variable of up to 28-29 effective digital, improves sinusoidal in formula
The return value precision of all functions such as function, cosine function, natural logrithm function is to decimal number variable level, at this point, ellipsoid
The areal calculation precision of face precise calculation greatly improves, and approximate formula is limited by giving up for the high power item of e,
It is unobvious to calculate error improvement, although approximate formula takes to e10, computational accuracy is also far below precise calculation.
High-precision function algorithm used in area formula is as follows:
Logarithmic function uses power series approximate calculation, specific formula are as follows:
Specific C# code is accomplished by
The function that extracts square root uses Newton iteration method.The square root for calculating positive number a is equivalent to accounting equation f (x)=x2- a=0
Solution, the classic algorithm for seeking Equation Approximate is Newton iteration method, iterative formula are as follows: Xn+1=(Xn+a/Xn)/2。
Newton iteration method principle is explained are as follows: sets the coordinate of certain point on function f (x) as (x0, y0), cross the tangent line side of the point
Journey is y=kx+b, enables y=0, then can calculate the intersection point abscissa x=-b/k of tangent line and x-axis, and b=y0-kx0, y0=
x0 2- a, k=2x0, intersection point abscissa can be x=(x with abbreviation0+a/x0)/2.Remember x1=(x0+a/x0)/2 continue to seek point (x1,
f(x1)) tangent line and x-axis intersection point abscissa x2=(x1+a/x1)/2, it is clear that x2Compare x1Closer to f (x) function curve and x
The intersection point of axis, it may be assumed that x2Compare x1Closer to the solution of Equation f (x)=0.Continue iteration, subduplicate computational accuracy steps up.
Enable x0The specific C# code of=a, the function that extracts square root are accomplished by
Taylor series, formula are used in SIN function calculating are as follows:
Its specific C# code is accomplished by
π constant has used 30 memory ballads after the decimal point of π value: " one temple of mountain top, one pot wine (3.14159), you are happy bitter
Very I (26535), fill a wine cup for and eat (897), wine kill you (932), kill not dead (384), sauntering, you saunter dead (6264), and fan fan is scraped (338), is fanned
Ear eats wine (3279)."
Definition π numerical constant is 3.141592653589793238462643383279M, totally 31 effective digitals.Program tune
Variable-value monitoring is shown during examination: the decimal number pai_math2 value of system actual use is
3.1415926535897932384626433833 totally 29 effective digitals.
The calculation process of new method are as follows:
(1) constructed fuction calTrap (L1, L2, B1, B2), calculate ellipsoid trapezoid area.The function uses ellipsoid ladder
Shape precise calculation:
Two warps and two wefts define that an ellipsoid is trapezoidal, might as well assume the trapezoidal lower-left angle point warp of the ellipsoid
Latitude coordinate is (L1, B1), upper right angle point latitude and longitude coordinates are (L2, B2), then (L2-L1) > 0, (B2-B1) > 0, ellipsoid is trapezoidal
Area formula calculated result is positive.This does not have any problem for calculating single trapezoidal segment area, still, any figure spot area meter
The top layer for calculating thinking needs the area algebraical sum of major trapezoidal segment, i.e. the trapezoidal area value of some in calculating process needs negative
Number, for this purpose, being examined closely again to the area formula.
If (L1, B1) and (L2, B2) be ellipsoid on any two points, as (L1-L2)*(B1-B2When)=0, only one warp
It is trapezoidal that line or weft can not really surround ellipsoid, i.e., ellipsoid trapezoid area is 0 value；As (L1-L2)*(B1-B2) > 0 when,
According to formula, ellipsoid is trapezoidal to obtain positive area value；And work as (L1-L2)*(B1-B2) < 0 when, ellipsoid is trapezoidal to obtain negative face
Product value.
(2) recursive function calPatch (L is constructed0, L1, B1, L2, B2, err), calculate big trapezoidal area.
With four parameter differences in ellipsoid trapezoid area calculating, ellipsoid segment areal calculation trapezoidal greatly needs six
Parameter, by taking Fig. 2 as an example, (B1, L1) it is A point coordinate, (B2, L2) it is B point coordinate, L0For the L in Fig. 20Warp (might as well assume L0Choosing
It selects in big trapezoidal left side, such as: the minimum longitude that figure layer range can be used makees L0), L0The two of warp and the side AB and excessively AB
Parallel circle surrounds ellipsoid ABB trapezoidal greatly1A1, err is precision controlling variable, by square metre as unit of (such as err=1m2)。
To this it is big it is trapezoidal carry out cutting benefits, be converted to that ellipsoid is trapezoidal, ellipsoid trapezoid area after conversion, which is greater than, controls variable err
When, recursive operation is carried out, conversely, terminating recurrence.Obviously, new two trapezoidal warps are respectively L0, (L1+L2)/2, two wefts
Respectively B1, B2, specific recursive function is as follows:
calPatch(L0, L1, B1, L2, B2, err)
{
// this function is used to calculate the area of ellipsoid segment trapezoidal greatly
Sum=calTrap (L0, (L1+L2)/2, B1, B2)；// calculate the big trapezoidal area estimated value of ellipsoid (ellipsoid ladder
Shape area)
If (Math.Abs (sum) < err) // specified when the absolute value of estimated value (the small triangle comprising cutting benefit) is less than
When value err, no longer recurrence returnsum；
Else//conversely, recursive calculation
return sum+calPatch((L1+L2)/2, L1, B1, (L1+L2)/2, (B1+B2)/2, err)+calPatch
((L1+L2)/2, (L1+L2)/2, (B1+B2)/2, L2, B2, err)；// big trapezoidal area true value=big trapezoidal area estimated value+yin
Shadow triangle area-Open triangle area.Wherein, calPatch ((L1+L2)/2, L1, B1, (L1+L2)/2, (B1+B2)/2,
It err) is Open triangle area, value is negative, calPatch ((L1+L2)/2, (L1+L2)/2, (B1+B2)/2, L2, B2, err)
For hatched triangle area, value is positive.
}
(3) sum function calAllPatch (L is constructed0, Ls, Bs), the algebraical sum for each big trapezoidal area of composition that adds up obtains
The ellipsoid area of any figure spot.When calculating, specified L is needed0Warp numerical value, and the latitude and longitude coordinates on any each vertex of figure spot
(Ls, Bs)。
calAllPatch(L0, Ls, Bs)
{
L1=Ls[0]；B1=Bs[0]；First apex coordinate of // figure spot and first big trapezoidal vertex
Sum=0；Any figure spot area of // storage
For (i=1；I < Ls.length；I++) // successively solve the big trapezoidal area of each composition and add up
{
L2=Ls[i]；B2=Bs[i]；// big second trapezoidal apex coordinate
Sum=sum+calPatch (L0, L1, B1, L2, B2, err) and // area accumulation
L1=Ls[i]；B1=Bs[i]；// big first trapezoidal apex coordinate
}
L2=Ls[0]；B2=Bs[0]；// the last one big second trapezoidal vertex, first apex coordinate of figure spot
Sum=sum+calPatch (L0, L1, B1, L2, B2, err) and // area accumulation
Sum=-sum；// ensuring that area of a polygon clockwise is positive, area of a polygon counterclockwise is negative
}
The advantage of new method:
New method provides the area computation method of Subgraph spot on a kind of ellipsoid of simple possible；This method calculating process
It is controlled on ellipsoid with can be convenient it can be readily appreciated that calculation method is easy to computer realization by controlling the size of variable err
The areal calculation precision of any figure spot.In addition, area of a polygon calculated result is positive clockwise, area of a polygon counterclockwise is calculated
As a result it is negative, slightly adjusts, this method can adapt to the calculating of the area of a polygon containing island.
Detailed description of the invention
Fig. 1: arbitrary polygon reference area (second national land investigation leading group, State Council office on ellipsoid
Room .2008)
Fig. 2: big trapezoidal segment reference area on ellipsoid
Fig. 3: Subgraph spot reference area example on ellipsoid
Specific embodiment
Specific implementation is divided into two parts, first part verify ellipsoid trapezoid area precise calculation precision (algorithm
Floor portions), second part provides the example (middle layer and top layer of algorithm) of any figure spot areal calculation.
First part: without loss of generality, using 80 ellipsoidal parameter of Xi'an, in calculated four warps of Shi Shouzheng etc. (2018)
The trapezoidal ellipsoid that latitude point surrounds is computing object, and specific coordinate is shown in Table 1.
The transit square degree of 1 Xi'an of table, 80 ellipsoid, 4 points
It is utilized respectively exact formulas and approximate formula calculates its ellipsoid area, in calculating process, approximate calculation is public
Formula is also all promoted to decimal number precision level (or 10 using self-editing function-28Precision level).Approximate formula is chosen
The product term formula that second national land investigation leading group, State Council is recommended, and calculate to e10, approximate formula
Calculated value is 2661732.9601182273781954774308m2, and the exact formulas calculated value of sample trapezoidal area is
2661732.9601182342785799189m2, from a high position to low level, 14 effective digitals are identical before the two, and and Shi Shouzheng
The calculated value of system function (its function argument and return value be type double precision variable) is used Deng (2018)
2661732.9601182m2It is identical.
From 15 effective digitals, there is difference in exact formulas and approximate formula.To judge the two precision difference, draw
Enter history and keep the positive improvement Rectangular Method formula for waiting (2018) to use, casts aside precise calculation and approximate formula, directly use
The approximate algorithm of definite integral calculates ellipsoid trapezoid area, and the function value in calculating process still uses high-precision self-editing letter
Number.Specific calculated result is shown in Table 2
Loxodrome area of 24 points of table on 80 ellipsoid of Xi'an
Note: it in ellipsoid trapezoid area calculating process, improves Rectangular Method and needs specified Δ B value size, product term and essence
True formula does not need then
Improving Rectangular Method is that reference Approximate calculation of definite integral thought is built-up, cuts small trapezoidal high Δ B by reducing,
Increase big trapezoidal cutting times, steps up computational accuracy.The calculated result phase of less cutting times and higher cutting times
Compare, the identical number of significant figures of the two is the reliable value of calculated result.As Δ B=10-3When (about 206 "), sample is trapezoidal not
It is cut；As Δ B=10-4When (about 20.6 "), sample is trapezoidal to be cut twice, is improved the reliable value that Rectangular Method calculates and is
2661732.9m2；As Δ B=10-5When, improving the reliable value that Rectangular Method calculates is 2661732.960；As Δ B=10-6When, change
The reliable value calculated into Rectangular Method is 2661732.960118m2；As Δ B=10-7When, improving the reliable value that Rectangular Method calculates is
2661732.96011823m2；As Δ B=10-12When, improving the reliable value that Rectangular Method calculates is
2661732.96011823427857992m2；In Δ B from 10-3Narrow down to 10-12During, calculated result gradually increases, and
As Δ B=10-13When, author computer program runtime is greater than 14 hours, and calculated result is compared with Δ B=10-12When reduce instead, no
Meet the increased trend of numerical value again, to find out its cause, it is general because the cumulative small trapezoidal area numerical value of individual event is too small, so as to cause precision damage
It becomes homeless cause.
Comparison improves the calculated result of Rectangular Method it is found that the sample trapezoidal area reliable value that exact formulas obtains is
2661732.9601182342785799m2, possess 23 effective digitals, far superior to the 14 of product term approximate formula are effective
Number.Although product term approximate formula (10) formula only takes to e10, increase e12Item, e14Its precision can be improved in item, but
It is that, really in this way, all parameters needs of approximate formula recalculate, can not show a candle to directly more convenient using exact formulas.
If referring to ESRI ArcGIS software, the default coordinate resolution value 0.0001m of projected coordinate system, corresponding figure spot face
Product, which calculates, to be needed to be accurate to 10-8m2, this correspond to the sample ellipsoid it is trapezoidal for, need to retain 15 effective digitals.But
In the areal calculation of Subgraph spot, when major trapezoidal segment areal calculation, need approximate to call the ellipsoid trapezoidal many times
Area formula with the cumulative area for obtaining big trapezoidal segment, and then adds up and obtains the area of any figure spot.Many times cumulative must
So the loss of precision is brought therefore in the areal calculation of Subgraph spot, to take to e10The ellipsoid trapezoid area approximate calculation of item
Formula can only ensure 14 effective digitals be it is inadequate, promoting precise calculation precision using high-precision function has to 23
Effect number is necessary.
Second part: without loss of generality, using 80 ellipsoidal parameter of Xi'an, in Shi Shouzheng etc. (2018) calculated ellipsoid
Increase No. 5 points other than trapezoidal greatly 2,4 points, forms triangle of three sides not along latitude and longitude direction and verified as example
(see Fig. 3).The specific coordinate values of triangle figure spot are shown in Table 3, if 2,4 two o'clocks are along weft direction in the meridian (L for crossing No. 5 points0
=116 ° 22 ' 00 ") on subpoint be respectively P '2、P’4, then algorithm top layer figure spot triangle P2P4P5Area can be expressed as
SFigure spot triangle P2P4P5=STrapezoidal P2P4P ' 4P ' 2+STrapezoidal P4P5P ' 4+STrapezoidal P5P2P ' 2
Wherein, latter two triangle, it is trapezoidal (Shi Shouzheng, 2018) to be considered as special ellipsoid, calculation method and one
As ellipsoid it is trapezoidal identical.
The transit square degree of 3 Xi'an of table, 80 ellipsoid, 3 points
Algorithm middle layer is big using recursive calculation using the thinking in the specification of number of patent application 201910046051.0
The area of trapezoidal segment.It is different from the trapezoidal high-precision surface product value for directly obtaining 23 effective digitals using formula of ellipsoid, it is ellipse
The big trapezoidal area calculating of spherical surface always needs to be cut into small trapezoidal to improve computational accuracy.In certain accuracy rating, ellipsoid
Big trapezoidal area is with the small trapezoidal high diminution of cutting, and precision improves and (in other words, reduces with the value of err variable, precision mentions
It is high, it is however generally that, err reduces 10 times, and area precision improves 1).On the other hand, with the small trapezoidal high diminution of cutting, journey
Sequence calculation amount increases.This example is calculating figure spot triangle P shown in Fig. 32P4P5Each composition big trapezoidal segment area when, if err
Variable-value 10-4m2, area is -3993066.3359626475022751453142m2；If err variable-value 10-5m2, area
For 3993066.3359626514660603331232m2.If calculated result is using the value that rounds up to 10-8, both for-
3993066.33596265m2.Consider computational efficiency and current mm grade coordinate setting precision, err variable-value 10-4m2, figure spot
Areal calculation precision uses enough.Err variable-value 10-4m2When, this example calculation is each in the process to form big trapezoidal segment area
Calculated result is shown in Table 4.
Subgraph spot reference area example unit on 4 ellipsoid of table: m2
Note: 1 area of a polygon counterclockwise is negative；2:err=0.0001m2
To sum up, The present invention gives the methods of any figure spot areal calculation, and are verified by example, practice
Show that this method is simple and clear, is easy to computer realization.
Claims (3)
1. the area computation method of Subgraph spot on a kind of earth ellipsoid face, it is characterised in that algorithm bottom is become using decimal number
Amount and function break the computational accuracy limitation of ellipsoid trapezoid area precise calculation, obtain it and be up to 23 significant figures
The area value of word ensure that the precision of any figure spot areal calculation from bottom.
2. determining L according to the precise calculation that right 1 requires0Under the premise of warp is on the left of figure spot, if on figure spot sideline
The order of adjacent two o'clock is from low latitudes to high latitude, then the areal calculation result of corresponding big trapezoidal segment is positive number, otherwise is
It is negative.
3. the areal calculation of the big trapezoidal segment required according to right 2 is as a result, polygon figure spot area is negative clockwise, to meter
It calculates result to negate, is adjusted to positive number, to meet routine, similarly, the final result of polygon figure spot area counterclockwise is adjusted to
Negative, to adapt to the areal calculation of the figure spot containing island polygon.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910587514.4A CN110334311A (en) | 2019-07-02 | 2019-07-02 | The area computation method of Subgraph spot on a kind of earth ellipsoid face |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910587514.4A CN110334311A (en) | 2019-07-02 | 2019-07-02 | The area computation method of Subgraph spot on a kind of earth ellipsoid face |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN110334311A true CN110334311A (en) | 2019-10-15 |
Family
ID=68144770
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910587514.4A Pending CN110334311A (en) | 2019-07-02 | 2019-07-02 | The area computation method of Subgraph spot on a kind of earth ellipsoid face |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110334311A (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111028286A (en) * | 2019-12-17 | 2020-04-17 | 中国人民解放军海军大连舰艇学院 | Precision-controllable pattern spot ellipsoid area calculation method |
| CN112131329A (en) * | 2020-09-11 | 2020-12-25 | 北京数字政通科技股份有限公司 | Method for calculating territorial space area based on ellipsoidal trapezoidal subdivision |
-
2019
- 2019-07-02 CN CN201910587514.4A patent/CN110334311A/en active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111028286A (en) * | 2019-12-17 | 2020-04-17 | 中国人民解放军海军大连舰艇学院 | Precision-controllable pattern spot ellipsoid area calculation method |
| CN112131329A (en) * | 2020-09-11 | 2020-12-25 | 北京数字政通科技股份有限公司 | Method for calculating territorial space area based on ellipsoidal trapezoidal subdivision |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Buchin et al. | Area-preserving simplification and schematization of polygonal subdivisions | |
| CN110909093B (en) | Method and device for constructing intelligent landmark control network | |
| CN110334311A (en) | The area computation method of Subgraph spot on a kind of earth ellipsoid face | |
| CN109671155A (en) | Surface mesh method for reconstructing, system and relevant device based on point cloud data | |
| Tabik et al. | High-performance three-horizon composition algorithm for large-scale terrains | |
| CN109085656B (en) | Feature-oriented high-precision gravity graph construction and interpolation method | |
| CN110879063B (en) | Rapid star map identification method based on triangle matching | |
| CN109345450A (en) | A kind of image mosaic method using geographical entity element information | |
| CN110995849B (en) | Transmission and storage method of vector data position information | |
| CN114387419B (en) | Three-dimensional geographic entity coding method and device based on multi-dimensional space subdivision | |
| CN104090945A (en) | Geographic space entity constructing method and system | |
| CN109813255A (en) | The area computation method of big trapezoidal segment on a kind of earth ellipsoid face | |
| CN109000650B (en) | Method and device for determining relative position of ground object | |
| CN112697136B (en) | Quick minimum area star map simulation method | |
| CN111506576B (en) | Land block coding method and device based on regional quadtree | |
| CN114138925A (en) | Location point belonging area retrieval method, device, electronic equipment, medium and product | |
| CN114490902B (en) | Multi-dimensional space self-adaptive subdivision and coding method and system for two-dimensional geographic entity | |
| Li et al. | Development of mobile GIS system for forest resources second-class inventory | |
| Li et al. | Multi-resolution representation of digital terrain models with terrain features preservation | |
| CN109581286A (en) | A kind of two dimension fingerprint distribution indoor positioning method and structure | |
| CN105761291B (en) | Dijkstra's algorithm based on Weight inlays line automatic generation method | |
| CN112031754B (en) | Method, device and equipment for determining seepage boundary | |
| TWI786874B (en) | Method of the digital grid model and system thereof | |
| JP7267380B2 (en) | METHOD, APPARATUS, AND ELECTRONIC DEVICE TO OVERLAY LASER POINT CLOUD ON HIGH-PRECISION MAP | |
| CN113220812B (en) | Data spatialization method and device based on multi-source map platform cross validation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20191015 |