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

Subject: Re: [geometry] Support for geographic coordinate system
From: Barend Gehrels (barend_at_[hidden])
Date: 2014-11-06 16:31:35

Hi Adam,

Thanks for your proposal!

Adam Wulkiewicz wrote On 6-11-2014 17:06:
> Hi,
>
> I'm planning to begin the implementation of the support for 
> geographical CS and would like to hear your opinion. Esspecially from 
> people already using the geographical extensions.
> Basically, to support this CS we need a distance strategy and side 
> strategy. They're already implemented (more or less) though in the 
> extensions so not officially released yet. To release them we need to 
> think about the interface first. The geographical strategies should 
> have an interface similar to the one used in spherical strategies. 

I don't know if that is so important. What is important is that we have 
one model or concept of the Earth, used in both distance strategies 
(andoyer, vincenty) and in the geographic projections.

That the interface is a bit different than Haversine, I don't know if 
that is a problem (maybe, but we can examine that). Or maybe they both 
can be implemented with a default constructor.

> Or the other way around, spherical strategies interface could be 
> extended somehow to match the geographical version. In all cases the 
> strategies should be backward compatible. Currently we have:
>
> CARTESIAN:
>
> - cartesian side has no ctor
>
> side_by_triangle<> s;
>
> - cartesian distance has no ctor
>
> pythagoras<> s;
>
> SPHERICAL:
>
> - spherical side has no ctor since it doesn't need to know the radius
>
> spherical_side_formula<> s;
>
> - spherical distance (haversine) requires a radius of type passed to 
> the strategy as template parameter, and by default radius = 1.0 is 
> used (unit sphere)
>
> haversine<RadiusType> s(radius);
> haversine<RadiusType> s; // radius = 1.0 by default
>
> so to calculate the distance on a sphere the user can pass a strategy 
> with radius or use the default one and multiply the result by the radius:
>
> distance(p1, p2, haversine<RadiusType>(radius));
> distance(p1, p2) * radius;
>
> GEOGRAPHICAL:
>
> - a spherical side strategy is just used by default which I guess 
> should be changed. I think the most precise would be a strategy 
> finding the side using the geographic courses found using the inverse 
> Vincenty's formula but this requires some testing

As far as I know, the spherical side strategy works for the geographic 
Earth too, please indicate if that is not the case.

See also
http://barendgehrels.blogspot.nl/2011/06/spherical-side-formula.html

But I might be wrong, it is good to research this.

> - for distance there are 2 strategies, both taking 
> detail::ellipsoid<RadiusType>, andoyer also has a ctor taking a 
> flattening of type RadiusType (like detail::ellipsoid<>). If 1 
> parameter is passed (flattening) then the radius is set to 1.0 but if 
> 2 parameters are passed then A and B radiuses are set. Furthermore by 
> default Earth's radiuses are used. So the situation looks like this:
>
> vincenty<RadiusType> s(detail::ellipsoid<RadiusType>()); // A = 
> 6378137.0, B = 6356752.314245
> vincenty<RadiusType> s(detail::ellipsoid<RadiusType>(f)); // A = 1.0, 
> B = [some fraction for oblate spheroid]
> vincenty<RadiusType> s(detail::ellipsoid<RadiusType>(a, b)); // A = a, 
> B = b
> andoyer<RadiusType> s; // A = 6378137.0, B = 6356752.314245
> andoyer<RadiusType> s(f); // A = 1.0, B = [some fraction for oblate 
> spheroid]
>
> It's confusing (this 1 and 2 paramers difference) and inconsistent 
> with spherical strategies.

Yes this should be harmonized with one "libary-universal" concept of the 
Earth model.

> Furthermore flattening should be stored using some Floating Point type 
> and currently it's just RadiusType. E.g. for integral RadiusType the 
> default ctor would initialize flattening with 0.
>
> PROPOSAL:
>
> 1. I propose to keep it simple for now, don't release 
> detail::ellipsoid and change the default parameters to make the 
> strategies close to the spherical ones:
>
> vincenty<RadiusType> s(); // A = 1.0, B = [fraction], F= 
> [fraction,default] - unit WGS84 spheroid
> vincenty<RadiusType> s(a); // A = a, B = [calculated from a using 
> default flattening], F= [fraction,default]
> vincenty<RadiusType> s(a, b); // A = a, B = b, F= [calculated from a 
> and b]
>
> This way the calls:
>
> distance(p1, p2, vincenty<RadiusType>(a));
> distance(p1, p2) * a;

I don't know - will anybody call Vincenty with the unit sphere?
I would prefer specifying the ellipsoid.

>
> would be equivalent (assuming that RadiusType is not integral type?).
>
> By default FlatteningType would be:
> - for floating point RadiusType <=> RadiusType
> - for integral RadiusType <=> promoted to FP of the same size (int-> 
> float, long long-> double)
> - for user defined types <=> UDT
> The FlatteningType could be optionally passed as a second template 
> parameter of a strategy however I'm not sure about it since it 
> wouldn't play nice with possible future improvement 3 (see below).
>
> ALTERNATIVES (FUTURE ADDITIONS)
>
> 2. To decrease the probability of a mistake the strategy could take 
> wrapped parameters (only?). This way it'd also be possible to create 
> an ellipsoid from A,B A,F or B,F pair:
>
> // below are functions returning wrapped values
> vincenty<RadiusType> s(major_radius(a), minor_radius(b));
> vincenty<RadiusType> s(major_radius(a), flattening(f));
> vincenty<RadiusType> s(minor_radius(b), inverse_flattening(if));
>
> 3. Desing concepts for models parameters. And then always pass a type 
> adapted to model's parameters concept to the strategy, instead of 
> RadiusType. Actually a single number could be adapted to e.g. sphere 
> parameters concept (radius) but why reserve it specifically for sphere?
>
> haversine< parameters::sphere<RadiusType> > 
> s(parameters::sphere<RadiusType>(radius));
> vincenty< parameters::ellipsoid<RadiusType> > 
> s(parameters::ellipsoid<RadiusType>(a, b));
> vincenty< parameters::ellipsoid_af<RadiusType> > 
> s(parameters::ellipsoid_af<RadiusType>(a, f));
> vincenty< parameters::ellipsoid_unit_flattening<RadiusType> > 
> s(parameters::ellipsoid_unit_flattening<RadiusType>(f));
> vincenty< parameters::ellipsoid_wgs84_1996<RadiusType> > 
> s(parameters::ellipsoid_wgs84_1996<RadiusType>());
> vincenty< parameters::ellipsoid_wgs84_2004<RadiusType> > 
> s(parameters::ellipsoid_wgs84_2004<RadiusType>());
> vincenty< parameters::ellipsoid_latest<RadiusType> > 
> s(parameters::ellipsoid_latest<RadiusType>());
>
> FYI, all names of namespaces, structs and functions are examples.
>
> For sphere we'd already have a small problem with this. However this 
> could be implemented in a backward compatible way. But then it could 
> be done this way in the future, for both coordinate systems.
>
> 4. 2 and 3 together:
>
> vincenty< parameters::ellipsoid<RadiusType> > 
> s(parameters::ellipsoid<RadiusType>(major_radius(a), minor_radius(b)));
>
> Do you have any thoughts or requests regarding the above?
>

Yes, I opt for a more simple model, one constructor and not many 
overloads, with the model of the Earth. As discussed earlier. So 
something like this, as you suggest:

vincenty<RadiusType> s(detail::ellipsoid<RadiusType>()); // A = 
6378137.0, B = 6356752.314245
andoyer<RadiusType> s(detail::ellipsoid<RadiusType>()); // A = 
6378137.0, B = 6356752.314245

(but the ellipsoid should not reside in detail)

This harmonization is the main reason why it is not moved earlier from 
extensions to the main library. As soon as this is fixed we can probably 
release them together. The main work will probably be in the projection 
part.

Regards, Barend