Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
PeskyGnat opened this issue Dec 16, 2015 · 3 comments
Labels
Milestone

Comments

@PeskyGnat
Copy link

@PeskyGnat PeskyGnat commented Dec 16, 2015

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?

@mguignard
Copy link

@mguignard mguignard commented Dec 17, 2015

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
Copy link
Contributor

@hobu hobu commented Dec 17, 2015

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?

@mguignard
Copy link

@mguignard mguignard commented Dec 24, 2015

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.

kbevers added a commit to kbevers/PROJ that referenced this issue Oct 19, 2016
@kbevers kbevers closed this in 2cac138 Oct 19, 2016
@kbevers kbevers added this to the 4.10.0 milestone Oct 19, 2016
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Linked pull requests

Successfully merging a pull request may close this issue.

None yet
4 participants
You can’t perform that action at this time.