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

Subject: Re: [geometry] Support for geographic coordinate system
From: Adam Wulkiewicz (adam.wulkiewicz_at_[hidden])
Date: 2014-11-27 21:26:34

Hi Barend,

Barend Gehrels wrote:
> Adam Wulkiewicz wrote On 27-11-2014 1:40:
>> 2014-11-26 23:22 GMT+01:00 Barend Gehrels <barend_at_[hidden] 
>> <mailto:barend_at_[hidden]>>:
>>
>>
>>     Adam Wulkiewicz wrote On 26-11-2014 21:25:
>>>
>>>     Of course in all formulas some approximation of a globe is used.
>>>     My assessment is based on a fact that SSF assume thet a globe is
>>>     a sphere 
>>
>>     Ah, so you did *not *read my blog, and apparently also not my
>>     statements that SSF can be used on a spheroid.
>>
>>
>> Of course I *did* read it.
>
> OK, good, then can we now discuss the spheroid properties of the SSF.
>
>>     So I repeat it again, there it goes:
>>
>>     http://barendgehrels.blogspot.nl/2011/06/spherical-side-formula.html
>>
>>     Quoting, a.o.:
>>     "Summary: Using some high-school mathematics I presented an
>>     algorithm and a formula to calculate at which side a point is
>>     with resepect to a segment, on a sphere or on the Earth"
>>
>>     Beside the summary, it has a whole paragraph about the Ellipsoid.
>>
>>     Because this is exact, and Vincenty is an approximation (close to
>>     reality, but still, an approximation), I now just assume that
>>     these results are correct and the Vincenty approximation is off
>>     within ranges of 0.6.
>>
>>     I assume that, until it is proven that it is the other way round...
>>
>>
>> The paragraph about the Earth, corresponding planes, etc. would be 
>> true only if a shape of spherical geodesic and ellipsoidal one was 
>> the same.
>
> I don't write about geodesic there. What is calculated is a plane 
> through the center of the Earth, and the two points involved. What is 
> determined is if it is left or right from the plane. That is exact, 
> also for a spheroid.
>

Ok, but do I understand correctly that since the plane is intersecting a 
sphere all intermediate points are lying on an arc on sphere's surface, 
not aspheroid surface, so a circle in 3D space?

> So what is written there is basically true, but... indeed the plane 
> going through the center of the Earth does not encompass the geodesic. 
> So if we want to have it comparable with the geodesic shortest 
> distance calculations, it cannot be used together.
>
> What is defined there is the "great elliptic arc" or "great ellipse" 
> (comparable to great circle, but on a spheroid) and  "The great 
> ellipse is not the shortest path between two points on the Earth's 
> surface".
>

AFAIU when SSF is used, a segment is not really an arc being a fragment 
of the intersection of a plane atached to the center of an ellipsoid and 
2 endpoints. It's an arc formed by coordinates projected into a sphere. 
It isn't an arc being an intersection of ellipsoid's surface and a plane 
perpendicular to the surface in 2 endpoints and intersecting the center 
of the CS. Still even if the shape of the arc in 3D is different than a 
circle fragment, it those may still be the same coordinates, I don't know.

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. 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.

(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?

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.

Anyway, I don't have anything against using SSF by default. I've done 
some testing (below). Unfortunately they didn't give a clear answer 
which method is better.

>
>>
>>     If SSF was only applicable for 100% spheres, the error would be
>>     much larger. Compare Haversine/Vincenty, it can be off many
>>     kilometers.
>>
>>
>> Yes, it's not surprising that the distance may be more influenced. 
>> I'm guessing that in the case of distance more important is the 
>> difference between major an minor axes than the actual shape of the 
>> shortest path. Still it doesn't mean that a side is calculated correctly.
>>
>>     Note: it is different from distance, where an exact mathematical
>>     formula for ellips or spheroid just does not exist.
>>
>>
>> AFAIU it's the same for ellipsoid. And greater the flattening, 
>> greater the difference. See the picture for flattening = 0.5.
>>
>>>     and Vincenty that it's a spheroid. So the approximation is
>>>     closer to reality. And the greater the flattening, the greater
>>>     the difference between SSF and Vincenty. I'm not considering
>>>     precision, numerical errors, stability etc.
>>>     As for the Vincenty formula here:
>>>     http://www.icsm.gov.au/gda/gdav2.3.pdf states:
>>>     "Vincenty's formulae (Vincenty, 1975) may be used for lines
>>>     ranging from a few cm to nearly 20,000 km, with millimetre
>>>     accuracy."
>>>     Calculation done with SSF is probably also precise but it uses a
>>>     sphere model, that's all.
>>
>>     Nope.
>>
>>
>>  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.

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?
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? 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

Regards,
Adam