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

Subject: Re: [geometry] Support for geographic coordinate system
From: Krzysztof Czainski (1czajnik_at_[hidden])
Date: 2014-11-06 13:18:46

2014-11-06 17:06 GMT+01:00 Adam Wulkiewicz <adam.wulkiewicz_at_[hidden]>:

> Hi,
>
> I'm planning to begin the implementation of the support for geographical
> CS and would like to hear your opinion.

Hi Adam,

I just wanted to say that I'm glad you're working on this. I used the
geographical projections from extensions in a project I finished a year
ego. The project used a Boost release 1.54.0 (I think) with a snapshot of
Boost.Geometry from svn trunk up front in include dirs list. To unit-test
our use of the projections, we generated some test data in Matlab, and the
tests confirmed accurate results of projections.

> 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. 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
> - 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.
> 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;
>

would be equivalent (assuming that RadiusType is not integral type?).
>

If I recall correctly, we didn't use the distance algorithm in the end. But
at some prototyping stage we used distance(p1, p2, vincenty<RadiusType>()),
expecting the distance on Earth. You're proposed change would break this
expectation, but I don't think this is a problem.

>
> 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?
>

Thanks for a wonderful Geometry library and keep improving it!
Kris