Bug in Oblique Mercator projection in Proj4

[iskander]

Dear Sir,

Within VITO-TAP we are using Proj4 in a earth observation data processing chain we are currently implementing.

When we define the following Oblique Mercator projection:

+proj=omerc +ellps=WGS84 +lon_1=62.581150 +lat_1=74.856102 +lon_2=53.942810 +lat_2=74.905884 +no_rot +units=km

we get the following results with cs2cs:

    >56.958381652832 72.8798
    -10018.75	-227.68 0.00          <<< incorrect !
    >56.9584 72.8798
    9985.16	-227.68 0.00

When we change in PJ_omerc.c the following part (which, by the way, is conform with what we found in "Map Projections - A working manual" by John P. Snyder):

    con = cos(P->bl * lp.lam);
    if (fabs(con) >= TOL) {
        us = P->al * atan((s * P->cosgam + vl * P->singam) / con) / P->bl;
        if (con < 0.)
            us += PI * P->al / P->bl;
    } else
        us = P->al * P->bl * lp.lam;

to:

    con = cos(P->bl * lp.lam);
    if (fabs(con) >= TOL)
        us = P->al * atan((s * P->cosgam + vl * P->singam) / con) / P->bl;
    else
        us = P->al * lp.lam;
    if (con < 0.)
        us += PI * P->al / P->bl;

we get:

    >56.958381652832 72.8798
    9985.16	-227.68 0.00          <<< correct !
    >56.9584 72.8798
    9985.16	-227.68 0.00

The modification followed from the observation that con --> 0 not only implies atan(...) --> pi/2, but also P->bl * lp.lam --> pi/2. Hence, we can do the following:

    us = al * atan(...) / bl = al * (pi/2) / bl = al * (bl * lam) / bl = al * lam. 

Furthermore, it should be clear that the correction for con < 0 also still has to be applied.

Regards,

W. Heyns
VITO-TAP
B-2400 Mol, Belgium.

[warmerdam]

In (r1817) the PJ_omerc.c code was largely replaced by code from Gerald's libproj4 and the quoted code above no longer appears to be the same. However, I observe that the (erroneous) results are the same so the problem persists.

I'm going to email Gerald to see what he thinks but I won't likely get this resolved before 4.8.0 is released.