tmerc "back folding" problems.

[proj4-bot]

Reported by warmerdam on 17 Jun 2008 15:46 UTC 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.

Migrated-From: https://trac.osgeo.org/proj/ticket/5

[proj4-bot]

Comment by warmerdam on 18 Jun 2008 00:56 UTC 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. ;-)

[proj4-bot]

Comment by warmerdam on 18 Jun 2008 01:16 UTC 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:

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.

[proj4-bot]

Comment by warmerdam on 18 Jun 2008 06:22 UTC

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?
         * 
         *  #5
         */
        if( lp.lam < -HALFPI || lp.lam > HALFPI )
        {
            xy.x = HUGE_VAL;
            xy.y = HUGE_VAL;
            pj_errno = -14;
            return xy;
        }

[proj4-bot]

Comment by evenden on 18 Jun 2008 17:44 UTC 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.

[proj4-bot]

Comment by heikok on 25 Aug 2010 09:03 UTC 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

[proj4-bot]

Comment by heikok on 25 Aug 2010 09:07 UTC 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

[proj4-bot]

Comment by heikok on 26 Aug 2010 07:55 UTC 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 (Krger and Lee) and even the simple Krger 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.

[proj4-bot]

Comment by heikok on 26 Aug 2010 09:12 UTC 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.

[proj4-bot]

Comment by mikrit on 12 Oct 2010 12:38 UTC 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-Krger 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-Krger can be approximated quite well by Gauss-Schreiber, at least near the latitude of origin:

"The differences in the various methods for the SPCS 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

[proj4-bot]

Comment by qulogic on 2 Mar 2015 01:34 UTC Replying to heikok:

The provide two implementations (Krger and Lee) and even the simple Krger 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.

Since !GeographicLib appears to be under the MIT/X11 License now, perhaps it is time to re-think this idea.

[proj4-bot]

Comment by mikrit on 2 Mar 2015 07:45 UTC Since Proj 4.8.0, there is an alternative Transverse Mercator implementation with an extended zone of accuracy, +proj=etmerc. See #97

I believe etmerc can be slower than tmerc, which is why etmerc was added as a separate projection, instead of as a refactoring of the old tmerc.

[QuLogic]

@cffk It is mentioned here that GeographicLib's implementation of TM is better than proj.4's. Since proj.4 is already using the geodesic portion of GeographicLib, perhaps you might be interested in also porting the TM implementation as well?

[cffk]

I extended etmerc to 6th order a year or so ago -- so now it gives full accuracy for a wide zone; the algorithms already match those of GeographicLib!

[QuLogic]

Ah, so I guess it's just a matter of documentation, then.

[cffk]

We still need to verify the fix to the "back folding" problem, which is somewhat of an independent issue compared with the straight accuracy of the projection. I'm on holiday now -- so I can't check this out at the moment; I'll be back the middle of next week.

[kbevers]

@cffk Did you ever get a chance to look at this?

[cffk]

No. But I'll take a look. There's maybe some discussion as to what the behavior will be. So at a minimum I'll document what the current behavior is. (I'll deal with this over the weekend.)

[kbevers]

Thanks, Charles. Since the problems presented here seem to be caused by using the projection way outside the defined area , i think documenting the behaviour would be enough.

[cffk]

This bug is fixed! The example in the bug report now returns

echo -149.66 64.237 |
cs2cs +proj=latlong +datum=WGS84 +to +proj=utm +zone=32 +datum=WGS84
==>
-519738.62 12698715.46 0.00

This is the correct result obtained by continuing the central meridan for zone 32 (lon = 9) to its anti-meridian (lon = -171). To see this repeat the calculation for the "opposite" zone:

echo -149.66 64.237 |
cs2cs +proj=latlong +datum=WGS84 +to +proj=utm +zone=2 +datum=WGS84
==>
1519738.62 7297214.42 0.00

The sum of these two results is

1000000.00 19995929.88

which is

2*(UTM-false-easting) 2*(UTM-central-scale)*(quarter-meridian)

as required.

The inverse projection also "works"

echo -519738.62 12698715.46 0.00 |
cs2cs +proj=utm +zone=32 +datum=WGS84 +to +proj=latlong +datum=WGS84 \
  -f "%.3f"
==>
-149.660 64.237 0.000

I haven't closely examined the code to see what its limits are. However I believe the evidence above is good enough to close the issue.

[kbevers]

That's handy! Do you have any idea with which changeset it was fixed? I've recently fixed a few other issues regarding tmerc, fixing this bug could very well be a side effect from that.

[cffk]

Probably it was a result of switching utm to use etmerc instead of tmerc. I would expect tmerc (which is an expansion in terms of longitude difference) to get this test case wrong. etmerc is an expansion in terms of flattening and so had much wider applicability.

[kbevers]

That's probably it then. I am just trying to figure out if this is something that needs to be mentioned in the changelog at the next release. But if it is fixed by the etmerc alias then it's a 4.9.3 thing...

[cffk]

Quite so. By the way, when tmerc is being documented, it should be clearly labeled as a "legacy" implementation. etmerc is recommended for all new applications. Perhaps now would be a good time to warn users that tmerc will (at some point in the future) be replaced be etmerc with the existing tmerc renamed to oldtmerc, or whatever.