MAP_CONTOUR Procedure
Draws a contour plot from longitude/latitude data stored in a 2D array.
Usage
MAP_CONTOUR, z[, x, y]
Input Parameters
z—A two-dimensional array containing the longitude/latitude values that make up the contour surface.
x—(optional) A vector specifying the longitude coordinates for the contour surface.
y—(optional) A vector specifying the latitude coordinates for the contour surface.
Keywords
Cartesian—When this keyword is set, MAP_CONTOUR first projects the grid of z values, grids the resulting projected points to a regular grid and performs the contouring in cartesian (data coordinate) space. This keyword allows you to get labeled contours and contours that can encircle a pole in a polar map. (Default: off). If Cartesian=1, you must also use the C_Labels keyword to get labeled contours.
If Cartesian=1 then FAST_GRID3 does the gridding, and if Cartesian=[nx, ny, r] then GRIDN does the gridding on a (nx, ny) grid with weighting function of order r (see GRIDN Function). Optionally, you can access the GRIDN t keyword with Cartesian=[nx, ny, r, t].
Cartesian has a third gridding option available only for those users with access to PV-WAVE IMSL Mathematics: Setting the Cartesian keyword equal to a PV-WAVE list variable invokes the RBF interpolator RBFIMSCL (see BAR3D Procedure) for the gridding step. This list must have three or four elements: The first two elements are integers designating the grid resolution. The next element is also an integer, with a nonzero value indicating that the RBF interpolant should include a linear term. The optional element is the name of a user-supplied basis function.
 
note
Performance Notes: If you use GRIDN, be aware that choosing an r value that is a FLOAT or an odd integer forces PV-WAVE to use floating point exponentiation.
To use the substantially faster integer exponentiation, set r to an even integer.
For larger data sets, the optional RBFIMSCL Procedure in the PV-WAVE IMSL Mathematics toolkit offers much better performance than GRIDN. RBFIMSCL is only available if you have PV-WAVE Advantage.
 
note
RUNNING PV-WAVE OVER A NETWORK: See the Fnam keyword.
Fnam—The name of a scratch file used by the third Cartesian gridding option. Defaults to rbfintrp.temp. For large datasets it is more efficient if Fnam refers to a file on a local disk. By default the scratch file is written to the PV-WAVE home directory. If this directory is on a disk remote from the machine running PV-WAVE and Fnam is not used, then scratch file I/O can be very time consuming for large datasets. In a situation like this, it is very important that Fnam be given the full path to some file on a disk local to the machine running PV-WAVE. Two instances of the gridding algorithm cannot simultaneously use the same scratch file, so unique values of Fnam should be ensured if the possibility of such a conflict exists. On exit the scratch file contents are erased, but the 0 byte file is not deleted.
File_Path—(string) The file pathname for a temporary file used by MAP_CONTOUR. This keyword is optional. See the Discussion section for more information.
Filled—If present and nonzero, fills the contours.
Iter or Nghbr—These keywords are used in conjunction with the Cartesian keyword and control parameters to the FAST_GRID3 function, which is used to perform gridding of the z grid data after applying the map projection. See the FAST_GRID3 documentation for a description of these parameters.
Pattern—A rectangular array of pixels giving the fill pattern.
Void —The name of a user supplied procedure describing voids in the lat-lon space. When the Cartesian keyword has enabled GRIDN, Void enables the GRIDN F keyword.
Standard Plotting Keywords
The following standard plotting keywords are used with MAP_CONTOUR. See Chapter 21: Graphics and Plotting Keywords for more information.
 
 
note
The C_Annotation, C_Charsize, C_Labels, and Font keywords only work when the Cartesian keyword is also specified.
Discussion
If x and y are provided, the contour is plotted as a function of the longitude and latitude locations specified by x and y. Otherwise, the contour is assumed to occupy the entire range specified in the last call to the MAP procedure.
Each element of x specifies the longitude coordinate for a column of z. For example, X(0) specifies the longitude coordinate for Z(0, *).
Each element of y specifies the latitude coordinate for a row of z.
In addition to line plot overlays, the Filled and Pattern keywords allow contours to be filled.
If you are not using the Cartesian keyword and your contours cross the 180th meridian or include either pole, your contours may not be drawn correctly. To correct this the contours must all be closed. This is usually accomplished by padding the data with zeros or some other value outside the range of the data, as the following code fragment illustrates, where m and n represent the dimensions of the original array.
; Create a background array that is two elements larger than
; the original array, array1 
new_array = REPLICATE(MIN(array1), m+2, n+2)
; Insert original data into the new array. 
new_array(1, 1) = array1
 
note
You cannot use the Filled keyword to fill contours when the projection is set to 99 (Satellite View); however, it is possible to place filled contours on a satellite projection. To do this, create a 2D image using CONTOUR and CONTOURFILL, or by using the CONTOUR2 procedure with the Fill keyword, and then use the Image keyword with the MAP procedure to wrap the image on the globe.
The MAP_CONTOUR procedure creates a temporary file when it is called that by default is named wvctmp.dat and placed in your current directory. You can override this filename and path using the File_Path keyword with MAP_CONTOUR. This allows you to specify a string containing the path to a temporary file. This file is deleted when MAP_CONTOUR is finished plotting.
Example 1
In this simple example, a 2D array of data is created and superimposed as contours on a map projection.
; Create a 2D array of data.
data = HANNING(20,20)
MAP, Projection=4, Range=[-150, 30, 30, 70]
MAP_CONTOUR, data, NLevels=10
Example 2
TEK_COLOR
c = WoColorConvert([1,2,3,4,5,6])
l = INDGEN(5)
z = DIST(40)
WINDOW, 0
MAP, Range=[-180,-90,180,90], /Exact, Projection=1, Color=c(0)
MAP_CONTOUR, z, Cartesian=[40,40,5], C_Colors=c(1:*), C_Labels=l
Example 3
This example demonstrates the use of the Void keyword.
PRO VOIDS, p, i, v
   ; void function used by GRIDN 
   COMMON data, lon, lat, datau, u, void
   ; 2d array of all lat-lon pairs
   t = TRANSPOSE(CPROD(LIST(lon, lat)))
   ; project these pairs
   MAP_PROJ, t 
   t = TRANSPOSE(t)
   x = TENSOR_SUB(t(*,0), p(*,0))
   y = TENSOR_SUB(t(*,1), p(*,1))
   z = MIN(x*x+y*y, D=0, i)
   ; for each grid point in p, find the index of closest data
   ; point in t 
   i = i MOD N_ELEMENTS(datau)
   i = WHERE(datau(i) NE void)
   v = u
END
 
PRO MAP_CONTOUR_VOIDS
   COMMON data, lon, lat, datau, u, void
   lon = INTERPOL([-180,180], 64)
   lat = INTERPOL([-90,90], 32)
   datau = DIST(64, 32)
   i = [0L]
   file = GETENV('RW_DIR') + $
      '/mapping-2_0/data/map_contour_voids.dat'
   s = DC_READ_FREE(file, i, /Column, Resize=[1])
   void = -9e30
   datau(i) = void
   v = WHERE(datau NE void)
   v = MIN(datau(v), Max=w)
   u = v - 10*(w-v)
   lv = INTERPOL([v,w], 10)
   MAP, Projection=1, /Exact, /Nodata
   MAP_CONTOUR, datau, lon, lat, Levels=lv, $
      Cartesian=[40,40,5], Void='voids' 
END
Example 4
This example demonstrates the RBF gridding option associated with the Cartesian keyword. It also shows how to mask any contour lines external to the map border and how to fill land masses with a solid color.
; Define a map range and a grid of points on that range
nx = 241
ny = 253
range = [118, 18, 150, 50]
lon = INTERPOL(range([0,2]), nx)
lat = INTERPOL(range([1,3]), ny)
 
; Define some data on the grid (1/4 of a hanning surface)
a = (HANNING(2*nx, 2*ny))(0:nx-1, 0:ny-1)
 
; Define 11 equally spaced contour levels
nv = 11
level = INTERPOL([MIN(a, Max=max), max], nv)
 
; Define two RGB triplets
color = [[0,0,255], [200,200,255]]
 
; Define RGB triplets for ocean color between each contour level
color = INTRP(color, 1, INTERPOL([0,1], nv+1))
 
; Convert the RGB triplets to longs representing 24 bit colors
color = TRANSPOSE(NINT(color)) # 256L^[0, 1, 2]
 
; Draw just the map boundary
pos = [0.1, 0.1, 0.9, 0.9]
MAP, Range=range, Projection=8, /Exact, Bound=2, Pos=pos, Col=0
 
; Find the regions inside and outside the boundary
mask = BLOBCOUNT(BYTSCL(TVRD(), Top=1), [0, 0])
 
; The larger region defines the pixels outside the boundary
junk = MAX([N_ELEMENTS(mask(0)), N_ELEMENTS(mask(1))], ibig)
mask = TRANSPOSE(mask(ibig))
 
; Define surface cover, and then extract the land pixels
covr = MAP_LAND(range, /Fast)
land = WHERE(covr(*, 2))
 
; Draw map contours, using the third Cartesian gridding option
MAP_CONTOUR, a, lon, lat, Level=level, /Filled, $
   C_Color=color, C_Labels=REPLICATE(1, nv), $
   Color=0, Cartesian=LIST(100, 100, 0)
 
; Draw the map, and mask any contour lines outside the boundary
MAP, Range=range, Projection=8, /Exact, Boundary=2, $
   Position=pos, Color='50aaff'x, Thick=2, /Noerase
PLOTS, mask, /Device, Psym=3, Color=0
 
; Fill the land with a solid color
MAP_PLOTS, covr(land,0), covr(land,1), Psym=3, Color='50aaff'x
See Also