tmerc "back folding" problems.

[warmerdam]

You can demonstrate the problem with the PROJ.4 commandline program cs2cs:

cs2cs +proj=latlong +datum=WGS84 +to +proj=utm +zone=32 +datum=WGS84

Give as input, this point in Alaska:

-149.66 64.237

And you get this:

522986.92       5029096.80 0.00

which is somewhere in Italy (near Milan).

I realize this result is an effect of using tmerc equations far from where they are valid, but it causes serious problems in MapServer? when trying to establish which images intersect a target map request. In this case we have images from around the world, some in relatively local projections like tmerc. We reproject the map request bounds into the target projection and check for overlap.

What I would like is if we could make the tmerc projection fail (ie return HUGE_VAL) for points so far from the central meridian that they are likely to start folding back on themselves.

[warmerdam]

Gerald writes (on the list):

etmerc has one limit level which is only a few degrees from del-lambda=90. The original code had several "levels" of error of with different return codes but because there was not anyway to do this multilevel return in libproj4 I just used the one at the furthest extent.

But the method of establishing a limit would appear to be done on the basis of Easting rather than a cutoff due to a lat-lon limit check. Perhaps do the calculation regarless of geographic coordinates and then reject based upon easting greater than a percentage of a or R or, of course, equation failure.

On inverse projection just reject if |easting| too large.

Allow an override option for those who want to play around. ;-)

[warmerdam]

Gerald I. Evenden wrote:

etmerc has one limit level which is only a few degrees from del-lambda=90. The original code had several "levels" of error of with different return codes but because there was not anyway to do this multilevel return in libproj4 I just used the one at the furthest extent.

Gerald,

I'm not sure what you mean above. I looked through the code for transverse mercator in PROJ.4, and libproj4 and I don't observe any obvious limit logic on lam or y.

But the method of establishing a limit would appear to be done on the basis of Easting rather than a cutoff due to a lat-lon limit check. Perhaps do the calculation regarless of geographic coordinates and then reject based upon easting greater than a percentage of a or R or, of course, equation failure.

I can understand this, but the problem is that the lon/lat to projected conversion is giving a plausible location near the central meridian in TM projected coordinates even though the lon/lat location was very very far from the central meridian. So I don't see how checking the easting after reprojection is going to be sufficient.

See:

 http://trac.osgeo.org/proj/ticket/5

On inverse projection just reject if |easting| too large. Allow an override option for those who want to play around. ;-)

I have taken the liberty of adding you as a cc: on the ticket. If you are interested in adding notes in the ticket, you would need to create an OSGeo login at:

 https://www.osgeo.org/cgi-bin/ldap_create_user.py

and then login using that.

I'm going to do a preliminary hack, but I'm hoping you will consider applying a more generalized/good solution upstream in libproj4.

[warmerdam]

As a preliminary step I have added the following logic in e_forward and s_forward in PJ_tmerc.c.

        /*
         * Fail if our longitude is more than 90 degrees from the 
         * central meridian since the results are essentially garbage. 
         * Is error -20 really an appropriate return value?
         * 
         *  http://trac.osgeo.org/proj/ticket/5
         */
        if( lp.lam < -HALFPI || lp.lam > HALFPI )
        {
            xy.x = HUGE_VAL;
            xy.y = HUGE_VAL;
            pj_errno = -14;
            return xy;
        }

[evenden]

Two major problems here:

1. grossly exceeding the defined logitude limits of UTM (nominally +-3.5 degrees from the central meridian and

2. grossly exceeding the limits of projection tmerc which starts to fall apart beyond 20 degrees from the central meridian.

I also feel that a great deal of the problem lies in the methodology of the MapServe? system. I would think that there would be failures of other projections that had to "stretch" their limits to meet the MapServer? coordinate needs.

Lastly, if one needs more extended TM coordinates use the either etmerc, ftmerc or ktmerc. Ktmerc will be available in the libproj4 library soon.

If one insists on the suggested patch to tmerc above, a more practical limit would be around 20 degrees rather than the 90 degrees suggest.

[heikok]

This problem seems to be general to tmerc with WGS-84 at > 7500000m north. I just run into the same problems in russia (proj 4.6.0 and proj 4.7.0):

$ proj +proj=tmerc +k=0.9996 +lon_0=9 +lat_0=0 +x_0=500000 +ellps=WGS84 63 63 2962014.55 8137552.85

$ proj -I +proj=tmerc +k=0.9996 +lon_0=9 +lat_0=0 +x_0=500000 +ellps=WGS84 2962014.55 8137552.85 29d57'1.587"E 66d22'34.344"N

Suddenly Russian points are placed on Greenland.

Doing the same with a spherical earth works well: $ proj +proj=tmerc +k=0.9996 +lon_0=9 +lat_0=0 +x_0=500000 +ellps=sphere +a=6371000 +e=0 63 63 2953651.67 8150383.02

$ proj -I +proj=tmerc +k=0.9996 +lon_0=9 +lat_0=0 +x_0=500000 +ellps=sphere +a=6371000 +e=0 2953651.67 8150383.02 63dE 63dN

[heikok]

And again, in readable Wiki syntax: Not working: WGS84

$ proj +proj=tmerc +k=0.9996 +lon_0=9 +lat_0=0 +x_0=500000 +ellps=WGS84
63 63
2962014.55 8137552.85

$ proj -I +proj=tmerc +k=0.9996 +lon_0=9 +lat_0=0 +x_0=500000 +ellps=WGS84
2962014.55 8137552.85
29d57'1.587"E 66d22'34.344"N 

Spherical example, working

$ proj +proj=tmerc +k=0.9996 +lon_0=9 +lat_0=0 +x_0=500000 +ellps=sphere +a=6371000 +e=0
63 63
2953651.67 8150383.02

$ proj -I +proj=tmerc +k=0.9996 +lon_0=9 +lat_0=0 +x_0=500000 +ellps=sphere +a=6371000 +e=0 
2953651.67 8150383.02
63dE 63dN

[heikok]

Checking even further, this seems to be a general problem of tmerc when latitude is > 10degree from the UTM-zone center, (x > 1000km).

At 30deg lon, error is approx.      7"
At 50deg lon, error is approx.   15'
At 60deg lon, error is approx. 1d47'

Latitude values don't affect the error. For UTM-zones, the tmerc is ok. But it cannot be used as general transverse Mercator.

If tested the (LGPL) code from the GeographicLib? from Karney:  http://geographiclib.sourceforge.net/tm.html (Transverse Mercator with an accuracy of a few nanometers) The provide two implementations (Krüger and Lee) and even the simple Krüger series expansion gives much better results than proj. Since the MIT and LGPL license don't match, I'm a bit hesitant to move that code to proj.

[heikok]

Sorry, I just saw that this issue has been discussed on the mailing list:  http://osgeo-org.1803224.n2.nabble.com/Transverse-Mercator-algorithm-with-good-accuracy-speed-trade-off-td2064100.html

I solved my problem by using the 'gstmerc' instead of the 'tmerc' projection, available since 4.7.

[mikrit]

Heikok wrote:

I solved my problem by using the 'gstmerc' instead of the 'tmerc' projection, available since 4.7.

That's fine when it works, but it could be risky as a general solution. Because gstmerc is the Gauss-Schreiber Transverse Mercator (in France known as Gauss-Laborde), which is a different projection. It is conformal, but the scale on the central meridian is only nearly constant, not exactly constant as in Gauss-Krüger Transverse Mercator. I would think that the differences between the two projections would accumulate over large areas. For a small area, though, it seems Gauss-Krüger can be approximated quite well by Gauss-Schreiber, at least near the latitude of origin:

"The differences in the various methods for the SPCS [U.S. State Plane Coordinate Systems] will NOT vary by more than a tenth of a foot or so - within the design area of the zone. --- The NAD27 Transverse Mercator was the Gauss-Schreiber math model, the NAD83 Transverse Mercator uses the Gauss-Krueger math model, etc." -- Clifford J. Mugnier,  http://newsgroups.derkeiler.com/Archive/Sci/sci.engr.surveying/2006-06/msg00093.html

Regards, Mikrit