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#L540 the 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 <[email protected]>
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 <[email protected]>
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

2012/4/6 blavet <[email protected]>
Dear QGis Developpers and Users


_______________________________________________________________________________________

Didier Blavet

Institut de Recherche pour le Développement - IRD (http://www.ird.fr)

UMR Eco&Sols - 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 : [email protected]

_______________________________________________________________________________________


_______________________________________________
Qgis-user mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/qgis-user




--

_______________________________________________________________________________________

Didier Blavet

Institut de Recherche pour le Développement - IRD (http://www.ird.fr)

UMR Eco&Sols - 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 : [email protected]

_______________________________________________________________________________________



--

_______________________________________________________________________________________

Didier Blavet

Institut de Recherche pour le Développement - IRD (http://www.ird.fr)

UMR Eco&Sols - 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 : [email protected]

_______________________________________________________________________________________

_______________________________________________
Qgis-user mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/qgis-user

Reply via email to