Omerc forward projection problem

mraspaud commented 4 years ago

We have observed strange "folding" effects of the omerc projection when resampling satellite data, and are realizing that two different input lon/lat can have the same projection coordinates

Example of problem

$ proj +proj=omerc +lat_0=0.377041113875403 +lonc=33.8250934444081 +alpha=8.16321575614333 +gamma=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs -f %.5f   
-146.23 -78.19112222222222
875464.92038    -11348616.85746

$ invproj +proj=omerc +lat_0=0.377041113875403 +lonc=33.8250934444081 +alpha=8.16321575614333 +gamma=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs -f %.5f
875464.92038    -11348616.85746
-145.02309  -78.19112

Problem description

The omerc projection should not exhibit such behaviour, and it works fine in most other cases we tested.

Expected Output

I expect a round-trip lon/lat->x/y->lon/lat to give really close input and output lon/lats

Environment Information

PROJ version Rel. 6.2.1, November 1st, 2019
RHEL8, 64 bits

Installation method

Running in a miniconda environment

rouault commented 4 years ago

For others: this is not a problem with the inverse implementation, but actually both"-146.23 -78.19112222222222" and "-145.02309217 -78.19112222" project to the same coordinate.

Reading http://www.epsg.org/Portals/0/373-07-02.pdf , page 58 "3.6 Oblique Mercator" : "Hotine projected the ellipsoid conformally onto a sphere of constant total curvature, called the ‘aposphere’, before projection onto the plane and then rotation of the grid to north" (just seeing the same is mentionned in https://proj.org/operations/projections/omerc.html )

proj=omerc implements the Hotine method (I've quickly compared the formulas between the EPSG guidance note and PROJ and don't see any obvious blunder), so it is probably? expected that at this point, which is at ~ 180 deg of longitude from +lon_c , we don't have nice properties due to the aposphere trick ?? I've tried your case by using a sphere of +R=6378137 instead of the ellipsoid, and I don't get those multiple solutions. So seems a limitation of the Hotine projection specific to points far away from its center and on the ellipsoid.

(Take this with a grain of salt. I'm not a map projection designer)

mraspaud commented 4 years ago

@rouault thanks you very much for taking a look at this. I can confirm that using ellps=sphere solves the problem. Maybe it would be a good thing to add a note about this in the documentation, that using a non-spherical earth might be problematic ? Anyway, feel free to close this if you think it's not worth it :)

OSGeo / PROJ