You're right. I was right thinking about it yesterday while committing the
changes.
*f *is misleading, and it should be substituded by *rf* in my formula.
Thanks for having reported your measurements. They helped a lot to figure
out the bug ;)
giovanni
2012/4/6 blavet didier.bla...@ird.fr
**
*Good evening*
Of course thank you again i fully agree with what you found which is
mathematically correct ,
it was not easy to find, and is essential for debugging the calculation
of ellipsoîdal distances
I just had a quick look at the code found at
http://trac.osgeo.org/qgis/browser/trunk/qgis/src/core/qgsdistancearea.cpp#L540
i'm not an informatician so i can just say few words
inverse flattening seems correctly defined in line 148 as
mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor )
therefater the bug that you found could be in line 153 (and may be in
some other lines of the code, i cannot say)
Therafter, it is just a detail (and a convention) but
just to try to help to avoid a potential confusion in the future code
maintenance. ,
about the letter f used in your message about inverse flattening i
think that it would may be better to use another letter than f for
designation of inverse flattening ... I say this just because in line
550 of the code found today at
http://trac.osgeo.org/qgis/browser/trunk/qgis/src/core/qgsdistancearea.cpp#L540the
letter f was already assigned as f = 1/mInvFlattening, indicating that
letter f in this code is already reserved for ... flattening)
So of course and again what you said is perfecly correct from a
mathematical point of view (as long as you say, by convention, that f =
inverse flattening) ,
Therefore the demonstration would may be a bit clearer just with
replacing this letter f by some other symbol or combination of symbols for
the designation of inverse flattening
(as this f is already used for something else in the code).
Would it be possible for instance to use the symbolic combination *rf*,
used in proj4 (for reciprocal flattening which is, as far as i know, in
english, equivalent to inverse flattening)
then this would give, from what you found this morning , and with exactly
the same meaning that what you said :
*In Qgis the semi-minor axis of the ellipsoid is calculated with:
b = a - (rf/a)
where
b = semi-minor axis
a = semi-major axis
rf = inverse flattening (= reciprocal flattening = a/(a-b) )
while it should be:
b = a - (a/rf)*
therefore, this rf symbol would also fit with the +rf proj4 parameter
defined as
+rf reciprocal of the ellipsoid flattening term (e.g. 298)
(also in proj4 +f is used for flattening
+f Flattening of the ellipsoid (often presented as an inverse,
e.g. 1/298)
(by the way , I guess, however, that the adjective reciprocal for
reciprocal numbers can have also some other meanings for the pure
mathematicians ..
... Therefore it would also be possible to use another symbol for inverse
flattening than rf (but preferentially not f (:-))
*Thank you again
Best regards
*
Le 06/04/2012 14:15, G. Allegri a écrit :
Not at all.
The change has already been committed to the development version.
giovanni
2012/4/6 blavet didier.bla...@ird.fr
Thank you very much for this very fast and efficient diagnosis
Best regards
Didier Blavet
Le 06/04/2012 12:01, G. Allegri a écrit :
Probably I've found the bug.
In Qgis the semi-minor axis of the ellipsoid is calculated with:
b = a - (f/a)
where
b = semi-minor axis
a = semi-majot axis
f = inverse flattening
while it should be:
b = a - (a/f)
In Qgis the WGS84 semi-minor axis is 6378136,xxx
while it should be 6356752.xxx
giovanni
2012/4/6 G. Allegri gioha...@gmail.com
AFAICS, the Qgis formula [1] to calculate distances on ellipssoid is
Vincenty's.
I've reproduced it for the first set of points and I can confirm tha
QGis gives 540.9973919726176 m
The online Vincenty calculator [2] gives 540.045 m
I suppose it's a rounding issue.
I will investigate it more...
giovanni
[1]
http://trac.osgeo.org/qgis/browser/trunk/qgis/src/core/qgsdistancearea.cpp#L540
[2] http://www.movable-type.co.uk/scripts/latlong-vincenty.html
2012/4/6 blavet didier.bla...@ird.fr
Dear QGis Developpers and Users
___
*Didier B**lavet*
Institut de Recherche pour le Développement - IRD (http://www.ird.fr)
UMR EcoSols - Ecologie Fonctionnelle Biogéochimie des Sols
Agroécosystèmes
(Montpellier SupAgro-CIRAD-INRA-IRD) (
http://www.montpellier.inra.fr/ecosols)
Bâtiment 12, 2 place Viala, F-34060 Montpellier Cedex 2 , France
Tel bureau : +33 (0)4 99 61 21 33
Secrétariat : +33 (0)4 99 61 21 01 / Fax : +33 (0)4 99 61 21 19
Courriel : didier.bla...@ird.fr
___
