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-07 05:03:55

Hi Adam,

Adam Wulkiewicz wrote On 7-11-2014 0:00:
> Hi Barend,
>
> Barend Gehrels wrote:
>> Adam Wulkiewicz wrote On 6-11-2014 17:06:
>>> 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.
>>
>
> Yes, this is what I'd like to know, gather some feedback, etc.
> The interface may be different but it should be consistent, i.e. 
> familiar for someone already knowing the spherical strategies, the 
> default parameters should be intuitive in both cases, etc.

Ideally yes. But I don't think there are many users who will use both of 
them. It's a different world.

>
> <snip>
>
>>
>>>
>>> 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.
>
> I feel that Vincenty's formula should give different results in some 
> edge cases because the shape of geodesic on ellipsoid is different 
> than on a sphere but I must perform some tests to be sure.

I'm curious about 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.
>
> I'm not sure if one "universal" model would be sufficient. Or do you 
> think about the model for geographic CS? 

Yes.

> E.g. Earth approximated as a sphere could have some radius "not 
> compatible" with the ellipsoid model. E.g. mean radius may be 
> calculated in various ways. So I'm guessing that we shouldn't use 
> ellipsoidal model for spherical CS.

Agreed - besides that it is an overkill. Spherical really only has the 
radius... But yeah, if we have an ellips Concept which has an equal 
major and minor axis... That might be usable.

>
>>
>>> 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.
>>
>
> I suggested it to be consistent with the spherical CS.
> We could do other way around. Make the spherical CS consistent with 
> the geographical CS taking care about the backward compatibility.

Yes, I think that is the way to go.

>
> However I don't have an idea how to do it in a "good" way. If we just 
> allowed passing a Sphere model besides a raw radius it could have some 
> different default value (spherical Earth radius). But then depending 
> on the type the default would be different so this would be 
> unintuitive and error-prone:
>
> haversine<double> s; // radius = 1
> haversine<sphere<double>> s // radius = ... e.g. 6367446.98883 or one 
> of the other mean radius values
>
> Though this would be similar in the case of geographical CS, where 
> defaults are Earth's semi- major and minor axes.
>
> But if those models had a word 'Earth' in the name explicitly then 
> it'd be obvious what default values there may be:
>
> // it's obvious that by default those'd have some "real" radius values
> haversine<spherical_earth<double>>
> vincenty<ellipsoidal_earth<double>>
>
>>
>>>
>>> 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
>>
>
> Are you saying that if a user wanted to use some different A or B he 
> should define his own model instead of passing some parameters?

More specificly, the parameters are for the Earth model, not for the 
strategy. That is not really defining an own model (though that would be 
possible too), more configuring the existing model.

>
>
> What do you think about passing a model as a Strategy tempalte 
> parameter instead of a RadiusType?
> It'd be less error-prone, since the radius type would be defined in 
> one place (the model).
> It'd also simplify things in the case is we needed more than 1 type in 
> some future model.
>
> So:
>
> vincenty<EllipsoidModel> s(EllipsoidModel())
> andoyer<EllipsoidModel> s(EllipsoidModel())
> haversine<SphereModel> s(SphereModel())
> // or rather
> haversine<SphereModelOrRadius> s(SphereModel())
> // for backward compatibility
>
> , e.g.:
>
> typedef strategy::ellipsoid<double> ellipsoid;
> strategy::vincenty<ellipsoid> s(ellipsoid(a, b));
>
> vs
>
> typedef strategy::ellipsoid<double> ellipsoid;
> strategy::vincenty<double> s(ellipsoid(a, b));

Sure, the first version looks much better. And ellipsoid follows the 
Ellipsoid Concept (or how it will be called)

>
> We could then pass a Sphere model to the haversine<> and handle it 
> differently than raw radius (as mentioned above).

Sure.

>
>
> Btw, I think Vincenty's formula could also be used in side calculation 
> so it could also be a name for a side strategy. But let's leave the 
> names for later.
>
>> (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. 
>
> Should we also provide similar model for spherical CS with some 
> default mean Earth radius?
> Another possibility would be to somehow extract this radius from an 
> ellipsoid<>:
> - just take semi-major axis
> - use one of those mean values: 
> http://en.wikipedia.org/wiki/Earth_radius#Mean_radii
> - or one of those: 
> http://math.wikia.com/wiki/Ellipsoidal_quadratic_mean_radius
> So I'd rather leave it to the user.
>
> Or do you think otherwise?

I think a model following our Ellipsoid concept,  implementing a sphere, 
would be perfect.

>
> Btw, some time ago we talked about convenient strategies specification 
> for Variant geometries and this would be "in line" with it, e.g.:
>
> boost::variant<SphPoint, SphLinestring, SphPolygon> var1, var2;
> // ...
> distance(var1, var2, strategy::default & 
> strategy::sphere<RadiusType>(radius))
>
> And this object of type strategy::sphere<> could go directly into the 
> Strategy. But this is another story... 

I don't know - I skip this for now.

But you will have noticed my other mail about the projections too ;-)

Regards, Barend