Geometry Boost List: Re: [geometry] Moving projections from extensions

Subject: Re: [geometry] Moving projections from extensions
From: Adam Wulkiewicz (adam.wulkiewicz_at_[hidden])
Date: 2016-06-23 19:00:30

Hi Barend,

Barend Gehrels wrote:
> Op 22-6-2016 om 18:39 schreef Adam Wulkiewicz:
>> Barend Gehrels wrote:
>>> Op 20-6-2016 om 17:59 schreef Adam Wulkiewicz:
>> <snip>
>>>
>>>>
>>>> So in general I propose to hide as much as we can and expose only 
>>>> crucial parts of the API.
>>>> For run-time projections use PIMPL idiom and wrap the pointer inside.
>>>> Require only to define CalculationType for projection and use 
>>>> double as default.
>>>> In transform strategy take projection as template parameter and use 
>>>> run-time projection as default.
>>>>
>>>> // compile-time
>>>> projections::robin_spheroid<>prj("+ellps=WGS84+units=m");
>>>> projections::robin_spheroid<>prj(srs::spheroid<>());
>>>> projections::forward(pt_ll,pt_xy, prj);
>>>>
>>>> // run-time
>>>> projections::projection<> prj("+proj=robin+ellps=WGS84+units=m");
>>>> projections::projection<> prj(epsg);
>>>> projections::forward(pt_ll,pt_xy, prj);
>>>>
>>>> // also runtime
>>>> projections::project_transformer<>projection("+proj=robin+ellps=WGS84+units=m");
>>>> transform(pt_ll,pt_xy,projection);
>>>
>>> This looks indeed simpler. And probably it is possible. This should 
>>> be changed in the generation.
>>>
>>
>> I was also thinking about an alternative API, based on projection tags:
>>
>>     // run-time
>> projections::projection<> prj("+proj=robin+ellps=WGS84+units=m");
>>
>> // compile-time
>> projections::projection<robin>prj("+ellps=WGS84+units=m");
>>
>> So robin (or robin_tag) would be a tag passed into projection<> as 
>> template parameter. For a default tag the projection would be runtime 
>> and for non-default tag the projection would be compile-time.
>
> This looks also good. The _tag suffix is mainly used internally, but 
> it is actually exposed to the user when registering new types, so I 
> think it could be used indeed. But without is also fine and more 
> conform with the initialization string.
>

Actually there is an error above, I should have used the proj4() 
wrapper, so e.g.:

projection<> prj(proj4("+proj=robin+ellps=WGS84+units=m"));

I'm not sure if this is a good thing that it's string-conformant. This 
is Proj4 specific string and the Proj4 names are very enigmatic. On the 
other they are compact and quite good for tags. See below.

>
>> And then it would be possible to pass only tag to transform strategy, 
>> so no need to pass the whole projection type.
>>
>> // run-time
>> projections::project_transformer<>projection("+proj=robin+ellps=WGS84+units=m");
>> transform(pt_ll,pt_xy,projection); // compile-time
>> projections::project_transformer<robin>projection("+ellps=WGS84+units=m");transform(pt_ll,pt_xy,projection);
>>
>> This could simplify the generic implementation. E.g. the constructors 
>> taking all of those possible input parameters directly in projection 
>> (mentioned below) could be implemented only once. Otherwise we'd be 
>> forced to implement all constructors in every static projection.
>>
>> On the other hand it could be harder to understand when projection is 
>> dynamic and when static.
>>
>> What do you think?
>
> It looks very good. I think the difference between static (specify in 
> source code) and dynamic (in initialization string) is clear enough.
>

I thought about it more and my conclusion is that not necessarily this 
has an advantage over the non-generic implementation in terms of lower 
ammount of code. It's because if we agreed that the common parameters 
could be ommited (because not all projections need the full set of 
parameters) then each projection could have its own parameters 
parsing/retrieval procedure. So we'd still end up writing 
projection-specific code.

Still I think that the API based on tags is more pleasant to use. It 
isn't required to pass the whole projection type into the transform 
strategy.

Which one do you prefer?

If someone else has a preference or any other idea it's a good time to 
express it. :)

>
>>
>>> Note that also possible should be:
>>>
>>> typedef boost::geometry::linestring<any point type> ls;
>>> ls line_ll, line_xy; // fill line_ll
>>> projections::forward(line_ls, line_xy, prj);
>>>
>>> So the line is automatically projected, where, if necessary, extra 
>>> points are added because straight lines will usually become curved 
>>> and vv is also possible. Same for most other geometry types 
>>> (multi_point excepted).
>>>
>>> This does not need to be implemented in the first release, but we 
>>> need to take this into account.
>>>
>>
>> Wouldn't it be the transform() function's purpose?
>
> Yes, that is true.

And it already returns bool so it's compatible with projections. I guess 
it's on purpose. :)

<snip>

>
>>
>>> Another way could be project_transformer automatically figuring out 
>>> the projection direction and explicit project_forward_transformer, 
>>> project_inverse_transformer strategies.
>>>>
>>>> I'm also confused about the names of projections. AFAIU the ones 
>>>> marked as xxx_ellipsoid map from ellipsoid of revolution or 
>>>> spheroid and the ones marked as xxx_spheroid map from a sphere. 
>>>> There is also e.g. cass_spheroid which AFAIU can work with 
>>>> ellipsoid of revolution. Are these simple inaccuracies or am I 
>>>> missing something?
>>> It is by the generation. I have to look this up, but if the proj4 
>>> source supports an ellipsoid for it (even if it does not make 
>>> sense), it is generated too. We can always adapt the generation.
>>>
>>>>
>>>> Currently the names are directly derived from Proj4 parameter names 
>>>> but in WKT they're different (PROJECTION["Transverse_Mercator"]) so 
>>>> users which doesn't know Proj4 could be confused.
>>>> What do you think about naming projections using full names, e.g.:
>>>>
>>>> projections::aea_ellipsoid -> projections::albers_equal_area
>>>> projections::cass_spheroid -> projections::cassini
>>>> projections::tmerc_ellipsoid -> projections::transverse_mercator
>>>
>>> I kept all the proj4 names and I think that is more convenient than 
>>> adding another translation table. For users, who usually know proj4 
>>> (it is quite well-known), it is also more convenient.
>>> Having said that, I agree that cassini looks better than just cass, 
>>> and if WKT has a complete map of all translations, we could consider 
>>> that too.
>>
>> And by "proj4 names" do you mean the names of functions used in the 
>> Proj4 sources?
>> Because the projection parameters has different names in Proj4 string 
>> format. I'd expect the users are more familiar with these names, not 
>> function names.
>>
>> E.g. for aea_ellipsoid (e.g. 
>> http://spatialreference.org/ref/epsg/2964/) I see that there are 
>> basically 3 formats/names:
>>
>> OGC WKT: PROJECTION["Albers_Conic_Equal_Area"]
>> Proj4: +proj=aea
>> ESRI WKT: PROJECTION["Albers"]
>
> OK

Let's compare some names:

Proj4 ESRIWKTOGCWKT
aea   AlbersAlbers_Conic_Equal_Area
aeqd Azimuthal_Equidistant   Azimuthal_Equidistant
cass Cassini                 Cassini_Soldner
cea   Cylindrical_Equal_Area Cylindrical_Equal_Area
lcc   Lambert_Conformal_Conic Lambert_Conformal_Conic_2SP
robin Robinson                Robinson

In general I find ESRI WKT names being some middle ground, they aren't 
as long as OGC WKT names but also not enigmatic as Proj4 ones. Still 
e.g. cylindrical_equal_area is quite long name.

project_transformer<cea>projection(proj4("+ellps=WGS84+units=m"));

project_transformer<cylindrical_equal_area>projection(proj4("+ellps=WGS84+units=m"));

We could come up with our own names, e.g. albers, cassini, robinson are 
quite good names. But not in all cases such good names could be created. 
I guess the author's name could always be full and the additional 
descriptors represented as single letters, e.g. albers_cea, aeqd (or 
aed, why "equal" is sometimes "e" and sometimes "eq"?), lambert_cc, 
lambert_aea, lambert_cca, cea, t_mercator (or mercator_t). But maybe 
it's a good solution to stick with Proj4 names because our projections 
are taken from this library.

Do you have a preference?

Regarding the names. There are xxx_ellipsoid and xxx_spheroid names even 
for the same projection, like in this case:
https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/extensions/gis/projections/proj/cea.hpp#L65
https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/extensions/gis/projections/proj/cea.hpp#L103
AFAIU the xxx_spheroid is a version optimized for sphere (e^2 == 0):
https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/extensions/gis/projections/proj/cea.hpp#L227
This is how this projection is implemented in Proj4, so the names are 
taken from Proj4 (e_forward, s_forward):
https://github.com/OSGeo/proj.4/blob/master/src/PJ_cea.c

So indeed if we want to allow to specify the exact version at compile 
time we need to have tags reflecting that, e.g.:

projection<cea_sphere>prj(/*...*/);
projection<cea_spheroid>prj(/*...*/);
projection<cea>prj(/*...*/);// same as spheroid

Where Proj4 s_ (Spheroidal) is BG _sphere and Proj4 e_ (Ellipsoidal) is 
BG _spheroid.

<TLDR>

or

projection<cea_spherical>prj(/*...*/);
projection<cea_geographic>prj(/*...*/);

or

// spherical_tag/geographic_tag or sphere/spheroid
projection<cea, spherical_tag>prj(/*...*/);
projection<cea, geographic_tag>prj(/*...*/); // default
projection<cea>prj(/*...*/); // same as geographic
// explicit version definition also for run-time projection?
projection<default_dynamic, spherical_tag>prj(/*...*/);
// otherwise it should be a different template
static_projection<cea, spherical_tag>prj(/*...*/);

Well, it'd also be possible to check eccentricity during the projection, 
each time but it'd add unnecessary overhead (plus there is additional 
initialization part for ellipsoidal version which woudl have to be run 
once, at the beginning). We could also pick the version based on the 
input Point types (but even in this case there is this initialization 
part which would have to be done once, so we'd still be forced to keep 
some is_initialized state). It'd also be possible to have a pointer to 
implementation but then it'd be similar to the run-time version so why 
bother with compile-time version. We could also simply allow to create 
only static ellipsoidal/spheroidal/geographic version and spherical 
optimization only for run-time projection, but this would be limiting.

On the other hand suppose the user creates the 
ellipsoidal/spheroidal/geographic static version but initializes it with 
+es=0 or bg::srs::sphere. It'd also be suboptimal. The general problem 
is that run-time parameters may be not "compatible" with compile-time 
projection/parameters. Maybe the solution is to NOT allow the user to 
initialize the compile-time projection with run-time parameters? Then 
the other parameters would have to be a part of a static projection 
type. Then we'd need separate template for static projections, taking 
additional static parameters, e.g.:

static_projection<cea, srs::spheroid<>/*, ...?*/>prj;
static_projection<cea, srs::sphere<>/*, ...?*/>prj;

static_projection<cea, srs::spheroid<> >prj(srs::spheroid<>(a, b)); // 
non WGS84 ellipsoid

Or separate type for each projection as it is now:

cea<srs::spheroid<> >prj;

But from here it's not far to:

cea_spheroid<> prj(srs::spheroid<>());

I don't see an elegant solution right now, at least more elegant than 
simple version-specific tags mentioned at the beginning and leaving the 
rest in the hands of the user.

</TLDR>

>
>>
>>>
>>>>
>>>> This would again be more clear and general. And wouldn't be 
>>>> confusing if we had tri-axial ellipsoid projections in the future.
>>>
>>> Thanks for your input.
>>
>> I'd like to do the necessary work now if you're still ok with that.
>
> Sure I am, it would be great to make new progress here. Do you also 
> want to update the conversion program?

If necessary, yes. Where is this program and what does it convert exactly?

Regards,
Adam