Some notes on IDL maps

This page gives some notes about using the IDL command map_set and related commands. map_set sets up the map coordinate system which is then used by the plots command and also cursor,/data. The first section gives the available map projections and also some map_set keywords not documented in the manual. The next section gives details on how IDL saves the map coordinate system information.

There are no links on this page. Color is used to highlight information.

Map projections and some undocumented keywords

Available map projections:

Code            Projection Name
----    --      --------------
  1     St      Stereographic
  2     Or      Orthographic
  3     LC      LambertConic
  4     La      LambertAzimuthal
  5     Gn      Gnomic
  6     Az      AzimuthalEquidistant
  7     Sa      Satellite
  8     Cy      Cylindrical
  9     Me      Mercator
  10    Mo      Mollweide
  11    Si      Sinusoidal
  12    Ai      Aitoff
  13    Ha      HammerAitoff
  14    Al      AlbersEqualAreaConic
  15    Tm      TransverseMercator
  16    Mi      MillerCylindrical
  17    Ro      Robinson
  18    Le      LambertConicEllipsoid*
  19    Go      GoodesHomolosine

  * Not a real case.  Use proj=3 with Ellipsoid.

There are 3 ways to set map projection in map_set:

• Use a projection keyword as listed in the manual.
  Ex: /orthographic, /mercator, ...
• Use the undocumented keyword PROJECTION=code keyword.
• Use the undocumented keyword NAME=name keyword.

The map_set call

A typical map_set call looks like:

   map_set, 40, -90, 30, /ortho

The 3 values in the call are: lat0, lon0, ang which are the central latitude, central longitude, and angle of the pole from vertical. The map projection is set by the keyword /ortho which sets Orthographic. There are many other keywords allowed, most are documented in the manual or online help. A major goal of this routine is to set up the map coordinate system. The coordinate system information is stored in a number of system variables, !map is one of the main ones, but the complete information is spread over a number of system variables as shown below.

Retrieving the projection information

Why would you want to retrieve the map projection information? One example is to allow the map projection to be set up automatically at a later time. One way is to simply save the map_set command for later use. This is not so easy in practice since many of the values may be given as variables. But the information needed is mostly stored in system variables as shown below.

The information needed to define the coordinate system is:

• The central latitude.
• The central longitude.
• The map rotation angle.
• The map projection.
• The map screen position.
• If the /ISOTOPIC keyword was used.
• If the /NOBORDER keyword was used.
• The map limits (4 or 8 elements).

In addition there are several special case parameters:

• STD_P: Standard Parallel or Parallels for Albers (proj=14) and Conic (proj=3).
• SAT_P: Satellite Parameters for Satellite (proj=7).
• ELLIPS: Ellipsoid Parameters for Transverse Mercator (proj=15) or Lambert Conic Ellipsoid (proj=18, but use proj=3 with ELLIPS).

To see if a map_set call has set up the latest coordinate system check that !x.type = 3. If not then no map projection is available or it has been lost by another plot command. If the map scaling is available the needed values may be extracted as described below.

First, the 3 positional parameters are as follows:

  lat0 = !map.p0lat
  lng0 = !map.p0lon
  ang  = !map.rotation

The map projection that was used is found from:

   code = !map.projection

where the code is given by the table in the first section (and may be used with the PROJECTION=code keyword.

Many other keywords allowed in map_set set other values. Not all are needed to set up the projection coordinate system at a later time (for example plot colors, linestyles, thicknesses, ...).

The map screen position is easy to find:

   pos = [!x.window(0),!y.window(0),!x.window(1),!y.window(1)]

It is possible to determine if the /ISOTROPIC keyword was used. The X and Y scale factors are set by map_set. If /ISOTROPIC was used these scale factors will have the same ratio as the X and Y size of the window. However a simple test for equality of this ratio will sometimes run into minor roundoff error, the following works ok:

   shape = float(!d.y_size)/!d.x_size
   sc_rat = !x.s[1]/!y.s[1]
   if abs((sc_rat-shape))/sc_rat*10000d0 lt 1 then iflag=1 else iflag=0

It seemed like there was no way to tell if the /NOBORDER keyword was used, but a routine, posbox, was written to make a very good attempt at determining if it was. Seems to work all the time. It also returns the color of the border so it may be replotted without harm. The call is:

   posbox,vis=vis,color=clr

The map LIMIT is sometimes given in !map.ll_box. But beware of special cases:

If all the values are 0 do not use the LIMIT keyword to reset the scaling.
Also if !map.ll_box EQ [-90,!map.p0lon-180,90,!map.p0lon+180] do not use the LIMIT keyword to reset the scaling.
Special Note: If !map.p0lon-180 is LT -180 then add 360 to both !map.p0lon-180 and !map.p0lon+180 as a test limit. If !map.ll_box EQ that test limit then do not use LIMIT.

The SATELLITE projection can use an 8 element form of LIMIT. It is not known if or where these 8 values are saved so the user should keep track of that form of the LIMIT:
4 elements: [latmin, lonmin, latmax, lonmax]
8 elements: [latLeft, lonLeft, latTop, lonTop, LatRt, lonRt, LatBot, LonBot]

Finally, there are the 3 special case sets of parameters mentioned above: Standard Parallels, Satellite Parameters, and Ellipsoid Parameters. These special parameters are saved in !map.p in various elements. These cases are described below.

Standard Parallels

Standard Parallel or Parallels for Albers (proj=14) and Conic (proj=3).

The standard parallel(s) keyword is:

   standard_parallels=[lat1,lat2]

Array may be one or two elements.

Retrieve values as follows:

   std_p = !map.p[[3,4]]/!dtor

!map.p has standard parallels in radians. For the single parallel case both elements will be the same.

Note: LIMIT must be handled differently for these projections. See notes below on this.

Satellite Projection

This is a somewhat complicated projection system (to use, easy to understand). A couple notes: the keyword /USA seems to not always work unless the 8 element LIMIT option is used (actually it mostly does). Also !MAP.LL_BOX will still contain only 4 elements, not sure where the limit information is saved if anywhere. !MAP.LL_BOX appears to give a valid 4 element limit, but if used in place of the original 8 element limit the map is different (but similar). So unless a location with the 8 elements is found this case must be handled by the user.

The satellite parameters are:

  sat_p = [F, Omega, Gamma]
  F = Dist in radii from earth center.
  Omega = downward tilt of camera.
  Gamma = rotation angle (put in ang0 by map_set).

The values may be retrieved as follows:

  sat_p = [ !map.p[0], !map.p[1]/!dtor, 0.]

Note, last 0. is view rotation, but !map.rotation is set to sat_p[2]+ang0 and is exactly the same, so use sat_p view rot=0 and all rot in ang0.

Ellipsoid

The map_set keyword is

   ellips = [a, e^2, k]

where a=equatorial radius, in meters, e^2: eccentricity squared, k: scale on the central meridian.
The default is ellips1 = [6378206.4d0, 0.00676866d0, 0.9996d0] (Clarke 1866)

Ellipsoids are allowed in two cases:

For Transverse Mercator (proj=15):

  !map.p=[ellips1[2],ellips1[1],e_2,ellips1[0],ellips1[1],0., 0.]
  So ellips = [ !map.p[3], !map.p[1], !map.p[0]]

For Lambert Conic Ellipsoid (proj will be 18. Extract ellipsoid and call with proj=3 but use ellipsoid):

  !map.p = [n, F, r0, s, ellips1]  (s=std par=[lt1,lt2] ellips=[e1,e2,e3]).
  So ellips = [ !map.p[5], !map.p[6], !map.p[7] ]

Using the extracted projection information

The information needed to set up the map projection at a later time can be stored somewhere and applied to a map image after it has been loaded. For non-lossy image formats this information may be embedded in the image pixels as ascii text (each character as a byte pixel). It is useful to include a magic number that flags the information as available. Also it is useful to include the LIMIT as the 8 element case, but the last 4 elements will be 0 for the 4 element case. The special case parameters can be included from !map.p elements 0 through 7. When the projection information is interpreted the special case values can be extracted based on the projection code found in the saved info.

The /NOBORDER keyword is one problem. It causes a scale difference of 1% so must be handled, but there is no clue in any of the system variables to tell if /NOBORDER was used in the map_set call. A new routine has been written, posbox, to look at the map on the screen and determine if /NOBORDER was used. It also returns the color of the map border. This gives the information on whether to use /NOBORDER for a later map_set call, and if so, what color to use so it won't change the map when the border is overplotted.

The LIMIT=... keyword is another problem. In many cases !map.ll_box gives the limits that were used in the map_set call. If no limit was used !map.ll_box is sometimes [0,0,0,0], but sometimes !map.ll_box = [-90,!map.p0lon-180,90,!map.p0lon+180] Both case indicate that LIMIT=... should not be used in the map_set call used to set up the scaling at a later time. Finally, no way is known to automatically deal with the 8 element case for the limit. This case is handled by the user giving the values to the map_put_scale routine below.

IDL map scaling Routines

Two IDL routines have been created:

map_put_scale extracts the needed info from the system variables and embeds it in the current window.

map_set_scale uses that info to reset the map scaling.
Still under construction.

Some notes for special cases

/CONIC and /ALBERS

The old form for this projection for /conic was:

   map_set,30,40,-90,/conic,/cont

where 30, 40 are the standard parallels, and -90 is the rotation angle which is the central longitude. All the values default to 0. Note the central lat is 30 deg in this case. LIMIT should override this. The equivalent new form is:

   map_set,30,-90,/conic,/cont,standard_par=[30,40]

This brings the map_set call back to standard form with lat0, lon0 (ang always 0). Last positional arg is ignored (can use 0). Note the central lat need no longer be the same as the first parallel.

Note: Central long and lat may be plotted as:

 ver,/norm,midv(!x.window) & hor,/norm,midv(!y.window)   ; Map center.
 plots,lon0,lat0,psym=2   ; Central long,lat.

LIMIT is a problem for /conic and /albers

For this case the limit gets set to some weird value not detected by map_put_scale. For example:

  map_set,20,-90,/conic,/cont,standard_par=[30,40]
  print,!map.ll_box
    gives
  -58.133053       14.589621       84.663124       357.84503

but no limit was used at all (exact same limit for /albers). How can that be detected???

map_set,20,-90,/conic,/cont,standard_par=[30,40],limit=[25,-120,55,-60]

works ok. So only the case where LIMIT was not used must be detected by map_put_scale somehow.

Here is a possible method to detect this case. When LIMIT is not used for the /CONIC or /ALBERS projection it appears that !map.uv_box always is [-2.0400000, -2.0400000, 2.0400000, 2.0400000]. This seems to hold if STANDARD_PARALLELS=[p1,p2] are used or not, and for the central point in any hemisphere. It does not seem to happen at all when LIMIT is used. Not sure about other projections.

---

Limit seems to now work for the special cases. That is handled in map_put_scale, not map_set_scale.

  Proj:   3, 14  Set up STD_P
  Proj:   7      Set up SAT_P
  Proj:  15, 18  Set up ELLIPS (for 18 use 3 with ELLIPS).

map_set_scale needs completed. Add all projections, also add list and return of all values (like set_scale).

---

It turns out there are a number of special cases where a LIMIT appears in !map.ll_box but /LIMIT must not be used. So far it seems these cases can be detected using !map.uv_box, testing it for special values depending on the map projection used. /NOBORDER seems to make a big difference in this area. For example:

t_idl> map_set,/cont,40,-90,30,/Stereographic          
t_idl> print,!map.uv_box,!map.ll_box,!map.projection
      -1.0200000      -1.0200000       1.0200000       1.0200000
       0.0000000       0.0000000       0.0000000       0.0000000
           1

t_idl> map_set,/cont,40,-90,30,/Stereographic,/noborder
t_idl> print,!map.uv_box,!map.ll_box,!map.projection   
      -1.0000000      -1.0000000       1.0000000       1.0000000
      -41.560764      -180.00000       90.000000       180.00000
           1

The map size is only changed by 1%. It is not at all clear why !map.ll_box changes so much.

Special cases:

• Proj 1: !map.uv_box = [-1,-1,1,1]
• Proj 2: !map.uv_box = [-1,-1,1,1]
• Proj 3: !map.uv_box = [-2,-2,2,2]
• Proj 4: !map.uv_box = [-2,-2,2,2]
• Proj 5: !map.uv_box = [-2,-2,2,2]
• Proj 6: Use Lon0 +/- 180 (already works).
• Proj 7: ?
• Proj 8: Use Lon0 +/- 180 (already works).
• Proj 9: Use Lon0 +/- 180 (already works).
• Proj 10: Use Lon0 +/- 180 (already works).
• Proj 11: Use Lon0 +/- 180 (already works).
• Proj 12: Use Lon0 +/- 180 (already works).
• Proj 13: Use Lon0 +/- 180 (already works).
• Proj 14: Works.
• Proj 15: Works.
• Proj 16: Works.

Broken cases

Problem case: /MILLER if /NOBORDER and /HORIZON used. since horizon plots in same place as border would.

Major problem

Use of the SCALE=sc keyword in map_set does not work with map_put_scale and map_set_scale. Not sure why yet.

Other cases should work. Expand map_text_scale.pro.