Add a dedicated proj=webmerc operation #925

rouault · 2018-04-07T10:30:56Z

Directly based on the auxiliary sphere type parameter of
http://desktop.arcgis.com/en/arcmap/10.3/guide-books/map-projections/mercator.htm

Implemented:

0: use semimajor axis or radius of the geographic coordinate system (WebMercator style)
1: use semiminor axis or radius
2: calculate and use authalic radius

Unimplemented:

3: use authalic radius and convert geodetic latitudes to authalic latitudes

rouault · 2018-04-07T10:32:27Z

Was also wondering if it would be worth to have a +proj=gmerc that would be an alias for +proj=merc +aux_sphere_type=0 ?

kbevers · 2018-04-07T11:02:15Z

This is nice and definitely a better solution to what we have now! This was on my list of things to do in the future so thanks for taking care of it.

Was also wondering if it would be worth to have a +proj=gmerc that would be an alias for +proj=merc +aux_sphere_type=0 ?

That would have been my approach. I think it is a cleaner solution although it is a little more work to add a new projection.
I would suggest not making it an alias but instead a "proper" projection in it's own right. As far as I can tell it would be fairly simple and mostly just requires an added PROJECTION setup function. I would add the auxiliary sphere stuff to that and keep it out of the Mercator setup. +lon_0=0, +lat_ts=0 and +k_0=1 could then also be fixed, simplifying the user interface. Maybe what you meant was along the lines of that solution - if so, please disregard the previous paragraph.

I think calling it webmerc or pseudomerc is more describing. Also comes with the benefit of being less tied to the name of a large company.

kbevers · 2018-04-07T11:03:29Z

src/PJ_merc.c

@@ -48,11 +49,36 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
    return lp;
 }

+/* Cf http://desktop.arcgis.com/en/arcmap/10.3/guide-books/map-projections/mercator.htm */
+#define AUX_SPHERE_TYPE_SEMI_MAJOR 0


Why not use a enum here? You could use a switch below instead of the else if's. I think that is cleaner.

rouault · 2018-04-07T11:15:18Z

I would add the auxiliary sphere stuff to that and keep it out of the Mercator setup

I really tried to implement ESRI stuff, which as far as I can see (from the doc ! not tested against Arc), is general mercator with all the the bells and whistles, and the ability to define how to compute the spherical radius

kbevers · 2018-04-07T11:19:49Z

Yeah, I figured that much. I am not opposed to it either, just trying to make sure we handle this as best as possible. Do you know of any other case, apart from the webmercator, where this is used?

kbevers · 2018-04-07T11:31:54Z

src/PJ_merc.c

+#define AUX_SPHERE_TYPE_SEMI_MAJOR 0
+#define AUX_SPHERE_TYPE_SEMI_MINOR 1
+#define AUX_SPHERE_TYPE_AUTHALIC_RADIUS 2
+/* Unimplemented #define AUX_SPHERE_TYPE_AUTHALIC_LATITUDE 3 */


healpix gives some hints as to how to implement this: https://github.com/OSGeo/proj.4/blob/master/src/PJ_healpix.c#L225-L246

rouault · 2018-04-07T11:45:26Z

Do you know of any other case, apart from the webmercator, where this is used?

No... Perhaps @melitakennedy could comment on the use cases for auxiliary_sphere_type=1/2/3 for the Mercator projection ?

kbevers · 2018-04-07T12:07:03Z

src/PJ_merc.c

@@ -76,3 +119,22 @@ PJ *PROJECTION(merc) {
    return P;
 }

+PJ *PROJECTION(webmerc) {
+
+    if( pj_param(P->ctx, P->params, "tk_0").i ) {


This behaviour is different to the rest of PROJ. Usually you are allowed to set a parameter even if it is not used. If it is not used by the projection it will just get ignored. proj will tell you if that is the case:

$ proj -v +proj=merc +lat_1=23 #Mercator # Cyl, Sph&Ell # lat_ts= # +proj=merc +ellps=WGS84 #--- following specified but NOT used # +lat_1=23

I think it would be better to keep this "tradition". Just make sure to override them in the setup with P->k0=1.0; P->lam0=0.0.

I kept the ones for k_0 and lat_0 since those are done in generic parts. I remove the one for lat_ts (which was wrongly spelled in my original commit)

I understand why you would want to do that but I still don't think it is the correct way to handle it. I think it will cause more confusion than is necessary. This would be the only place in the code base where using one of the default projection parameters would cause an error. There are plenty examples of overriding default parameters in the setup without raising errors.

I think it would be awesome if we had a system that caused errors when using parameters that can't be used with a given operation but that is a problem for another day (and major version release :-) ).

OK, checks reverted and parameters overriden

kbevers · 2018-04-07T12:09:28Z

test/gie/4D-API_cs2cs-style.gie

+accept      487147.594520173    4934316.46263998    0
+expect      -10370728.80        5552839.74          0
+-------------------------------------------------------------------------------
+


How about some tests for aux_sphere_type == 1 and 2?

Added, by the way while adding them I had originally issues with the z component being modified, which I think was wrong. I fixed that issue with commit ac78787 which I believe is a bugfix that should probably be applied to the 5.0 branch

Added, by the way while adding them I had originally issues with the z component being modified, which I think was wrong.

Do you mind writing that up in an issue so it can be referred to in the release notes later on?

I fixed that issue with commit 1a575b6 which I believe is a bugfix that should probably be applied to the 5.0 branch

I don't plan on doing more bug fix release for the 5.0 branch. 5.1 will be out in two months. I think that will be fine.

kbevers · 2018-04-07T13:08:13Z

Another thing, if the epsg file is updated to use +proj=webmerc we can re-instate the fix to #22. Would it cause problems downstream if we change the epsg-file in PROJ? Wil that have an effect on GDAL, QGIS, etc?

rouault · 2018-04-07T13:43:03Z

Would it cause problems downstream if we change the epsg-file in PROJ? Wil that have an effect on GDAL, QGIS, etc?

I'm not completely sure. GDAL should probably not parse itself the file, and QGIS probably not too. But the epsg file is generated by GDAL that for now will still generate the old hacky definition of web mercator. And we should expect that old hacky definition to be still passed for years !

rouault · 2018-04-07T13:59:26Z

By the way: this had apparently no impact on spherical Mercator but pj_calc_ellipsoid_params() behaviour is quite surprising: it doesn't recompute e, f and b is they are originally non zero ! I see that in ellps_spherification() there is a 'workaround' for that which I copied, but this feels messy

… original ellipsoid parameters (before any projection mess with them)

kbevers · 2018-04-07T14:56:11Z

I'm not completely sure. GDAL should probably not parse itself the file, and QGIS probably not too. But the epsg file is generated by GDAL that for now will still generate the old hacky definition of web mercator. And we should expect that old hacky definition to be still passed for years !

Alright. Do you plan on updating GDAL so it generates the epsg file using webmerc?
Well yes, but we can't keep supporting that hack forever! It would be nice if we could issue a warning when the hack is used. Not sure if we have a smart way of doing that at the moment though...

By the way: this had apparently no impact on spherical Mercator but pj_calc_ellipsoid_params() behaviour is quite surprising: it doesn't recompute e, f and b is they are originally non zero ! I see that in ellps_spherification() there is a 'workaround' for that which I copied, but this feels messy

It does indeed. I wonder what the reason for that is. Can you answer that @busstoptaktik ? It doesn't affect the tests (apart from one minor detail in the gie pj_cart_selftest function) if changed to not respect the original values.

rouault · 2018-04-07T15:54:16Z

Do you plan on updating GDAL so it generates the epsg file using webmerc?

Potentially, the exportToProj4() code could, once proj 5.1 is released and detected at runtime, export EPSG:3857 as webmerc.

we can't keep supporting that hack forever!

The hack is only in the pj_transform() legacy interface where WGS84 has always (and will remain until we remove it) the pivot datum, so I don't think that's a big issue.

kbevers · 2018-04-07T16:10:00Z

The hack is only in the pj_transform() legacy interface where WGS84 has always (and will remain until we remove it) the pivot datum, so I don't think that's a big issue.

That's true. I was thinking about the +nadgrids=@null part of the hack but we seem to have that covered in the cs2cs emulation functions so that is not a problem.

melitakennedy · 2018-04-09T16:10:08Z

Although I believe "web Mercator" was the instigator, Esri added several "Auxiliary Sphere" versions of existing projections. At the time, we were trying to match results from NGA's GEOTRANS and ERDAS. They both had several projections that were implemented as sphere-only, but Esri had ellipsoidal versions. One or the other (I don't remember which) would internally calculate the authalic radius and use that. Esri's default is to use the semimajor axis.

The developer, David Burrows, pointed out that for an equal-area sphere-only projection, calculating the authalic radius is not quite enough. The geodetic latitudes should be converted to authalic latitudes as well. As far as I know, no one using any of this except Mercator aux sphere.

kbevers · 2018-04-09T17:52:01Z

@melitakennedy Thanks for your insight.

It sounds to me the "Auxiliary Sphere" doesn't bring anything new to PROJ. We already have both ellipsoidal and spherical versions of the Mercator projection and with the introduction of webmerc we can also handle the Webmercator satisfactorily.

@rouault Based on the above I suggest removing the aux. sphere stuff. Do you agree with that sentiment?

rouault · 2018-04-09T18:10:09Z

Based on the above I suggest removing the aux. sphere stuff. Do you agree with that sentiment?

OK, I've just kept the new webmerc and dropped the aux_sphere_type option of merc

kbevers · 2018-04-09T20:38:09Z

OK, I've just kept the new webmerc and dropped the aux_sphere_type option of merc

Thanks!

kbevers reviewed Apr 7, 2018

View reviewed changes

rouault force-pushed the aux_sphere_type branch 2 times, most recently from 51df307 to 6596bad Compare April 7, 2018 11:42

rouault changed the title ~~Mercator projections: add a aux_sphere_type=0/1/2 parameter (cleaner modeling of WebMercator)~~ Mercator projection: add a aux_sphere_type=0/1/2 parameter (cleaner modeling of WebMercator), and add a dedicated proj=webmerc operation Apr 7, 2018

kbevers reviewed Apr 7, 2018

View reviewed changes

rouault force-pushed the aux_sphere_type branch from 6596bad to fd131c7 Compare April 7, 2018 13:36

rouault force-pushed the aux_sphere_type branch from fd131c7 to 25f65e1 Compare April 7, 2018 13:54

rouault force-pushed the aux_sphere_type branch from 25f65e1 to 646cfa0 Compare April 7, 2018 14:02

Pipeline: make sure geocentric/cartesian space transform is done with…

ac78787

… original ellipsoid parameters (before any projection mess with them)

rouault force-pushed the aux_sphere_type branch 2 times, most recently from c90bb74 to 7072ba4 Compare April 7, 2018 14:44

rouault force-pushed the aux_sphere_type branch from 7072ba4 to c7380c4 Compare April 7, 2018 16:20

rouault force-pushed the aux_sphere_type branch from c7380c4 to ee2d782 Compare April 9, 2018 18:08

rouault changed the title ~~Mercator projection: add a aux_sphere_type=0/1/2 parameter (cleaner modeling of WebMercator), and add a dedicated proj=webmerc operation~~ Add a dedicated proj=webmerc operation Apr 9, 2018

Add webmerc projection

3a2ddb6

rouault force-pushed the aux_sphere_type branch from ee2d782 to 3a2ddb6 Compare April 9, 2018 18:51

kbevers merged commit 875d516 into OSGeo:master Apr 9, 2018

kbevers added this to the 5.1.0 milestone Apr 10, 2018

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add a dedicated proj=webmerc operation #925

Add a dedicated proj=webmerc operation #925

rouault commented Apr 7, 2018

rouault commented Apr 7, 2018

kbevers commented Apr 7, 2018

kbevers Apr 7, 2018

rouault Apr 7, 2018

rouault commented Apr 7, 2018

kbevers commented Apr 7, 2018

kbevers Apr 7, 2018

rouault commented Apr 7, 2018

kbevers Apr 7, 2018

rouault Apr 7, 2018

kbevers Apr 7, 2018

rouault Apr 7, 2018

kbevers Apr 7, 2018

rouault Apr 7, 2018 •

edited

Loading

kbevers Apr 7, 2018

kbevers commented Apr 7, 2018

rouault commented Apr 7, 2018

rouault commented Apr 7, 2018

kbevers commented Apr 7, 2018 •

edited

Loading

rouault commented Apr 7, 2018

kbevers commented Apr 7, 2018

melitakennedy commented Apr 9, 2018

kbevers commented Apr 9, 2018

rouault commented Apr 9, 2018

kbevers commented Apr 9, 2018

Add a dedicated proj=webmerc operation #925

Add a dedicated proj=webmerc operation #925

Conversation

rouault commented Apr 7, 2018

rouault commented Apr 7, 2018

kbevers commented Apr 7, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rouault commented Apr 7, 2018

kbevers commented Apr 7, 2018

Choose a reason for hiding this comment

rouault commented Apr 7, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rouault Apr 7, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kbevers commented Apr 7, 2018

rouault commented Apr 7, 2018

rouault commented Apr 7, 2018

kbevers commented Apr 7, 2018 • edited Loading

rouault commented Apr 7, 2018

kbevers commented Apr 7, 2018

melitakennedy commented Apr 9, 2018

kbevers commented Apr 9, 2018

rouault commented Apr 9, 2018

kbevers commented Apr 9, 2018

rouault Apr 7, 2018 •

edited

Loading

kbevers commented Apr 7, 2018 •

edited

Loading