UV_Box Ranges Are Incorrect in Map_Proj_Init

QUESTION: It appears to me that the UV ranges (also known as the projected XY range) returned in the uv_box field of map structure created by Map_Proj_Init are incorrect. For example, consider this.

   IDL> mapStruct = Map_Proj_Init('Equirectangular', Center_Longitude=90.0)
   IDL> Print, mapStruct.uv_box[[0,2]]
        -20015077.       18903129.

These values should be exactly the same, with opposite signs. In fact, since I am specifying a longitude range of -180 to 180 (by not using a LIMIT keyword), they should map to the same point. But they don't. Consider this.

   IDL> ll = Map_Proj_Inverse(mapStruct.uv_box[[0,2]], mapStruct.uv_box[[1,3]], $
             Map_Structure=mapStruct)
   IDL> Print, ll[0,0], ll[0,1]
         -90.000000      -100.00000

You see the values are off by 10 degrees.

ANSWER: Yes, this appears to be a bug in Map_Proj_Init when the center longitude is any value other than 0.0. (See another problem with the uv_box ranges in the next question below.) You can see the problem by running this example code.

    FOR j=0,10 DO BEGIN
       center_lon = 0 > j*36.0d < 360.0
       mapStruct = Map_Proj_Init(117, CENTER_LONGITUDE=center_lon)
       Print, center_lon, mapStruct.uv_box[[0,2]]
       
       ll = Map_Proj_Inverse(mapStruct.uv_box[[0,2]], $
              mapStruct.uv_box[[1,3]], MAP_STRUCTURE=mapStruct)
       Print, center_lon, ll[0,0], ll[0,1]
       Print, ""
    ENDFOR
    END

Here are some representative results.

      0.00000000      -20015077.       20015077.
      0.00000000       180.00000      -180.00000

       108.00000      -19792688.       19125518.
       108.00000      -70.000000      -80.000000

       180.00000      -20015077.       18903129.
       180.00000      0.00000000      -10.000000

       252.00000      -19125518.       19792688.
       252.00000       80.000000       70.000000

       360.00000      -20015077.       20015066.
       360.00000       180.00000       179.99990

You can see that the results are only correct when the center longitude is equal to zero.

This problem is corrected in the Coyote Graphics map projection routine, cgMap, which is a smart wrapper for the Map_Proj_Init command.

QUESTION: But, I think the uv_box is wrong even if I use a LIMIT keyword. If, as the IDL documentation states, the uv_box of the map structure "contains the map range in meters", then surely this is coincident with the limits of the map projection. But, it doesn't appear to be. Consider this example.

   limit =  [42.4000d, -92.9000d, 47.1000d, -86.8000d]
   mapStruct = Map_Proj_Init('Albers Equal Area', /GCTP, $
      Limit=limit, Standard_Par1=40, Standard_Par2=-39, Ellipsoid=8)
   xy = Map_Proj_Forward(limit[[1,3]], limit[[0,2]], MAP_STRUCTURE=mapStruct)
   Print, 'X Range: ', mapStruct.uv_box[[0,2]]
   Print, 'X Range: ', Reform(xy[0,*])
   Print, 'Y Range: ', mapStruct.uv_box[[1,3]]
   Print, 'Y Range: ', Reform(xy[1,*])

Here are the results. You can see clearly that the uv_box results are different from the projected limit results in the latitude values.

   X Range:       -7930061.5      -7404415.9
   X Range:       -7930061.5      -7404415.9

   Y Range:        5598899.5       6088932.0
   Y Range:        5604417.8       6083417.4

ANSWER: Yes, I think this incorrect, too. You can confirm this by inverse projecting the uv_box values back to latitude/longitude values.

   ll = Map_Proj_Inverse(mapStruct.uv_box[[0,2]], mapStruct.uv_box[[1,3]], MAP_STRUCTURE=mapStruct)
   Print, 'Lons: ', Reform(ll[0,*])
   Print, 'Lats: ', Reform(ll[1,*])

Here are the results.

   Lons:       -92.899293      -86.800663
   Lats:        42.348513       47.155763

Clearly, the latitude values are wrong by a significant amount.

In the cgMap object, I assume that the forward map projection code is working correctly (or else all is lost!) and I correct the map structure's uv_box field to correspond to the map limit. Here is an example.

   limit =  [42.4000d, -92.9000d, 47.1000d, -86.8000d]
   map = Obj_New('cgMap', 'Albers Equal Area', /GCTP, $
      Limit=limit, Standard_Par1=40, Standard_Par2=-39, Ellipsoid=8)
   s = map -> GetMapStruct()
   Print, 'X Range: ', s.uv_box[[0,2]]
   Print, 'Y Range: ', s.uv_box[[1,3]]

And, here are the results.

   X Range:       -7930061.4      -7404415.8
   Y Range:        5604418.1       6083417.7

And, when back projected:

   ll = map -> Inverse(s.uv_box[[0,2]], s.uv_box[[1,3]])
   Print, 'Lons: ', Reform(ll[0,*])
   Print, 'Lats: ', Reform(ll[1,*])

You get this:

   Lons:       -92.900000      -86.800000
   Lats:        42.400000       47.100000

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