Geometry Boost List: Re: [geometry] Support for geographic coordinate system

Subject: Re: [geometry] Support for geographic coordinate system
From: Barend Gehrels (barend_at_[hidden])
Date: 2014-11-28 04:55:23

Hi Adam,

[snipped]

Adam Wulkiewicz wrote On 28-11-2014 3:26:
>
>
> I'm not sure if there is something like "great elliptic arc" or "great 
> ellipse". AFAIU shortest paths on spheroid aren't fragments of spheres 
> or ellipses. 

http://en.wikipedia.org/wiki/Great_ellipse

> According to wikipedia 
> (http://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid), geodesics 
> which supposed to be the shortest paths between points, in general 
> aren't closed.

That is right.

>
> (From wikipedia)
> If the Earth is treated as a sphere, the geodesics are great circles 
> (all of which are closed) and the problems reduce to ones in spherical 
> trigonometry. However, Newton (1687) showed that the effect of the 
> rotation of the Earth results in its resembling a slightly oblate 
> ellipsoid and, in this case, the equator and the meridians are the 
> only closed geodesics. Furthermore, the shortest path between two 
> points on the equator does not necessarily run along the equator.
>
> Though I'm not really sure if the side of this geodesic can be checked 
> comparing azimuths calculated using Vincenty formula. And how precise 
> it is.
>
>>
>>> In other words, it would be true if all intermediate points of a 
>>> segment was going through the same coordinates. Are you sure that 
>>> this is the case? I'm not. My test seems to prove that this is not 
>>> true or at least that a method (Vincenty) which is known for giving 
>>> precise results for an ellipsoid gives different results than SSF. 
>>> Of course this is more visible for greater flattening.
>>
>> OK. Vincenty is precise but very slow (maybe you have to compare 
>> performance too). At least distance calculations and comparisons I 
>> did in the past, comparing with Andoyer, Vincenty was much slower.
>>
>> If we are going to use side-information on ellipsoids we have to have 
>> a reasonably fast algorithm, the side calculations are used very very 
>> often. The SSF is already slower of course than the cartesian 
>> calculation.
>>
>> If the difference is very small, and the results of SSF-only are 
>> mutually consistent, and we don't use it together with distances (I 
>> don't think we do - but we use fractions which have to be comparable) 
>> it is not impossible that we can still use it (needs more research). 
>> It is just a different method to calculate the side, and if we can 
>> use that for turn-calculations and point-in-polygon calculations, it 
>> might give the correct results. Other libraries or packages first 
>> project to a Cartesian coordinate system, do the operation and 
>> convert back - you will also loose some information there too.
>> I don't state that it is possible, we just can research the 
>> possibilities.
>>
>> Furthermore the formulas given (on that blog) and probably also what 
>> is implemented for SSF are about spheres. For the ellipsoid x,y,z 
>> have to be calculated differently, taking axis lengths into account. 
>> That might (most probably will) result in different results. It has 
>> to be implemented. That is the phase where we are.
>>
>
> Hmm, do you mean by this that we could place the plane with more 
> precision on a spheroid, as perpendicular to the surface in 2 
> endpoints? Which maybe could result in a plane oriented differently? 
> But still use the SSF for calculation, so assuming that the segment is 
> a fragment of a sphere?

That segment is a fragment of an ellipse, the cut of a plane with the 
spheroid is an ellipse.
But for the rest, sure, I mean that. We should do that.

>
> So the difference would be the way how the orientation of the plane 
> was calculated. It might as well give the same result because he 
> definition of geographical coordinates is the same as for spherical - 
> latitude is an arc between XY plane and a line perpendicular to the 
> surface. So it's possible that this plane would have the same 
> orientation. But it's definietly something worth checking.

But the xyz coordinates are different so the plane will be different.

>>>
>>>  Do you have an idea how we could verify this?
>>>
>>
>> No... or basically it is probably already clear now. The methods are 
>> different.
>
> I came out with somethign like this. Please say if I'm missing 
> something or if it's completly wrong.
>
> A hypothesis:
>
> If we gathered points of geodesic (coordinates where left side changes 
> to right side between the segment's endpoints), form a path, smooth it 
> e.g. by running bg::simplify(), and then calculate the length by 
> suming all distances between the path's points, we should be able to 
> check which curve is the shortest when measured on spheroid and 
> sphere. This way we could check which one better represents the 
> shortest path on both models.

It is known that the path over a great ellipse is not the shortest path. 
See in "Geometrical Geodesy: Using Information and Computer Technology" 
(google gives a preview, search for great elliptic arc", quote: "The 
great Elliptic arc is not the shortest possible connection between two 
points".

>
> Furthermore if a hypothesis about a different coordinates (shape) of 
> geodesic on spheroid and sphere was true, when calculating the length 
> of geographical geodesic on spheroid we should get smaller result than 
> when calculating the length of spherical geodesic. And the other way 
> around spherical geodesic should be shorter on sphere.
>
> And finally if we also took a naiive cartesian-like segment between 
> endpoints in the coordinates space, e.g. the result of a cartesian 
> side_by_triangle. The length of this segment "projected" on the 
> surface of a sphere or spheroid should be different than the "true" 
> distance between endpoints calculated using haversine or vincenty. 
> What is more, it should be longer than the spherical segment on a 
> sphere and geographical segment on a spheroid. It's because the 
> shortest path e.g. on a sphere is defined by a great circle so any 
> other path should be longer.
>
> I've implemented this in test-geoside and this is what I've got:
>
> All results are for a segment (-31 -31, 31 31).
> distance is the result of a bg::distance(s1, s2, strategy)
> length is the sum of distances between points of a geodesic/curve of 
> left/right edge points.
> length geo is the length of a curve of side change according to 
> azimuths calculated with Vincenty formula
> length sph is the length of a curve of side change according to SSF
> length car is the length of a curve of side change according to 
> side_by_triangle (also verified algebraically)
> A letter in brackets is a strategy used to calculate distances
> (V) - vincenty, (A) - andoyer : so assuming that points are on spheroid
> (H) - haversine mean radius = (2a+b)/3, (M) - haversine max radius = a 
> : assuming that points are on sphere
>
> 1. a=1 b=1 - unit sphere
>
> distance (V) = 1.4910384559056726
> distance (A) = 1.4910384559056726
> distance (H) = 1.4910384559056729
> distance (M) = 1.4910384559056729
> length geo (V) = 1.4904379394565843
> length sph (V) = 1.4904379394565843
> length car (V) = 1.4897397736057592 // the shortest?
> length geo (A) = 1.4904379394565841
> length sph (A) = 1.4904379394565841
> length car (A) = 1.4897397736057592 // ?
> length geo (H) = 1.4904379394565841
> length sph (H) = 1.4904379394565841
> length car (H) = 1.4897397736057592 // ?
> length geo (M) = 1.4904379394565841
> length sph (M) = 1.4904379394565841
> length car (M) = 1.4897397736057592 // ?
>
> All geo and sph distances are the same, ok.
>
> But the length of cartesian segment is the smallest. It's shorter than 
> the distance calculated using directly. It's quite late so it's 
> possible that I'm missing something. But it seems that either my 
> assumptions are bad or the method is wrong or not precise enough.
>
> But let's go forward...
>
> 2. a=1 b=1-1/300 - unit Earth
>
> distance (V) = 1.4867431728537253
> distance (A) = 1.4867328035365481
> distance (H) = 1.4893817465102219
> distance (M) = 1.4910384559056729
> length geo (V) = 1.4856765746677782 // smaller than sph, ok
> length sph (V) = 1.4861482692673844
> length car (V) = 1.4854692836326653 // the shortest again?
> length geo (A) = 1.485666275142842  // also smaller than sph
> length sph (A) = 1.4861379754374746
> length car (A) = 1.485459072084804  // ?
> length geo (H) = 1.4883117296448467 // but this shouldn't be smaller 
> than sph
> length sph (H) = 1.4887818973016325
> length car (H) = 1.4880813951276908 // ?
> length geo (M) = 1.4899672488101912 // also this shouldn't be smaller
> length sph (M) = 1.4904379394565841
> length car (M) = 1.4897366580811147 // ?
>
> So the length of geographical segment seems to be smaller no matter 
> how distances are calculated. This is rather impossible since the 
> shortest path on a sphere is a great arc fragment. Ok, maybe the 
> method's error is too big in this case. Maybe the flattening is too 
> small to get meaningful results.
>
> Cartesian the shortest again.
>
> 3. a=1 b=0.9 - unit Saturn
>
> distance (V) = 1.3716008213860109
> distance (A) = 1.3618688848319354 // now the difference between A and 
> V is greater
> distance (H) = 1.4413371740421503
> distance (M) = 1.4910384559056729
> length geo (V) = 1.3705358784007811 // smaller than sph, ok
> length sph (V) = 1.3711112408432529
> length car (V) = 1.3709916809430027 // also this is greater than geo, 
> but the closest to distance
> length geo (A) = 1.3608547786450773 // the smallest!
> length sph (A) = 1.3614390188832899
> length car (A) = 1.3614090781917894
> length geo (H) = 1.4403203633519746 // damn, this shouldn't be smaller 
> than sph
> length sph (H) = 1.4407566748080314
> length car (H) = 1.4400787694784105 // and this is the smallest now
> length geo (M) = 1.4899865827779046
> length sph (M) = 1.4904379394565841
> length car (M) = 1.4897366580811147 // the smallest
>
> So again the geo segment seems to be smaller than sph in all cases, 
> which is AFAIU wrong.
>
> Cartesian still confusing. E.g. the smallest for Vincenty but the 
> length is the closest to directly calculated distance.
>
> 3. a=1 b=0.75
>
> distance (V) = 1.2321397094677933
> distance (A) = 1.1681145282213292 // note the difference between V and A
> distance (H) = 1.3667852512468668
> distance (M) = 1.4910384559056729
> length geo (V) = 1.2294812729423881 // the smallest
> length sph (V) = 1.2316859225469647
> length car (V) = 1.2319096796713431 // finally, the greatest but the 
> closest to the distance
> length geo (A) = 1.165467037015419  // here ok too
> length sph (A) = 1.1679406380233488
> length car (A) = 1.168887231887985
> length geo (H) = 1.3649878832500875 // but this should be greater than sph
> length sph (H) = 1.3662347778352022
> length car (H) = 1.3655723783015628 // and this too should be greater 
> than sph
> length geo (M) = 1.4890776908182768 // same as above
> length sph (M) = 1.4904379394565841
> length car (M) = 1.4897153217835226
>
> Cartesian confusing again.
>
> 4. a=1 b=0.5
>
> distance (V) = 1.1102958385138535
> distance (A) = 0.8451906005369858 // andoyer can be ignored I guess
> distance (H) = 1.2425320465880607
> distance (M) = 1.4910384559056729
> length geo (V) = 1.1064353530360909 // smaller than sph, but the 
> furthest from distance
> length sph (V) = 1.1146694840023121
> length car (V) = 1.1087899062397581 // this is the closest
> length geo (A) = 0.83893092464406038 // smaller than sph
> length sph (A) = 0.84544333659011361
> length car (A) = 0.84807843141588624
> length geo (H) = 1.2432889622300805 // greater than sph, still may be 
> a coincidence
> length sph (H) = 1.2420316162138203
> length car (H) = 1.241438095270917 // the smallest
> length geo (M) = 1.4919467546760965 // greater than sph
> length sph (M) = 1.4904379394565841
> length car (M) = 1.4897257143251006
>
> Similar story, if the length is the smallest it's the furthest from 
> the distance calculated directly.
>
>
> So I'm guessing that the length and the distance depends mainly on the 
> algorithm used for distance calculation. The variance in the result of 
> calculation using various distance strategies seems to be so big that 
> the possible error based on the method (gathering pixels + 
> simplification) influences the result a lot more than the shape of a 
> path. Hence, we wouldn't check anything like this.
>
> Or do you have some ideas what may be happening?

Not yet.

> The assumption about the length of geodesics should at least be good 
> for sphere.
> Any help would be welcome.
>
>
> Vincenty distance is always shorter than haversine which I guess is good.
> Should Vincenty give good results for different flattening than Earth's?
> Should Andoyer give good results only for Earth?

Might be, it is AFAIK designed for Earth...

> It's less precise approximation but it uses flattening and radius but 
> still gives worse results than haversine. Is this a bug or a 
> limitation to flattenings near Earth's?
>
> I also noticed that the greater the flattening, the longer it takes to 
> use Vincenty formula. Maybe we should calculate the threshold based on 
> the flattening?
>
> Furthermore at flattening 0.8 the formula falls into an infinite loop 
> in which lambda and previous_lambda (also sigma) swapps between 2 
> values all of the time but the difference is greater than the 
> threshold. So we could add a test for this effect. So add some test 
> for local minimum. And maybe some iterations counter just in case. If 
> you're curious, it's for (in radians):
>
> lon1 = -0.54105206811824214
> lat1 = -0.54105206811824214
> lon2 =  0.89709923552509452
> lat2 = -1.0646508437165410
> a    =  1.0000000000000000
> b    =  0.20000000000000001

Thanks a lot. This requires more study, will come back later.

Regards, Barend