
Geometry : 
Subject: Re: [geometry] Support for geographic coordinate system
From: Barend Gehrels (barend_at_[hidden])
Date: 20141128 04:55:23
Hi Adam,
[snipped]
Adam Wulkiewicz wrote On 28112014 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 sideinformation 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 SSFonly 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 turncalculations and pointinpolygon 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 cartesianlike 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 testgeoside 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=11/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
Geometry list run by mateusz at loskot.net