UTM Map Projection Produces Incorrect Results

QUESTION: I am trying to compare features on a GeoTiff file in a UTM map projection, using a WGS84 datum, in both IDL and ENVI, but I get different results. I believe the IDL results are wrong. Do you know anything about this?

ANSWER: You are correct. The IDL UTM map projection, set up with Map_Proj_Init, produces incorrect results.

Consider, for example, the location with an Easting value of 519279.36 and a Northing value of 7443030.98 in Zone 4 of a UTM projection, using the WGS84 datum. In IDL we would set this map projection up like this to calculate the corresponding latitude and longitude values.

   mapStruct = Map_Proj_Init('UTM', ELLIPSOID='WGS 84', ZONE=4)
   ll = Map_Proj_Inverse(519279.36, 7443030.98, Map_Structure=mapStruct)
   Print, 'Longitude: ', ll[0]
   Print, 'Latitude:  ', ll[1]

The values we obtain are these.

   Longitude:       -158.57325
   Latitude:         66.105565

If you put these same location values into UTM converters found at the Research Coordination Network, the The Engineering ToolBox, or the Chuck Taylor Toolbox, or into the proj4 Cartographic Projections Library, or even into ENVI itself, you get the following answer.

   Longitude:  -158.55594   
   Latitude:    67.103896

If you transform these values (which by overwhelming consensus we must consider correct) back into projected XY coordinates in IDL, you can get a sense of how far the IDL values are off.

   xy = Map_Proj_Forward(-158.55595, 67.103896, Map_Structure=mapStruct)
   Print, 'Easting Error:  ', Abs( 519279.36 - xy[0]), ' meters'
   Print, 'Northing Error: ', Abs(7443030.98 - xy[1]), ' meters'

We find these results.

   Easting Error:       0.14 meters
   Northing Error:    185.97 meters

You can see that this is a significant error.

Source of the UTM Projection Error

The source of the error appears to be that the IDL UTM map projection does not allow any datum except for a sphere, despite the fact that the documentation says otherwise and that other ellipsoids can be obviously be specified in the call to Map_Proj_Init.

You can see that this is so by looking at the A and E2 fields of the map structure. These fields hold the value for the semi-major axis and the squared value of the eccentricity (a measure of the flattening of the ellipsoid), respectively. For a WGS ellipsoid these fields should be set to 6378137.0 and 0.0066943800. But, as you can see, these are clearly set to the same values they would be set to if I had specified a spherical datum.

   IDL> Print, mapStruct.A, mapstruct.E2
        6370997.0    0.00000000
   IDL> m = Map_Proj_Init('UTM', ELLIPSOID='SPHERE', ZONE=4)
   IDL> Print, m.A, m.E2
        6370997.0    0.00000000

This problem arises from incorrect information for the UTM projection in the program map_proj_init_common.pro (available in the IDL distribution) on line 776. The line reads like this:

   '101 UTM  |clon|clat|zone|    |    |    |  |  |',$

When it should read like this:

   '101 UTM  |smaj|smin|zone|    |clon|clat|  |  |',$

Also, line 1204 on map_proj_init.pro should be changed from this:

   'UTM': projParams[0:1] = 0

To this:

   'UTM': 

But, unfortunately, even making these changes doesn't change the final result to the correct result. Which leads me to conclude that the problem must reside in the GCTP internal code (seems unlikely) or in the IDL interface to the CGTP code. I haven't received word (see the work-around section below) as to exactly what the problem is.

This is a problem because most satellite data uses the WGS84 datum to location latitude and longitude coordinates. If IDL cannot use this datum, then it becomes impossible to accurately navigate an image produced in this way. A datum transformation, which is not available in IDL itself, would be required to use this type of data.

This problem, it seems to me, simply points out how old the GCTP software is and how badly IDL's map projection software is in need of an update. Someone (hopefully ITTVIS) needs to implement the proj4 Cartographic Projections Library in IDL to bring the map projection software up to modern standards.

Technical Support Response and Work-Around

I have heard from the technical support folks this morning who acknowledge that there is a problem with the UTM projection when used with at WGS84 datum. They tell me this problem will be fixed in the next version of IDL (presumably IDL 8.2). They suggest using the Walbeck ellipsoid, which is nearly identical to the WGS84 ellisoid. Indeed, in my tests, the final results are off by less than a quarter of a meter.

   mapStruct = Map_Proj_Init('UTM', ELLIPSOID='WALBECK', ZONE=4)
   ll = Map_Proj_Inverse(519279.36, 7443030.98, Map_Structure=mapStruct)
   Print, 'Longitude: ', ll[0]
   Print, 'Latitude:  ', ll[1]

The values we obtain are these.

   Longitude:       -158.555938
   Latitude:         66.1038971

We can estimate the error with the known correct values.

   xy = Map_Proj_Forward(-158.555946, 67.103895, Map_Structure=mapStruct)
   Print, 'Easting Error:  ', Abs( 519279.36 - xy[0]), ' meters'
   Print, 'Northing Error: ', Abs(7443030.98 - xy[1]), ' meters'

We find these results, which appear to be acceptable.

   Easting Error:  0.020 meters
   Northing Error: 0.164 meters

This and other problems with map projection routines in IDL have been corrected in the dozen or so map projection routines available for the Coyote Graphics Library.

Version of IDL used to prepare this article: IDL 8.1.