Choose Wisely When Choosing Mercator

QUESTION: I have a GeoTIFF file in a Mercator map projection. But when I attempt to geolocate the corners of the image in latitude and longitude coordinates, the latitudes are off by nearly 20 degrees from what I know to be correct. I don't want to be picky, but this is rather substantial. It's the difference between hitting Moscow and hitting Baghdad, if you know what I mean. You can see how it might cause problems for me.

Is there something going on here that I'm not aware of?

ANSWER: Ah, my young knight. You must choose wisely when choosing the Mercator projection!

It seems there are two Mercator map projections in IDL. One gives the correct results, and the other, uh, does not. I strongly recommend against the latter.

Here is a partial listing of the listgeo output for the file in question.

   Geotiff_Information:
      Version: 1
      Key_Revision: 1.0
      Tagged_Information:
         ModelTiepointTag (2,3):
            0                0                0
            -167346.019      3677892.56       0
         ModelPixelScaleTag (1,3):
            463.312717       463.312717       0
         End_Of_Tags.

   Projection Method: CT_Mercator
      ProjNatOriginLatGeoKey: 20.000000 ( 20d 0' 0.00"N)
      ProjNatOriginLongGeoKey: -96.500000 ( 96d30' 0.00"W)
      ProjScaleAtNatOriginGeoKey: 1.000000
      ProjFalseEastingGeoKey: 0.000000 m
      ProjFalseNorthingGeoKey: 0.000000 m
   GCS: 4326/WGS 84
   Datum: 6326/World Geodetic System 1984
   Ellipsoid: 7030/WGS 84 (6378137.00,6356752.31)
   Prime Meridian: 8901/Greenwich (0.000000/  0d 0' 0.00"E)
   Projection Linear Units: 9001/metre (1.000000m)

   Corner Coordinates:
   Upper Left    ( -167346.019, 3677892.557)
   Lower Left    ( -167346.019, 2982923.482)
   Upper Right   ( 1222592.131, 3677892.557)
   Lower Right   ( 1222592.131, 2982923.482)

Using the information in the GeoTIFF file, we can calculate the corner values of this 3000x1500 GeoTIFF file ourselves.

   xscale  = 463.312717d
   yscale  = 463.312717d
   tp      = [-167346.019d, 3677892.56d ]
   s = [3000, 1500] ; Size of image
   
   ; Origin in upper-left in TIFF file. I want to move to bottom left
   ; for my own sanity.
   xOrigin = tp[0]
   yOrigin = tp[1] - (yscale * s[1])
   xEnd = xOrigin + (xscale * s[0])
   yEnd = tp[1]
      

If we print our results, they agree completely with the listgeo output from above.

   Corner Coordinates:
   Upper Left:  (-167346.019, 3677892.560)
   Lower Left:  (-167346.019, 2982923.485)
   Upper Right: (1222592.132, 3677892.560)
   Lower Right: (1222592.132, 2982923.485)

There are now two different Mercator map projections you can use in IDL. The old one, associated with Map_Set, and a newer one from the GCTP map projection library. The only apparent difference between them is that the former uses a CENTER_LATITUDE keyword to express the latitude at the center of the map projection, and the latter uses a TRUE_SCALE_LATITUDE keyword to indicate the latitude of true scale in the map projection.

Since the GeoTIFF file contains a tag named ProjCenterLatGeoKey, you might be tempted to choose the Mercator map projection that uses the CENTER_LATITUDE keyword to indicate the center latitude. You would have chosen poorly.

Here is the code you chose to use to convert opposite corners of the image into latitude and longitude coordinates.

    map = MAP_PROJ_INIT('Mercator', $
       DATUM=8, $
       CENTER_LATITUDE=20.0, $
       CENTER_LONGITUDE=-96.5)
    lonlat = MAP_PROJ_INVERSE([xorigin, xend],[yend, yorigin], $
        MAP_STRUCTURE=map)

If we print these results, we find this.

   Corners in Lon/Lat in "Normal" Mercator Projection
   Upper Left:  (-98.555, 51.338)
   Lower Right: (-82.367, 45.408)

While the longitude values (first value in each pair) are correct. The latitude values are wildly wrong. In fact, you should have chosen the GCTP Mercator projection, like this.

    map = MAP_PROJ_INIT('Mercator', /GCTP, $
       DATUM=8, $
       TRUE_SCALE_LATITUDE=20.0, $
       CENTER_LONGITUDE=-96.5)
    lonlat = MAP_PROJ_INVERSE([xorigin, xend],[yend, yorigin], $
        MAP_STRUCTURE=map)

Printing these results, we find these to be correct.

   Corners in Lon/Lat in GCTP Mercator Projection
   Upper Left:  (-98.099, 33.305)
   Lower Right: (-84.817, 27.554)

I am currently investigating why these results are so different from one another. I have a call into ITTVIS, too. I'll post the results of what I learn.

The word I get from ITTVIS technical support is that they don't understand the situation any better than I do, but they second my recommendation to use the GCTP map projections. Good advice, I would say. :-)

Version of IDL used to write this article: IDL 7.0.1.

Google
 
Web Coyote's Guide to IDL Programming