Geodesic algorithms: an experimental study
This is a modal window.
The media could not be loaded, either because the server or network failed or because the format is not supported.
Formal Metadata
| Title |
| |
| Title of Series | ||
| Number of Parts | 295 | |
| Author | ||
| Contributors | ||
| License | CC Attribution 3.0 Germany: You are free to use, adapt and copy, distribute and transmit the work or content in adapted or unchanged form for any legal purpose as long as the work is attributed to the author in the manner specified by the author or licensor. | |
| Identifiers | 10.5446/43353 (DOI) | |
| Publisher | ||
| Release Date | ||
| Language |
Content Metadata
| Subject Area | ||
| Genre | ||
| Abstract |
| |
| Keywords |
FOSS4G Bucharest 2019230 / 295
15
20
28
32
37
38
39
40
41
42
43
44
46
48
52
54
57
69
72
75
83
85
87
88
101
103
105
106
108
111
114
119
122
123
126
129
130
131
132
137
139
140
141
142
143
144
147
148
149
155
157
159
163
166
170
171
179
189
191
192
193
194
195
196
197
202
207
212
213
214
215
216
231
235
251
252
263
287
00:00
Coordinate systemCurvatureProjektive GeometrieAlgorithmGeodesicBitDistanceEllipsoidSurface of revolutionEndliche ModelltheorieWorkstation <Musikinstrument>Web pagePhysical systemBenchmarkComputational geometryLecture/Conference
00:47
MultiplicationSphereGeodesicGoodness of fitCoordinate systemState of matterSurface of revolutionSpacetimeInformationEndliche ModelltheorieHyperplaneAnalogyCurveObject (grammar)Rhumb lineAngleAlpha (investment)NavigationCommutatorDistanceComputational geometryLevel (video gaming)SphereApproximationsalgorithmusSeries expansionComplex (psychology)Numerical integrationPoint (geometry)Trigonometric functionsInverse elementWell-formed formulaPythagorean theoremProjektive GeometrieLine (geometry)Great circleEllipsoidINTEGRALImage resolutionIterationSet (mathematics)AzimuthOrder (biology)BenchmarkSeries (mathematics)Observational studyGeodesicAlgorithmPresentation of a groupQuicksortService (economics)Video projectorForcing (mathematics)Multiplication signCASE <Informatik>Arithmetic meanCondition numberProduct (business)Valuation (algebra)Regulator geneThermal expansionMeasurementOvalFault-tolerant systemBridging (networking)PhysicalismLecture/Conference
06:17
AlgorithmGeodesicGeometryIterationSeries (mathematics)ApproximationSeries (mathematics)Numerical analysisDisintegrationEllipseInverse elementCartesian productLocal ringAxonometric projectionPoint (geometry)ApproximationAreaPopulation densitySimilarity (geometry)Hausdorff spaceRevision controlBenchmarkEllipsoidWikiLengthWell-formed formulaCASE <Informatik>Goodness of fitCodeProjektive GeometrieAlgorithmSoftwareSeries (mathematics)Line (geometry)Point (geometry)Computational geometryResultantDistanceBenchmarkGeometryImplementationEllipseSimilarity (geometry)ApproximationsalgorithmusBlock (periodic table)String (computer science)BuildingInterpolationDivergenceGeodesicNumerical analysisInverse elementAzimuthAreaImage resolutionFlow separationDirection (geometry)State of matterMeasurementOpen sourceSeries expansionInverse problem1 (number)Arc (geometry)Black boxOrder (biology)Limit (category theory)Thermal expansionForm (programming)TelecommunicationNetwork topologyPolygonSet (mathematics)TrailMedical imagingService (economics)Visualization (computer graphics)Open setGraphics tabletSummierbarkeitProcess (computing)JSONXML
13:25
GeodesicAverageTable (information)GeometryCurvatureSphereAzimuthDistanceHopf algebraApproximationSeries (mathematics)Point (geometry)Ordinary differential equationAreaCartesian productAxonometric projectionEquals signOperations researchMilitary operationGroup actionCylinder (geometry)Conic sectionPlot (narrative)OracleScaling (geometry)Execution unitGoodness of fitCASE <Informatik>Automatic differentiationRevision controlFile formatGrass (card game)Newton's law of universal gravitationBitProjektive GeometrieTable (information)AreaResultantDistancePoint (geometry)Connected spaceDifferent (Kate Ryan album)Error messageEqualiser (mathematics)Computational geometryBenchmarkMetreParameter (computer programming)Optical disc driveSoftware testingInsertion lossDenial-of-service attackLogic gateSaddle pointDreizehnMultiplication signProcess (computing)Limit (category theory)Digital photographyCartesian coordinate systemSet (mathematics)AzimuthWell-formed formulaFlow separationGeodesicEllipsoidAcoustic shadowTrapezoidCurvatureApproximationsalgorithmusWikiHierarchySeries (mathematics)Order (biology)GeometryAlgorithmLibrary (computing)AuthorizationInverse problemImplementationComputer animation
20:33
Basis <Mathematik>Game theoryDifferent (Kate Ryan album)AreaBitPoint (geometry)CollisionMetreSquare numberError messageOperator (mathematics)Computational geometrySeries (mathematics)Computer animationLecture/Conference
21:47
Mach's principleAzimuthDemosceneDistancePoint (geometry)Process (computing)Roundness (object)Workstation <Musikinstrument>Multiplication signWeb pageMachine visionState of matterPresentation of a groupLecture/Conference
Transcript: English(auto-generated)
00:08
All right, good afternoon, everyone. Thank you for coming to the session on geodesic algorithms. I'd like to welcome our next speaker to present his research. Thank you. Thanks.
00:21
So today we're going to talk a bit about how we compute distances on ellipsoid revolution and some benchmarks. So let me start with some basics. So there are several models of the Earth and coordinate system.
00:41
So the simplest one, if you get some projection to the two dimensions, is to model the Earth as flat and then have some bad approximation in general, but good if you are in some local place. And then a better level of approximation
01:01
is to consider the Earth as a sphere. And then you could have log-latitude coordinates. As far as I understand, the state of the art for GIS computations is the ellipsoid of a revolution.
01:20
So the Earth is not exactly a sphere. And then you could also have geodes. So I'm not a space analyst in this, so maybe you have geophysics here can tell us more. So you have more information about the altitude also and valleys and mountains. So I'm not considering this model here,
01:44
at least in this presentation. So what is a geodesic? It's the shortest path between two points. So in 2D, it's just a straight line. In the sphere, it's a great circle. So it's an intersection of a hyperplane passing
02:03
through the center of the sphere. In ellipsoid, it's not the same. So if you intersect the ellipsoid with a hyperplane passing from the center, you will get a great ellipse, which is analog to great circle, which is not the shortest path. So the shortest path is a curve
02:21
that goes around your object. So since there are a few talks here about Mercator projections, it's nice to see that, of course, there the straight lines are not shortest paths.
02:43
They're called loxodromes that keep the azimuth, which is the angle from the point that you are to the north or south pole. And of course, as many of you know better than me, the shortest paths are curves on this projection.
03:06
So in general, there are two main problems, two main geodesic problems. The one is called the direct and the other the inverse. So in the first, we are given the point, the azimuth, and the distance. And we won't compute the new point at that distance.
03:22
So here, we are given a, and then lowercase alpha 1, which is the azimuth, and the distance, s12. And we want to compute b. So it's like navigating. We are in some point of an angle
03:41
and a value of the distance. And then we're going to go there. So this is a direct geodesic problem. And we also have the inverse, actually the distance commutation problem. So we have two points, a and b. And we won't compute s12. And also as a byproduct, we can compute the azimuth alpha 1.
04:06
So first, sorry for this slide. This is just to illustrate the complexities between those two models. I will not explain the formulas. The first is self-explanatory.
04:20
It's the Pythagoras theorem in 2D. There is an ancient closed formula. For the sphere, the shortest distance between two points is a trigonometric function. It's readily evaluated. But on ellipsoid, the shortest distance
04:41
is given by a set of integrals. So you have to do some more complicated computations. This is not rocket science. It's something that is solved. But here we like to study benchmarks on performance and accuracy.
05:01
So this is really old, of course. But compared to the sphere, this is not a closed formula. So what you usually do is some evaluation, either by numerical integration or series expansion. And then you also need some iteration
05:24
to convert to the solution. So no more technical details here. So this is a set of, I tried to collect a set of algorithms for geodesic problems. So I'm open to suggestions for more if you have.
05:42
The first one is a series approximation, which goes back to Bessel. And it's rediscovered, maybe re-engineered by Carney, together with an iterative method to give a very accurate solution up to series order
06:05
8 to the two geodesic problems that I described before. Another method is a numerical integration, how we studied the literature.
06:20
So we are not considering this here. And the third one is, as far as I see, the most widely used. It's within the iterative method for computing the distance. And there are some old school ones that can be retrieved from papers from the 60s and 70s.
06:44
There are series approximation of order one and two by Thomas and Doyer. And then you have the spherical one, which is the Haberstein formula that we showed. And there are some, again, in this series of papers
07:00
from the 60s, there are also some inverse elliptic arc length formulas. Again, those are by series expansions. And the last thing that you can do is to take some projection that maybe fits on your data and do some Cartesian computation.
07:23
So the first case is implemented in geographic leap. We still have a work from boost geometry. So we have a pull request from a good number of code project that we, in the future, we also
07:41
add this algorithm in boost geometry. So two is not, I don't even know if it's open source or I don't know if any implementation is available. And the rest, we support them in one or the other way
08:02
in boost geometry. So I'm just saying this to declare which software I will use for the benchmark. So it's basically geographically bent implementations in boost geometry. OK, so this is the first example. I want to compute the distance between Athens and Brussels. These are given in longitude latitude.
08:23
So the first method of this high series expansions give this value for this particular problem up to this accuracy. It visually seems the same as Vincenti. And then you can see that by choosing
08:41
the other less accurate methods, we have some divergence in the results. So Thomas is a few centimeters. I assume that the first one is the most accurate. And then I have Andoia and these elliptic formulas.
09:02
OK, so apart from Vincenti, all these are these 70s old school formulas that we use in boost geometry because they are fast. But as we can see, they also give some nice results. OK, I will systematically explain this later.
09:24
And then we have also a spherical approach. And here if we project and compute the distance, so this is one kind of projection, which I am not citing here, you get something really bad
09:46
because the distance is too long. OK, so distance computation and the direct problem are building blocks for other algorithms too.
10:01
So if you want to compute the distance between geometries, in general, so for example, between a line string and a multi-polygon, then the building block will be to compute the distance between two points or a point of the segment. So it is important to build on top
10:22
all the or many geographic algorithms. So when you want to compute the area of a polygon, then you need to compute the azimuth. So again, you need to solve the inverse problem. In other algorithms too, like a line interpolate points
10:43
when you interpolate points on a line string or if you densify a geometry, then inside you want to solve both direct and inverse problems. So they are almost everywhere in the whole set
11:02
of geographic algorithms. Also the last one, it's the nearest case of algorithms. So if you want, for example, you have tracks, GPS tracks that is given by line strings and you
11:22
want to compute the similarity, there are some similarity measures like Hausdorff or Frézé distance. Again, you have to compute your black box will be distance computation. So you have to solve the inverse problem. So OK, one more example here.
11:43
So this is a building block for distance between geometries. So I have a line string, so this is just a segment. And I have a given point, and I want to compute the distance. So this will solve, at least in both geometry,
12:01
how it's implemented will solve several inverse and direct geodesic problems. And yes, bi-spherical have a very simple formula. I'm getting this result. And I can use Andoyer, Thomas, and Vinceti algorithms to get those results.
12:22
And then out of curiosity, we can also use a PostGIS, which is the state of the art and also gives some most widely used and also give two versions, like spherical, which is the same as in both geometry. So we are OK here.
12:41
But geographic gives a different number. So yes, I will leave it to you to think of how is it possible to check if these results are correct or not, and which result is more accurate. Because it's not visually, from the first sight,
13:07
possible to do it. So I'll continue with a benchmark. So that was some illustrative examples. So we use a data set, which is open in Zenodo.
13:22
It contains 100,000 geodesics computed by uniform distributed points on the ellipsoid. And it also contains 50,000 short geodesics that are less than 10 kilometers.
13:45
And using this data set, I will just run tests with the previous described algorithms, just to see which is faster, which is more accurate. So regarding accuracy, this data set
14:02
also provides an accurate solution for the direct and the inverse problem. So this is computed in geographic lib by the author of geographic lib, with some exact arithmetic,
14:21
so no floating points, and supposed to be very exact. So we can use it to test the accuracy. And there is also a wiki for more details. So looking at performance for distance, OK, flat.
14:40
It's very fast, but as we saw, even for simple examples, it does not give any good results. Spherical, which is the first approximation that you can find, is 10 times faster than Andoyer, which is our first working geographic approach.
15:02
And then this is almost 10 times faster than geographically boosted, which is the most accurate. So more or less, we have a hierarchy of performances. So the stars are using the implementations in boost geometry, and the other is geographically.
15:25
So then I want to test accuracy. These are some remarks. So for long geodesics, the faster the method, the lower the accuracy. This is something expected, and we wanted to show it also in benchmarks.
15:42
So for long geodesics, geodesics can almost cover the globe, but are not covering the whole globe. I mean, you don't have two antipodal points, because then you have issues. These are the occurrences that we get from our algorithm.
16:01
So like a spherical, you can have errors in kilometers, and you go up to nanometers for the high series approximation. So Andoyer is the one that in boost geometry we use by default, has a few hundred meters
16:25
accuracy in that long geodesics. And flat approximation has extremely low accuracy, but it's only suitable for very, very specific applications when your data sets are only inside one city.
16:49
Yes, the only method that actually works for points that are antipodal is the method with high series approximation.
17:00
So everything else is broken when the two points are nearly antipodal on the globe. So the result is the formulas are going really, really wrong. So this graph, it's log scale, so those two
17:21
are supposed to be very inaccurate. And you can see the scaling of the algorithms while the distance increases. The first is almost one meter, so it's very short. The last cases are the near-adipodal points. So you can see that apart from a geographic leap, which
17:43
are those versions here that are stable, all the other algorithms at some point, they explode and produce very wrong results. We use several versions of geographic leap. You can compile it with a smaller order series.
18:05
And yes, it's interesting that even with four order series, we get very accurate results. OK, but even three is better than Vincente, so OK. So then we run some benchmarks with Azimuth and ARIA.
18:24
So these are a bit connected, because ARIA is using Azimuth. But yes, and here, the difference is that with very short geodesics, we have very large errors. And also here, we can test some equal area projections.
18:43
OK, let's see the equal area projections. OK, this is accuracy. This is accuracy of Azimuth. Then we can see that in short geodesics, the computed Azimuth is really off. You have to go to around 100 meters to maybe 10 meters
19:03
to get something which is accurate. Geographically, of course, it's performing, and Vincente, relatively good. OK, but still, geographically, it has large errors in small geodesics.
19:20
And then ARIA, it's a bit more interesting. So here, we also test those guys here that are projections, equal area projections, for small distances. So actually, we compute the area of the geodesic and the shadow on the equatorial.
19:42
So we have something like a trapezoid. This is the areas we compute. And yes, you can see that in small distances, projection makes sense. But then here, there is an explosion in the area. So it's really nonsense, the results that you get.
20:03
Again, we have this hierarchy of geographically, which is the most accurate. Then Vincente, and again, this time is over. Yes, and then Vincente, and then Thomas, and Andoja.
20:25
OK, and that's all. Thank you.
20:41
Question on the area calculation, what was the reference point for calculating the error there for the area calculations? What do you mean, the reference point? So you were measuring in accuracy on the y-axis, no? Yes, it's square meters.
21:00
The errors are in square meters. Yeah, compared to what? What's the basis for that? Again, you have an exact computation of the area with the geographic lib, and this exact arithmetic operations.
21:25
So what was the performance of the high series compared to Vincente? Sorry? Performance of the high series method compared to Vincente method? What's the difference? Yes, could you go back to the slide? In performance, OK.
21:43
Yes, it's a bit more than two times, so it's not that. And this is in distance. In azimuth, it's the same. So actually, what you do in Vincente,
22:02
intuitively, is that somehow you kill the process early. But even if you leave it iterate if the points are on the portal, then there is a strong issue. So geographically, by taking care of this,
22:21
probably introduce some performance. Any other last questions before we thank him? No? All right, thank you. One last round of applause. Unfortunately, the third speaker is not here today.
22:41
So if anyone was here for that presentation, our sincere apologies. But other than that, I guess the session is closed. So thank you, everyone, for attending. And enjoy the rest of your time here. See you at the AGM. Thank you.