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

Subject: [geometry] Support for geographic coordinate system
From: Adam Wulkiewicz (adam.wulkiewicz_at_[hidden])
Date: 2014-11-06 11:06:49

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

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?

Adam