
Geometry : 
Subject: Re: [geometry] Support for geographic coordinate system
From: Barend Gehrels (barend_at_[hidden])
Date: 20141130 14:26:21
Hi Adam,
Adam Wulkiewicz wrote On 28112014 16:34:
> Hi Barend,
>
> Barend Gehrels wrote:
>> 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.
>
> Ok, so an intersection of a plane containing 2 endpoints and a
> spheroid's center is indeed an ellipse. Even though it's fragment is
> not the shortest path AFAIU it's mapping into a sphere is commonly
> used because it's a great simplification. Still AFAIU an
> approximation, but probably good enough.
>
>>
>>>
>>> (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.
>>
>
> Ok, AFAIU our discussion is about various mappings of the coordinates
> into a sphere.
Hmm, was that the central point? I thought it is not about mapping
coordinates into a sphere (why should we?) but side calculations.
What I said above is indeed not true, xyz coordinates are different
(true) but the plane will be the same, as you point out below. So
indeed, no need to change the SSF (as basically was also stated on my
blog). Sorry for the confusion.
> Because at this point, when we had the sphere we'd be able to handle
> the arc as a part of a great circle, so e.g. use SSF. In the case of
> side calculation the radius of a sphere is not important, only angles.
That is right.
> Furthermore the longitude would be the same on a sphere so only the
> latitude matters. On the wikipedia page mentioned by you above there
> are 3 methods of mapping mentioned.
>
>
>
>
>
> When you mentioned about a great ellipse I assumed that you'd like to
> use the 2nd one ? (theta), so to use latitude being an angle
> corresponding to the line between the center and the point (the green
> one on the picture above). So to project an elliptic arc radially into
> the surface of a sphere. Did I missunderstood you?
>
> Anyway, I wanted to preserve the latitude ? (fi) (use plane
> perpendicular to the surface, the red one) which is the same as the
> 3rd mapping (the blue one). I felt that the arc would be closer to the
> shortest one. The difference is that they're using radius = a and I
> wanted to calculate it from the endpoints to better represent the
> original arc, but in the case of side calculation on a sphere the
> radius is not important.
>
> So AFAIU the last mapping preserves lon and lat, so it's basically the
> same as the original SSF.
Right, that is indeed the intention.
> And since it states (on the wiki) that:
>
> The last method gives an easy way to generate a succession of
> waypoints on the great ellipse connecting two known points A and B.
> Solve for the great circle between (\phi_1,\lambda_1) and
> (\phi_2,\lambda_2) and find the waypoints on the great circle. These
> map into waypoints on the corresponding great ellipse.
>
> I think we don't need to do anything additional with the SSF. We could
> think about a totally different formula but we probably won't improve
> the SSF. Still I don't know if we need to improve it.
>
>
> But, I feel that we could play with the distance calculation. I.e. use
> one of the mappings, calculate the position and the radius of a
> corresponding sphere and then use haversine. Actually I'm thinking
> about the 3rd mapping (or rather the red one). So we could take the
> coordinates, calculate their position in 3D, then calculate the 3D
> plane perpendicular to the surface. Then somehow calculate the mid
> point. Using those 3 points in 3D calculate the radius of a
> circle/sphere containing all 3 points and the position of the center.
> Then the angle between the endpoints in 3D. Then the length of the arc
> would be the distance. Assuming that the above could be done I'd say
> that it'd be a method corresponding to SSF. This would also be an
> approximation but most certainly faster than Vincenty distance. And
> most certainly good enough in some situations.
Do I understand that this might be an alternative distance calculation,
faster than Vincenty and still better than Haversine?
I did not really plan or intend to change distance calculations.
However, if we have a distance calculation over the spheroid which is
compatable with the SSF over the spheroid, that would be welcome of
course, I've the feeling that could be of use.
>
>>
>>>>>
>>>>> 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".
>>
>
> Ok, so AFAIU if we used it, it'd be an approximation. I'm not saying
> that a bad approximation. So ok, I think we agree.
OK good ;)
>
>>>
>>> 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.
>
> Ok, for now I don't have any other ideas ;)
I did not have the time yet to dive deep into this. I'm not sure if you
mean to check distance or side calculations, or both...
The side calculations (likewise, segmentsegment intersections) are
crucial to all the overlay algorithms (union, intersection,
(sym)difference, buffer) we want to implement, as is also
segment/segment intersections. They should be mutually consistent, we
should use them together. If distance belong to them, I don't know yet,
that is to be figured out. Even if SSF over the spheroid is a little
(but only a little) different, I mean does not follow the geodesic, it
can possibly be still useable.
Regards, Barend
Geometry list run by mateusz at loskot.net