Silent Substitution Affects Map Projection Results

QUESTION: I've read the IDL on-line documentation and I am sure I am setting up the Albers Equal Area map projection with Map_Proj_Init exactly as shown there. But, I am getting results with Map_Proj_Forward and Map_Proj_Inverse that are at odds with other software I have at my disposal. Can you explain why this is the case? Here is how I am create my map structure and print out the projected XY locations of Fort Collins, Colorado.

  IDL> mapStruct = Map_Proj_Init('Albers Equal Area', $
          ELLIPSOID='WGS 84', STANDARD_PAR1=21.0, STANDARD_PAR2=-19.0)
  IDL> Print, Map_Proj_Forward(-105.1, 40.6, Map_Structure=mapStruct)
      -10860690.0       4606730.4

In fact, I know these values to be wrong. The correct values are -10865540.0 and 453741.7, respectively. What am I doing wrong?

ANSWER: You are not doing anything wrong. Rather, something has been done to you without your knowledge.

While it is true that the documentation indicates that you can select the Albers Equal Area map projection in this way, if fails to mention that it will not select this map projection unless you also set the GCTP keyword. If you fail to set this keyword, then it is actually the Albers Equal Area Conic map projection (map projection index number 14, rather than 103) that is used. This map projection can only use a spherical datum, so the ELLIPSOID keyword value in your Map_Proj_Init call is silently ignored, and a sphere with a radius set to the length of the semi-major axis of the WGS 84 ellipsoid is silently used in its place.

In other words, this command is silently substituted for the command you used above.

  IDL> mapStruct = Map_Proj_Init('Albers Equal Area Conic', $
          SPHERE_RADIUS=6378137.0, STANDARD_PAR1=21.0, STANDARD_PAR2=-19.0)
  IDL> Print, Map_Proj_Forward(-105.1, 40.6, Map_Structure=mapStruct)
      -10860690.0       4606730.4

To get the correct values, set the GCTP keyword.

  IDL> mapStruct = Map_Proj_Init('Albers Equal Area', /GCTP, $
          ELLIPSOID='WGS 84', STANDARD_PAR1=21.0, STANDARD_PAR2=-19.0)
  IDL> Print, Map_Proj_Forward(-105.1, 40.6, Map_Structure=mapStruct)
      -10865540.0       4583741.7

The distance these values are off is not cataclysmic. It is roughly the same as locating IDL headquarters not in downtown Boulder, but at the bottom of the Boulder Reservoir, outside of town. You wouldn't be doing any damage to the enemy if you were using these numbers to target artillery, but they could be embarrassing, nevertheless. In fact, I think it's entirely possible that the first time you realize these numbers are wrong would be after that important paper you submitted using them has gone to print.

There are other map projection analogues between the old Map_Set projection routines and the more recent GCTP routines. Be sure you are choosing the right one by setting the GCTP keyword and avoid the pernicious problem of silent substitution.

These and other problems 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.