Opened 10 years ago

Closed 10 years ago

#3623 closed defect (fixed)

UTM Zone Computation Stability Issue

Reported by: warmerdam Owned by: ilucena
Priority: normal Milestone:
Component: OGR_SRS Version: 1.7.2
Severity: normal Keywords:
Cc: warmerdam

Description

Dear Mr. Warmerdam,

for some technical reasons, I took a look into the GDAL library, where I found the following issue I ran into during testing:

File Id: $Id: ogrspatialreference.cpp 18544 2010-01-14 06:47:29Z warmerdam $

Starting in line 4705, you have the following code to compute a UTM zone:

    double      dfZone = (dfCentralMeridian+183) / 6.0 + 0.000000001;

    if( ABS(dfZone - (int) dfZone) > 0.00001
        || dfCentralMeridian < -177.00001
        || dfCentralMeridian > 177.000001 )
        return 0;
    else
        return (int) dfZone;

This is numerically unstable, since the integer jump is taking place at the zones' centers, which you compensate by adding 0.000000001. When placing the integer jump at the zones' borders, my problems went away:

    double      dfZone = (dfCentralMeridian + 186.0) / 6.0; /* <--
changed */

    if( ABS(dfZone - (int) dfZone - 0.5) > 0.00001 /* <-- changed */
        || dfCentralMeridian < -177.00001
        || dfCentralMeridian > 177.000001 )
        return 0;
    else
        return (int) dfZone;

With best regards from Munich, Yours Henning Lorch

Attachments (1)

getutmzone.cpp (1.1 KB) - added by ilucena 10 years ago.
unit test fro GetUTMZone change

Download all attachments as: .zip

Change History (6)

comment:1 Changed 10 years ago by warmerdam

Henning writes:

as a temporary solution and workaround around some other business solutions, I was porting a piece of your code into a different programming language, though still using double precision. There I got a roundoff error which was negative and above 0.000000001 (what I didn't expect...).

That doesn't mean that the problem has to show up when using a C-compiled library. Anyway I thought I'd rather tell you.

comment:2 Changed 10 years ago by warmerdam

Cc: warmerdam added
Owner: changed from warmerdam to ilucena

Ivan,

I'd appreciate your looking into this when you get some time. I suspect the code in GetUTMZone() is only used with actual zone central meridians so it might be fine as is, but if it can be easily made more robust that might be a nice idea.

Changed 10 years ago by ilucena

Attachment: getutmzone.cpp added

unit test fro GetUTMZone change

comment:3 Changed 10 years ago by ilucena

Frank Warmerdam wrote:

Lucena, Ivan wrote:

Frank,

I wrote this "unit test" code here to check the difference of the two methods:

...

If that doesn't raise the red flag I guess we can do the change.

What do you think?

Ivan,

Have you confirmed that the various longitude values are being assigned to the proper zones now, and were not before?

-117 is the central meridian of zone 11. So it is appropriate that -118 and -119 are also assigned to zone 11, while -113 and -112 clearly belong to zone 12. -120 and -114 are "on the line" and could go either way.

So it seems the change is reasonable.

Please check that the end condition and handling across zero is clean.

This is a piece of the output from the test code (on attachment):

ctm,old,new,diff = -120,10,11,1 ctm,old,new,diff = -119,10,11,1 ctm,old,new,diff = -118,10,11,1 ctm,old,new,diff = -117,11,11,0 ctm,old,new,diff = -116,11,11,0 ctm,old,new,diff = -115,11,11,0 ctm,old,new,diff = -114,11,12,1 ctm,old,new,diff = -113,11,12,1 ctm,old,new,diff = -112,11,12,1 ctm,old,new,diff = -111,12,12,0 ctm,old,new,diff = -110,12,12,0 ctm,old,new,diff = -109,12,12,0 ... ctm,old,new,diff = -6,29,30,1 ctm,old,new,diff = -5,29,30,1 ctm,old,new,diff = -4,29,30,1 ctm,old,new,diff = -3,30,30,0 ctm,old,new,diff = -2,30,30,0 ctm,old,new,diff = -1,30,30,0 ctm,old,new,diff = 0,30,31,1 ctm,old,new,diff = 1,30,31,1 ctm,old,new,diff = 2,30,31,1 ctm,old,new,diff = 3,31,31,0 ctm,old,new,diff = 4,31,31,0 ctm,old,new,diff = 5,31,31,0

It looks like the problem is fixed so I will apply that change.

D'accord?

comment:4 Changed 10 years ago by ilucena

 ctm,old,new,diff = -120,10,11,1
 ctm,old,new,diff = -119,10,11,1
 ctm,old,new,diff = -118,10,11,1
 ctm,old,new,diff = -117,11,11,0
 ctm,old,new,diff = -116,11,11,0
 ctm,old,new,diff = -115,11,11,0
 ctm,old,new,diff = -114,11,12,1
 ctm,old,new,diff = -113,11,12,1
 ctm,old,new,diff = -112,11,12,1
 ctm,old,new,diff = -111,12,12,0
 ctm,old,new,diff = -110,12,12,0
 ctm,old,new,diff = -109,12,12,0
 ...
 ctm,old,new,diff = -6,29,30,1
 ctm,old,new,diff = -5,29,30,1
 ctm,old,new,diff = -4,29,30,1
 ctm,old,new,diff = -3,30,30,0
 ctm,old,new,diff = -2,30,30,0
 ctm,old,new,diff = -1,30,30,0
 ctm,old,new,diff = 0,30,31,1
 ctm,old,new,diff = 1,30,31,1
 ctm,old,new,diff = 2,30,31,1
 ctm,old,new,diff = 3,31,31,0
 ctm,old,new,diff = 4,31,31,0
 ctm,old,new,diff = 5,31,31,0

comment:5 Changed 10 years ago by ilucena

Resolution: fixed
Status: newclosed

Committed revision r19841

Note: See TracTickets for help on using tickets.