Unexpected Oblique Mercator (omerc) inverse proj result when alpha is between 90 and 270

[KodingForKittehs]

PROJ.4 version : v4.9.2, 08 September 2015

I was working on a module to allow a user to specify a local grid using the 'omerc', basically following what was suggested here:

http://gis.stackexchange.com/questions/83861/using-customized-coordinate-system-for-archaeological-site-data

This appears to work fine when the alpha angle is 270 to 360 and 0 to 90. If I specify an alpha greater than 90 but less than 270, I get strange results, see my test case below:

invproj +proj=omerc +lat_0=45 +lonc=-80 +x_0=0 +y_0=0 +alpha=45 +gamma=0 +k_0=1
0 0
80dW    45dN  (this is expected)

invproj +proj=omerc +lat_0=45 +lonc=-80 +x_0=0 +y_0=0 +alpha=89 +gamma=0 +k_0=1
0 0
80dW    45dN  (this is expected)

invproj +proj=omerc +lat_0=45 +lonc=-80 +x_0=0 +y_0=0 +alpha=91 +gamma=0 +k_0=1
0 0
77d10'18.699"W  45dN  (this is not expected)

I would expect that proj should be able to handle this because reversing the angle from 91 to 271 does return the correct answer.

Is this a bug or is it expected behavior? is there a way to adjust the parameters to get the behavior I'm looking for?

[cazimuth]

I compared the implementation with the reference at http://www.ihsenergy.com/epsg/guid7.pdf, on page 34, there is

In general
uC = (A / B) atan[(D2
 – 1)0.5 / cos (αC) ]. SIGN(ϕC

In PJ_omerc.c, the code is using atan2 which is giving different results when alpha is greater than 90 degrees, I tried replacing the atan2 with the atan from the reference and I now get correct results in my test case (I've not tested anything else)

        //P->u_0 = fabs(P->ArB * atan2(sqrt(D * D - 1.), cos(alpha_c)));
        P->u_0 = fabs(P->ArB * atan(sqrt(D * D - 1.)/ cos(alpha_c)));

I'm not sure why atan2 is being used here, I think it could be an issue, but changing it might break something else relying on this behavior.

[hobu]

Is this a bug or is it expected behavior? is there a way to adjust the parameters to get the behavior I'm looking for?

Probably a bug. Can you bring this issue up to the mailing list and see if anyone can provide some guidance for us?

[cazimuth]

Sorry that I haven't looked at this for a bit, I did look through the archives and found a message that is essentially the same issue that I ran into:

http://lists.maptools.org/pipermail/proj/2012-June/006331.html

I have a workaround for the moment, I subtract 180 degrees from the alpha and inverse the sign on easting/northing.

[rich9403]

I believe this issue was not correctly resolved. My recent attempts to perform omerc transforms with pyproj when +alpha is between 90 and 270 degrees are not working correctly (see issue https://github.com/pyproj4/pyproj/issues/1247). After investigating the issue with OSGeo via the shell, I believe there to be an existing problem within OSGeo/PROJ that pyproj is inheriting.

Varying the +alpha value an equal offset from the 90 or 270 degrees provides the same coordinate transform.

cs2cs +init=epsg:4326 +to +proj=omerc +lat_0=0 +lonc=0 +alpha=260 +x_0=0 +y_0=0 +gamma=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0

[rich9403]

I have attempted to search code to figure out how the "omerc" transform is calculated, but the all searches in the code of this repository do not show up with a github search. Any idea what is going on?

I also attempted to find the file PJ_omerc.c where this bug was located in an earlier comment on this subject, but I could not located this file in the repository. Was this file possibly removed with a version change and the same problem re-inherited from a new source?

[rouault]

I have attempted to search code to figure out how the "omerc" transform is calculated,

==> https://github.com/OSGeo/PROJ/blob/master/src/projections/omerc.cpp

[busstoptaktik]

I have attempted to search code to figure out how the "omerc" transform is calculated, but the all searches in the code of this repository do not show up with a github search. Any idea what is going on?

Yes, you're right - the search function does not work at all on the PROJ repo, which is quite annoying. Unfortunately, I do not have the slightest idea of why, but perhaps the widespread wrapping of everything into the anonymous namespace has to do with this? @rouault could this be the reason or do you have any other suggestions? (perhaps I should open a separate issue on this)

[rouault]

the widespread wrapping of everything into the anonymous namespace has to do with this?

?!? I don't see the connection.
Even looking for the main function used by utilities doesn't work (https://github.com/OSGeo/PROJ/search?q=main) and it is definitely not wrapped in a anonymous namespace

perhaps I should open a separate issue on this

to github itself. I don't see what we can do on our side.

[busstoptaktik]

to github itself. I don't see what we can do on our side.

But this appears to be specific to the PROJ repo, so not a generic Github thing, but something relating to "whatever is unique to PROJ", so prior to opening an issue with github herself, it might be nice to have an list of PROJ specifics that might trigger the problem

[rouault]

it might be nice to have an list of PROJ specifics that might trigger the problem

basically nothing works. I've tried "main", "proj_create", "proj_trans", etc.

[busstoptaktik]

I've tried "main",

So have I, and I'm puzzled: Apparently this only happens in our repo, thus my suspicion it had to do with something PROJ specific

[busstoptaktik]

I have opened a bug report on https://support.github.com/contact/bug-report, cf. https://support.github.com/ticket/personal/0/2069797

[mwtoews]

The GitHub "search" feature seams to work now, so it must have been a glitch.

A better "search" in a local git clone is with git-grep, e.g. git grep omerc.

[busstoptaktik]

so it must have been a glitch

In that case it must have been a many-months-long glitch :-) I first noticed it some time last fall. But nice to hear it works now!

[busstoptaktik]

I first noticed it some time last fall

Sorry, @mwtoews - should have been "last spring" at your latitudes. Probably I should have said "some time in October or November of 2022 current epoch"