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 <[email protected]> > ** > *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 <[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 >>> >>> >>> [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 <[email protected]> >>> >>>> Dear QGis Developpers and Users >>>> >>>> >>>> >>>> _______________________________________________________________________________________ >>>> >>>> *Didier B**lavet* >>>> >>>> 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 B**lavet* >> >> 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 B**lavet* > > 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
