STARTING TO PROGRAM IN IDL


Written by Virginia Bluth

Last updated: April 4, 2003



Writing basic IDL programs

In this first part, you will start writing simple IDL programs. If you have any experience with Fortran or C, you are that much ahead of the game, but if you do not, you will catch on quickly enough. IDL is very user-friendly. It has a very useful and thorough help file. And there are some very useful links on the web, most notably, Coyote's Guide to IDL Programming at

http://www.dfanning.com

or try the links at

http://idlastro.gsfc.nasa.gov/other_url.html

or just try typing in IDL at google.com.


You will try the interactive mode and the programming mode of IDL. You will write a program which uses the building blocks of imperative program - commands, IF-statements, DO-loops, FOR-loops, etc.. Mostly, you should try to have fun with all of the neat options that IDL allows you to use in your plots.


To Get Started:

The important thing to remember is that there is nothing horrendous that you can do the computer as long as you stop yourself before you hit the computer with something hard. Probably the worst thing you could do is delete the entire program you just spent 8 hours writing and debugging. This can be easily avoided by getting into the habit of saving your work regularly (like every 15 minutes - 5 minutes is usually recommended, but probably not worth the effort). You might also want to keep a backup. You can do this most easily by typing:

cp lastversion.pro lastversion.copy

at the screen prompt at the beginning of every session. ('lastversion.pro' would be the current name of the program. 'Lastversion.copy' would be the new name - I would eliminate the pro suffix in your copy, since it will not compile if the filename is different from the program name. More on that later.)

So go ahead and try things. Find things on the web and paste them into a program and see how they work. Try typing in a command you think might work. I also recommend that when you are trying out something new, try it in interactive mode, or as a short program. Later when you are sure it works, you can add it to your full program.


So, here goes...

1. At unix prompt, type

idl -c

which puts you in the IDL Developmental Environment.

2. If you ever need the IDL help, on the IDL prompt line, type

?

wait a minute and it will boot up. I find that I can find things pretty fast by going

to the index then doing a search (try typing in a word you think might be the name of a command). For now, try the help file by typing in contour. You will see a list of information topics. The top one will be Contour. Click on describe. You should be overwhelmed with more information than you could possibly want. Note that the help file has examples at the bottom of the description. You will need to scroll way down to see them, but they are usually worth the effort. For now, exit help.

3. To run an IDL program:

under the file command, click on save program

under run command, click on compile program

under run command, click on run program


4. When you are done, to leave IDL, at IDL prompt type

exit

or click under file (in lefthand corner) on exit.



A few things to know which might be helpful:

1. IDL is case insensitive, except for commands that interact with the operating system and when it is performing string comparisons. This means that the word,

else = ELSE = ElSe = eLsE = ...etc.

BUT "else" does not equal "ELSE" because words inside of the quotes are strings.


2. Variable names must start with a letter. They can include other letters, digits,

underscore characters, and dollar signs. A variable name may have up to 255

characters.


SYNTAX :

the way programs and commands must be written in order for them to work


For the most part, your programs will consist of a sequence of commands. There is a format which must be followed in order for IDL (I mean the IDL compiler, but I'm just going to personify IDL to cut down on typing) to understand what you are trying to do.


Here is a line of code which you might use:


Contour, peak, lon, lat, Xstyle=1, Ystyle=1, /Follow, $

Levels=vals, C_Labels=[1,0,1,0, 0, 1, 1, 0] ;


The first word, Contour, is the name of a command which IDL recognizes, in other words, Contour (or contour, or CONTOUR, etc.) is in IDL's library as long as it is in the correct format. If you look up the contour command in the IDL help file1, you will find that IDL recognizes the contour command when it is in one of these 2 formats:

Contour, z

or

Contour, z, x, y

Our command follows the second format (Below I will show how we know it's the 2nd and not the first format).


The command is broken down as follows:


Contour is the name of the command. This command tells IDL to make a contour plot using the following data (z, or z, x and y). A command name must be spelled out in its entirety (so the command,

cont, peak

will not be recognized by IDL as the contour command because the word contour does not appear).


The next 3 words, peak, lon, and lat are called positional parameters because they must appear in an exact order. Peak, lon, and lat would be 3 variables which the programmer defined before typing the contour command. So, their names could really be anything. When IDL sees the command, contour, it will know to check if you defined them and if they are of the correct type. If you did not define them, or they are not of the type that is expected for this command then your program will not compile and you will get an error message to that effect.


Because there are 2 defined formats for contour, IDL knows to check if there is 1 variable after the word, contour, or 3. Any other number of variables after contour will result in an error message.


IDL expects the first positional parameter to be either:

1. a two-dimensional array with evenly spaced data which ranges from the indices of z (= peak), or

2. a two-dimensional array which is dependent on the data in the next 2 positional parameter x and y, named lon and lat in our command.


You must position the name of the array with the z values first, then the x values, then the y, or the plot will not look right, or the program might not even run.


All of the terms following the positional parameters are called keyword parameters. Keyword parameters are optional. They can appear in any order, even in the middle of the positional parameters (the compiler will figure out which is which and not get confused). Keyword parameters have special names which IDL recognizes. To know which keyword parameters (and their specific values) can be used for a given command, it is best to consult the help file. Again, look at all the different keyword parameters which can be used with contour. These parameters make your plot look a lot nicer. You can also sop up a huge amount of time playing with them to make your plot look exactly the way you want it to.


Notes on syntax:

/Keyword is identical to (means the same as) the syntax Keyword=1.


So could have been written:

Contour, peak, Levels=vals, lon, /Xstyle, /Ystyle, $

/Follow, Levels=vals, lat, C_Labels=[1,0,1,0,0,1,1,0]

and still get the same contour plot.



The Contour help file:


The CONTOUR procedure draws a contour plot from data stored in a rectangular array or from a set of unstructured points. Both line contours and filled contour plots can be created. Note that outline and fill contours cannot be drawn at the same time. To create a contour plot with both filled contours and outlines, first create the filled contour plot, then add the outline contours by calling CONTOUR a second time with the OVERPLOT keyword.


Contours can be smoothed by using the MIN_CURVE_SURF function on the contour data before contouring.


Using various keywords, described below, it is possible to specify contour levels, labeling, colors, line styles, and other options. CONTOUR draws contours by searching for each contour line and then following the line until it reaches a boundary or closes.


Smoothing Contours


The MIN_CURVE_SURF function can be used to smoothly interpolate both regularly and irregularly sampled surfaces before contouring. This function replaces the older SPLINE keyword to CONTOUR, which was inaccurate and is no longer supported. MIN_CURVE_SURF interpolates the entire surface to a relatively fine grid before drawing the contours.


Syntax


CONTOUR, Z [, X, Y] [, C_CHARSIZE=value] [, C_CHARTHICK=integer] [, C_COLORS=vector] [, C_LABELS=vector{each element 0 or 1}] [, C_LINESTYLE=vector] [{, /FILL | , /CELL_FILL} | [, C_ANNOTATION=vector_of_strings] [, C_ORIENTATION=degrees] [, C_SPACING=value]] [, C_THICK=vector] [, /CLOSED] [, /DOWNHILL] [, /FOLLOW] [, /IRREGULAR] [, LEVELS=vector] [, NLEVELS=integer{1 to 60}] [, MAX_VALUE=value] [, MIN_VALUE=value] [, /OVERPLOT] [{, /PATH_DATA_COORDS, PATH_FILENAME=string, PATH_INFO=variable, PATH_XY=variable} | , TRIANGULATION=variable] [, /PATH_DOUBLE] [, /XLOG] [, /YLOG] [, /ZAXIS] [, /ZLOG]


Graphics Keywords: Accepts all graphics keywords accepted by PLOT except for: LINESTYLE, PSYM, SYMSIZE. See Graphics KeywordsAccepted.


Arguments


Z


A one- or two-dimensional array containing the values that make up the contour surface. If arguments X and Y are provided, the contour is plotted as a function of the (X, Y) locations specified by their contents. Otherwise, the contour is generated as a function of the two-dimensional array index of each element of Z.


If the IRREGULAR keyword is set, X, Y, and Z are treated as vectors. Each point has a value of Z(i) and a location of (X(i), Y(i))


This argument is converted to double-precision floating-point before plotting. Plots created with CONTOUR are limited to the range and precision of double-precision floating-point values.


X


A vector or two-dimensional array specifying the X coordinates for the contour surface. If X is a vector, each element of X specifies the X coordinate for a column of Z (e.g., X[0] specifies the X coordinate for Z[0,*]). If X is a two-dimensional array, each element of X specifies the X coordinate of the corresponding point in Z (i.e., Xij specifies the X coordinate for Zij).


Y


A vector or two-dimensional array specifying the Y coordinates for the contour surface. If Y a vector, each element of Y specifies the Y coordinate for a row of Z (e.g., Y[0] specifies the Y coordinate for Z[*,0]). If Y is a two-dimensional array, each element of Y specifies the Y coordinate of the corresponding point in Z (Yij specifies the Y coordinate for Zij).


Keywords


C_ANNOTATION


The label to be drawn on each contour. Usually, contours are labeled with their value. This parameter, a vector of strings, allows any text to be specified. The first label is used for the first contour drawn, and so forth. If the LEVELS keyword is specified, the elements of C_ANNOTATION correspond directly to the levels specified, otherwise, they correspond to the default levels chosen by the CONTOUR procedure. If there are more contour levels than elements in C_ANNOTATION, the remaining levels are labeled with their values.


Use of this keyword implies use of the FOLLOW keyword.


Note - This keyword has no effect if the FILL or CELL_FILL keyword is set (i.e., if the contours are drawn with solid-filled or line-filled polygons).


Example


To produce a contour plot with three levels labeled "low", "medium", and "high":


CONTOUR, Z, LEVELS = [0.0, 0.5, 1.0], $

C_ANNOTATION = ['low', 'medium', 'high']


C_CHARSIZE


The size of the characters used to annotate contour labels. Normally, contour labels are drawn at 3/4 of the size used for the axis labels (specified by the CHARSIZE keyword or !P.CHARSIZE system variable. This keyword allows the contour label size to be specified directly. Use of this keyword implies use of the FOLLOW keyword.


C_CHARTHICK


The thickness of the characters used to annotate contour labels. Set this keyword equal to an integer value specifying the line thickness of the vector drawn font characters. This keyword has no effect when used with the hardware drawn fonts. The default value is 1. Use of this keyword implies use of the FOLLOW keyword.


C_COLORS


The color index used to draw each contour. This parameter is a vector, converted to integer type if necessary. If there are more contour levels than elements in C_COLORS, the elements of the color vector are cyclically repeated.


Example


If C_COLORS contains three elements, and there are seven contour levels to be drawn, the colors c0, c1, c2, c0, c1, c2, c0 will be used for the seven levels. To call CONTOUR and set the colors to [100,150,200], use the command:


CONTOUR, Z, C_COLORS = [100,150,200]


C_LABELS


Specifies which contour levels should be labeled. By default, every other contour level is labeled. C_LABELS allows you to override this default and explicitly specify the levels to label. This parameter is a vector, converted to integer type if necessary. If the LEVELS keyword is specified, the elements of C_LABELS correspond directly to the levels specified, otherwise, they correspond to the default levels chosen by the CONTOUR procedure. Setting an element of the vector to zero causes that contour label to not be labeled. A nonzero value forces labeling.


Use of this keyword implies use of the FOLLOW keyword.


Example


To produce a contour plot with four levels where all but the third level is labeled:


CONTOUR, Z, LEVELS = [0.0, 0.25, 0.75, 1.0], $

C_LABELS = [1, 1, 0, 1]



C_LINESTYLE


The line style used to draw each contour. As with C_COLORS, C_LINESTYLE is a vector of line style indices. If there are more contour levels than line styles, the line styles are cyclically repeated. See LINESTYLE for a list of available styles.


Note - The cell drawing contouring algorithm draws all the contours in each cell, rather than following contours. Since an entire contour is not drawn as a single operation, the appearance of the more complicated linestyles will suffer. Use of the contour following method (selected with the FOLLOW keyword) will give better looking results in such cases.


Example


To produce a contour plot, with the contour levels directly specified in a vector V, with all negative contours drawn with dotted lines, and with positive levels in solid lines:


CONTOUR, Z, LEVELS = V, C_LINESTYLE = (V LT 0.0)


C_ORIENTATION


If the FILL keyword is set, this keyword can be set to the angle, in degrees counterclockwise from the horizontal, of the lines used to fill contours. If neither C_ORIENTATION nor C_SPACING are specified, the contours are solid filled.


C_SPACING


If the FILL keyword is set, this keyword can be used to control the distance, in centimeters, between the lines used to fill the contours.


C_THICK


The line used to draw each contour level. As with C_COLORS, C_THICK is a vector of line thickness values, although the values are floating point. If there are more contours than thickness elements, elements are repeated. If omitted, the overall line thickness specified by the THICK keyword parameter or !P.THICK is used for all contours.


CELL_FILL


Set this keyword to produce a filled contour plot using a "cell filling" algorithm. Use this keyword instead of FILL when you are drawing filled contours over a map, when you have missing data, or when contours that extend off the edges of the contour plot. CELL_FILL is less efficient than FILL because it makes one or more polygons for each data cell. It also gives poor results when used with patterned (line) fills, because each cell is assigned its own pattern. Otherwise, this keyword operates identically to the FILL keyword, described below.


Tip - In order for CONTOUR to fill the contours properly when using a map projection, the X and Y arrays (if supplied) must be arranged in increasing order. This ensures that the polygons generated will be in counterclockwise order, as required by the mapping graphics pipeline.


Warning - Do not draw filled contours over the poles on Cylindrical map projections. In this case, the polar points map to lines on the map, and the interpolation becomes ambiguous, causing errors in filling. One possible work-around is to limit the latitudes to the range of -89.9 degrees to + 89.9 degrees, avoiding the poles.


CLOSED


Set this keyword to a nonzero value to close contours that intersect the plot boundaries. After a contour hits a boundary, it follows the plot boundary until it connects with its other boundary intersection. Set CLOSED=0 along with PATH_INFO and/or PATH_XY to return path information for contours that are not closed.


DOWNHILL


Set this keyword to label each contour with short, perpendicular tick marks that point in the "downhill" direction, making the direction of the grade readily apparent. If this keyword is set, the contour following method is used in drawing the contours.


FILL


Set this keyword to produce a filled contour plot. The contours are filled with solid or line-filled polygons. For solid polygons, use the C_COLOR keyword to specify the color index of the polygons for each contour level. For line fills, use C_ORIENTATION, C_SPACING, C_COLOR, C_LINESTYLE, and/or C_THICK to specify attributes for the lines.


If the current device is not a pen plotter, each polygon is erased to the background color before the fill lines are drawn, to avoid superimposing one pattern over another.


Contours that are not closed can not be filled because their interior and exterior are undefined. Contours created from data sets with missing data may not be closed; many map projections can also produce contours that are not closed. Filled contours should not be used in these cases.


Note - If the current graphics device is the Z-buffer, the algorithm used when the FILL keyword is specified will not work when a Z value is also specified with the graphics keyword ZVALUE. In this situation, use the CELL_FILL keyword instead of the FILL keyword.


FOLLOW


In IDL version 5, CONTOUR always uses a line-following method. The FOLLOW keyword remains available for compatibility with existing code, but is no longer necessary. As in previous versions of IDL, setting FOLLOW will cause CONTOUR to draw contour labels.


IRREGULAR


Set this keyword to indicate that the input data is irregularly gridded. Setting IRREGULAR is the same as performing an explicit triangulation. That is:


CONTOUR, Z, X, Y, /IRREGULAR


is the same as


TRIANGULATE, X, Y, tri ;Get triangulation

CONTOUR, Z, X, Y, TRIANGULATION=tri


ISOTROPIC


Set this keyword to force the scaling of the X and Y axes to be equal.


Note - The X and Y axes will be scaled isotropically and then fit within the rectangle defined by the POSITION keyword; one of the axes may be shortened. See POSITION for more information.


LEVELS


Specifies a vector containing the contour levels drawn by the CONTOUR procedure. A contour is drawn at each level in LEVELS.


Example


To draw a contour plot with levels at 1, 100, 1000, and 10000:


CONTOUR, Z, LEVELS = [1, 100, 1000, 10000]


To draw a contour plot with levels at 50, 60, ..., 90, 100:


CONTOUR, Z, LEVELS = FINDGEN(6) * 10 + 50


MAX_VALUE


Data points with values above this value are ignored (i.e., treated as missing data) when contouring. Cells containing one or more corners with values above MAX_VALUE will have no contours drawn through them. Note that the IEEE floating-point value NaN is also treated as missing data. (See Special Floating-Point Values for more information on IEEE floating-point values.)


MIN_VALUE


Data points with values less than this value are ignored (i.e., treated as missing data) when contouring. Cells containing one or more corners with values below MIN_VALUE will have no contours drawn through them. Note that the IEEE floating-point value NaN is also treated as missing data. (See Special Floating-Point Values for more information on IEEE floating-point values.)


NLEVELS


The number of equally spaced contour levels that are produced by CONTOUR. If the LEVELS parameter, which explicitly specifies the value of the contour levels, is present, this keyword has no effect. If neither parameter is present, approximately six levels are drawn. NLEVELS should be a positive integer.


OVERPLOT


Set this keyword to make CONTOUR "overplot". That is, the current graphics screen is not erased, no axes are drawn and the previously established scaling remains in effect. You must explicitly specify either the values of the contour levels or the number of levels (via the NLEVELS keyword) when using this option, unless geographic mapping coordinates are in effect.


PATH_DATA_COORDS


Set this keyword to cause the output contour positions to be measured in data units rather than the default normalized units. This keyword is useful only if the PATH_XY or PATH_FILENAME keywords are set.


PATH_DOUBLE


Set this keyword to indicate that the PATH_FILENAME, PATH_INFO, and PATH_XY keywords should return vertex and contour value information as double-precision floating-point values. The default is to return this information as single-precision floating-point values.


PATH_FILENAME


Specifies the name of a file to contain the contour positions. If PATH_FILENAME is present, CONTOUR does not draw the contours, but rather, opens the specified file and writes the coordinates of the contours, into it. The file consists of a series of logical records containing binary data. Each record is preceded with a header structure defining the contour as follows:


If the PATH_DOUBLE keyword is not set:


{CONTOUR_HEADER, TYPE:0B, HIGH:0B, LEVEL:0, NUM:0L, VALUE:0.0}


If the PATH_DOUBLE keyword is set:


{CONTOUR_DBL_HEADER, TYPE:0B, HIGH:0B, LEVEL:0, NUM:0L, VALUE:0.0D}


The fields are: Field Description

TYPE A byte that is zero if the contour is open, and one if it is closed.

HIGH A byte that is 1 if the contour is closed and above its surroundings, and is 0 if the contour is below. This field is meaningless if the contour is not closed.

LEVEL A short integer with value greater or equal to zero (It is an index into the LEVELS array).

NUM The longword number of data points in the contour.

VALUE The contour value. If the PATH_DOUBLE keyword is not set, this is a single-precision floating-point value; if the PATH_DOUBLE keyword is set, this is a double-precision floating-point value.


Following the header in each record are NUM X-coordinate values followed by NUM Y-coordinate values. By default, these values are specified in normalized coordinates unless the PATH_DATA_COORDS keyword is set.


PATH_INFO


Set this keyword to a named variable that will return path information for the contours. This information can be used, along with data stored in a variable named by the PATH_XY keyword, to trace closed contours. To get PATH_INFO and PATH_XY with contours that are not closed, set the CLOSED keyword to 0. If PATH_INFO is present, CONTOUR does not draw the contours, but rather records the path information in an array of structures of the following type:


If the PATH_DOUBLE keyword is not set:


{CONTOUR_PATH_STRUCTURE, TYPE:0B, HIGH_LOW:0B, $

LEVEL:0, N:0L, OFFSET:0L, VALUE:0.0}


If the PATH_DOUBLE keyword is set:


{COUNTOUR_DBL_PATH_STRUCTURE, TYPE:0B, HIGH_LOW:0B, LEVEL:0,

N: 0L, OFFSET:0L, VALUE:0.0D}


The fields are: Field Description

TYPE A byte that is zero if the contour is open, and one if it is closed. In the present implementation, all contours are closed.

HIGH_LOW A byte that is 1 if the contour is above its surroundings, and is 0 if the contour is below.

LEVEL A short integer indicating the index of the contour level, from zero to the number of levels minus one.

N A long integer indicating the number of XY pairs in the contour's path.

OFFSET A long integer that is the offset into the array defined by PATH_XY, representing the first XY coordinate for this contour.

VALUE The contour value. If the PATH_DOUBLE keyword is not set, this is a single-precision floating-point value; if the PATH_DOUBLE keyword is set, this is a double-precision floating-point value.


See the examples section below for an example using the PATH_INFO and PATH_XY keywords to return contour path information.


PATH_XY


Set this keyword to a named variable that returns the coordinates of a set of closed polygons defining the closed paths of the contours. This information can be used, along with data stored in a variable named by the PATH_INFO keyword, to trace closed contours. To get PATH_XY and PATH_INFO with contours that are not closed, set the CLOSED keyword to 0. If PATH_XY is present, CONTOUR does not draw the contours, but rather records the path coordinates in the named array. If the PATH_DOUBLE keyword is not set, the array will contain single-precision floating point values; if the PATH_DOUBLE keyword is set, the array will contain double-precision floating point values. By default, the values in the array are specified in normalized coordinates unless the PATH_DATA_COORDS keyword is set.


See the examples section below for an example using the PATH_INFO and PATH_XY keywords to return contour path information.


TRIANGULATION


Set this keyword to a variable that contains an array of triangles returned from the TRIANGULATE procedure. Providing triangulation data allows you to contour irregularly gridded data directly, without gridding.


XLOG


Set this keyword to specify a logarithmic X axis.


YLOG


Set this keyword to specify a logarithmic Y axis.


ZAXIS


Set this keyword to draw a Z axis for the CONTOUR plot. CONTOUR draws no Z axis by default. This keyword is of use only if a three-dimensional transformation is established.


Graphics Keywords Accepted


See Graphics Keywords for the description of graphics and plotting keywords not listed above.


BACKGROUND, CHARSIZE, CHARTHICK, CLIP, COLOR, DATA, DEVICE, FONT, NOCLIP, NODATA, NOERASE, NORMAL, POSITION, SUBTITLE, T3D, THICK, TICKLEN, TITLE, [XYZ]CHARSIZE, [XYZ]GRIDSTYLE, [XYZ]MARGIN, [XYZ]MINOR, [XYZ]RANGE, [XYZ]STYLE, [XYZ]THICK, [XYZ]TICKFORMAT, [XYZ]TICKINTERVAL, [XYZ]TICKLAYOUT, [XYZ]TICKLEN, [XYZ]TICKNAME, [XYZ]TICKS, [XYZ]TICKUNITS, [XYZ]TICKV, [XYZ]TICK_GET, [XYZ]TITLE, ZVALUE.


Examples


Example 1


This example creates a contour plot with 10 contour levels where every other contour is labeled:


;Create a simple dataset to plot:

Z = DIST(100)

;Draw the plot:

CONTOUR, Z, NLEVELS=10, /FOLLOW, TITLE='Simple Contour Plot'


Example 2


This example shows the use of polygon filling and smoothing.


;Handle TrueColor displays:

DEVICE, DECOMPOSED=0

;Create a surface to contour (2D array of random numbers):

A = RANDOMU(seed, 5, 6)

;Smooth the dataset before contouring:

B = MIN_CURVE_SURF(A)

;Load discrete colors for contours:

TEK_COLOR

;Draw filled contours:

CONTOUR, B, /FILL, NLEVELS=5, C_COLOR=INDGEN(5)+2

;Overplot the contour lines with tickmarks:

CONTOUR, B, NLEVELS=5, /DOWNHILL, /OVERPLOT


Alternatively, we could draw line-filled contours by replacing the last two commands with:


CONTOUR, B, C_ORIENTATION=[0, 22, 45]

CONTOUR, B, /OVERPLOT, NLEVELS=5


Example 3


The following example saves the closed path information of a set of contours and plots the result:


; Create a 2D array of random numbers:

A = RANDOMU(seed, 8, 10)

; Smooth the dataset before contouring:

B = MIN_CURVE_SURF(A)

; Compute contour paths:

CONTOUR, B, PATH_XY=xy, PATH_INFO=info

FOR I = 0, (N_ELEMENTS(info) - 1 ) DO BEGIN

S = [INDGEN(info(I).N), 0]

; Plot the closed paths:

PLOTS, xy(*,INFO(I).OFFSET + S ), /NORM

ENDFOR


Example 4


This example contours irregularly-gridded data without having to call TRIGRID. First, use the TRIANGULATE procedure to get the Delaunay triangulation of your data, then pass the triangulation array to CONTOUR:


;Make 50 normal X, Y points:

x = RANDOMN(seed, 50)

y = RANDOMN(seed, 50)

;Make the Gaussian:

Z = EXP(-(x^2 + y^2))

;Get triangulation:

TRIANGULATE, X, Y, tri

;Draw the contours:

CONTOUR, Z, X, Y, TRIANGULATION = tri


See Also


IMAGE_CONT, SHADE_SURF, SHOW3, SURFACE




GENERAL NONRULES ON FORMATTING

These are rules which can make your program a lot more readable.

1. Use comments. These are for you, or anyone else who might want to understand your code. Anywhere the symbol, ; , is placed, IDL will ignore the rest of that line.

2. Any word which cannot be another word, e.g., a command name, or a keyword, is capitalized. This is not mandatory (hence, it's placement in the non-rules). IDL sure doesn't care, but gives an idea of the structure of your program at a glance.

3. Indent blocks of code (one tab or 4 spaces is standard). This is standard in all languages I have written in. Again, IDL does not care, but you will. It makes DO-loops, FOR-loops, IF-statements, etc. pop out. You will be able to quickly examine these blocks of code - and usually these will be where all your problems are. One of the most common mistakes in programming is to forget to close a block of code with an END statement. This is very easy to find, if you have indented.

Additional note: loops inside of loops, or expressions inside of loops, or expressions inside of expressions, etc. should be further indented. Don't worry, you won't go off the edge of the page.


Here's a program (actually a function, but not a huge difference. Functions return something, often a value, when they are done. Procedures do not, they just perform a task, for example, sorting an array ), with proper documentation, description, and spacing. A bit of a pain, but much better for others to read, and more helpful when you come back to it a few months later. I have added a few comments on what I have actually had to do or not do for some of my own code.


FUNCTION integrate_table, a, b

;

; NAME:

: INTEGRATE_TABLE

;

; PURPOSE:

; Integrate the area under a curve

;

; CATEGORY: optional

; Math

;

; CALLING SEQUENCE: how you would call this function in your program

; Result = INTEGRATE_TABLE(A, B)

;

; INPUTS:

; A: 1-D vector of independent variables (ie: x)

; May be any format, converted to FLOAT before operation

;

; B: 1-D vector of dependent variables (ie: f(x))

; May be any format, converted to FLOAT before operation

;

;

; OUTPUT:

; Floating point value representing area under the curve f(x)

;

; COMMON BLOCKS: optional

; None

;

; SIDE EFFECTS: sometimes your program will change something (like change the value of a variable)

; while it is doing what it is supposed to do. Other people using this function will want to know.

; None

;

; RESTRICTIONS:

; COMPLEX entries will be converted to FLOAT before operation

;

; MODIFICATION HISTORY: critical

; Written, VJR, March, 1995

;


area1 = 0.0

area2 = 0.0

x = FLOAT(a)

y = FLOAT(b)

final = N_ELEMENTS(x) - 1


FOR i=0, final-1 DO BEGIN

; add up

xdist = (x(i+1) - x(i))/2.0

area1 = area1 + (y(i)*xdist)


; add down

bck = final - i

xdist = (x(bck) - x(bck-1))/2.0

area2 = area2 + (y(bck)*xdist)

ENDFOR


This function returns a value

RETURN, area1+area2

END


Here's the function with everything that IDL does not need cut out.


function integrate_table, a, b

area1 = 0.0

area2 = 0.0

x = float(a)

y = float(b)

final = n_elements(x) - 1

for i=0, final-1 do begin

xdist = (x(i+1) - x(i))/2.0

area1 = area1 + (y(i)*xdist)

bck = final - i

xdist = (x(bck) - x(bck-1))/2.0

area2 = area2 + (y(bck)*xdist)

endfor

return, area1+area2

end



Basic information about IDL


1. The $ lets a user continue on the next line when the command gets too long.

2. Basic Data Types (defined in IDL)

From help file:


The IDL language is dynamically typed. This means that an operation on a variable can change that variable's type. In general, when variables of different types are combined in an expression, the result has the data type that yields the highest precision.

For example, if an integer variable is added to a floating-point variable, the result will be a floating-point variable.


In IDL there are twelve basic, atomic data types, each with its own form of constant. The data type assigned to a variable is determined either by the syntax used when creating the variable, or as a result of some operation that changes the type of the variable.


· Byte: An 8-bit unsigned integer ranging in value from 0 to 255. Pixels in images are commonly represented as byte data.

· Integer: A 16-bit signed integer ranging from -32,768 to +32,767. In IDL an integer is only 2 bytes long. So the biggest integer in IDL is 32,768. Therefore, in many loops is is better to use a long which is 4 bytes.

· Unsigned Integer: A 16-bit unsigned integer ranging from 0 to 65535.

· Long: A 32-bit signed integer ranging in value from approximately minus two billion to plus two billion.

· Unsigned Long: A 32-bit unsigned integer ranging in value from 0 to approximately four billion.

· 64-bit Long: A 64-bit signed integer ranging in value from -9,223,372,036,854,775,808 to +9,223,372,036,854,775,807.

· 64- bit Unsigned Long:A 64-bit unsigned integer ranging in value from 0 to 18,446,744,073,709,551,615.

· Floating-point: A 32-bit, single-precision, floating-point number in the range of ±10^38, with approximately six or seven decimal places of significance.

· Double-precision: A 64-bit, double-precision, floating-point number in the range of ±10^308 with approximately 14 decimal places of significance.

· Complex: A real-imaginary pair of single-precision, floating-point numbers. Complex numbers are useful for signal processing and frequency domain filtering.

·Double-precision complex: A real-imaginary pair of double-precision, floating-point numbers.

Note - In previous versions of IDL prior to version 4, the combination of a double-precision number and a complex number in an expression resulted in a single-precision complex number because those versions of IDL lacked the DCOMPLEX double-precision complex data type. Starting with IDL version 4, this combination results in a DCOMPLEX number.

· String: A sequence of characters, from 0 to 32,767 characters in length, which is interpreted as text.


Complex Data Types (these are for object-oriented programming/widgets)


· Structures: Aggregations of data of various types. Structures are discussed in Structures.

· Pointers: A reference to a dynamically-allocated heap variable. Pointers are discussed in Pointers.

· Object References: A reference to a special heap variable that contains an IDL object structure. Object references are discussed in Object Basics.


3. Arrays:

First a note about arrays (or vectors = 1-dimensional matrix or array). In the beginning, a user will need to generate a range of input values in order to generate a plot. To do this, choose a variable name, such as x, and assign it an array of floats or integers or whatever, which, by default, start at 0 and go up to whatever value the user assigns.

For example, to create an array of integers called arr from 0 to 500, a user would type the command

arr = indgen(500)

same thing for floats is

arr=findgen(500)

(To create arrays of the following types, for the range 0 to 500:

unsigned integers uindgen(500)

longs lindgen(500)

unsigned longs ulindgen(500)

doubles dindgen(500)

complex #s cindgen(500)

strings sindgen(500)

a few others on p 27 of text)

But what if want range to be less than 1 or 0?

Example:

To create an array from -500 to +200, type

x=findgen(700)-500

Example:

To create an array from 1E-5 to 1E-2, type

x=findgen(1000)+1 ; have now created an array from 1 to 1000

x=x/1E8 ; divide all values by 1E8


4. Fonts:

The default font for IDL is Roman, and it's pretty hard to read, so can type something like

!p.font=1

before setting titles, and this resets the default font.

If want to change the font on 1 axis or one part of a title then put one of the following symbols immediately before the part you want to change: (from text p 253)

!3 Roman font

!4 Greek font

!6 complex Roman font

!X revert to entry font

!C begin new line

!D select subscript level and character size

!U select superscript level and character size

!N select normal level and character size

There are others, see help file.


Problem One

1. create an x-y plot in interactive mode. Select variations on symbols, fonts, sizes. Overplot a second set of data.

Plot F vs C using the first the standard equation

F = 9/5(C) + 32 (1)

and then the approximation

F = 2C + 32 (2)



Some keywords to play around with:

color

thick, xthick, ythick, charthick, charsize (making all of these = 2 looks nice)

xrange, yrange

plot, plot_oo (log-log), plot_io, plot_oi (can also be achieved by /xlog, /ylog, etc.)

psym

xstyle, ystyle

linestyle

position (positions placement of plot. See page 236)

(see p. 242 for others)


Some useful tables:


1. Plot-related system variables (p. 235)

Variable Purpose

!p general plot properties

!x x-axis properties

!y y-axis properties

!z z-axis properties

!d graphics window/device properties (read-only)


You would use these when you want to change a plot characteristic for all your plots. For example, if you wanted all your plots to be dashed lines instead of solid, you would type

!p.linestyle = 2

But if only want to change it for one plot, it would be better to include the linestyle as part of the plot command, such as,

plot, indgen(100), linestyle=2

2. Symbol codes (p. 230) [psym = # ]

Symbol Code

none 0

plus 1

asterisk 2

dot 3

diamond 4

triangle 5

square 6

cross 7

user defined 8



3. Linestyles (from help) [linestyle = # ]

linetype code

solid 0

dotted 1

dashed 2

dash dot 3

dash dot dot 4

long dashes 5


4. Colors:

See text pp248ff.


Answer:

1. Type in the following sequence of commands (without the comments):

x=findgen(200)-100 ; create an array of floats between -100 and 100

y=9/5*x+32

plot, x, y ; plot x vs y, -100<=x<=100, -100<=y<=150

; y limits chosen by idl

2. to overplot the second equation, now type in

y=2*x+32

oplot, x, y

(or might be better to use a new variable such as y2, so type

y2=2*x+32 ; keep x the same

oplot, x, y2)

Now there should be 2 axes ( -100<=x<=100, -100<=y<=150) with 2 lines intersecting at (0, 32).

3. Playing with symbols, fonts, and sizes. To assign the following titles:

Title of graph = Degrees C vs F

x-axis title = Degrees Fahrenheit

y-axis title = Degrees Celsius


need to type in the following commands

title='Degrees C vs F' ; notice the quotes

xtitle='Degrees Fahrenheit'

ytitle='Degrees Celsius'


then retype the plot command with these modifications

plot, x, y, title=title, xtitle=xtitle, ytitle=ytitle

(or could put all in 1 command line this way:

plot, x, y, title='Degrees C vs F',$

xtitle='Degrees F', ytitle="Degrees C'




The text gives a neat example:

title='!3CO!D2!N Spectral Absorption Feature!X'

(CO2 Spectral Absorption Feature)

Also for Planck's function (later), notice the use of ''

title='!7Planck''s Function'

(Planck's Function)

The double '' allows the symbol, ', to be printed.


Colors:

See text pp248ff.

I have tried this series of commands

device, decomposed=0

tvlct, 0, 255, 0, 1

plot, expression, color=1

and I got a green plot.


There is also a short procedure in the book called loadcolors which might be useful, but haven't tested it yet because would need to type procedure into my directory first (and too lazy).


Also, tried the following commands

device, decomposed=1

color = 'FF0000'XL ; or any other color from table p. 251

plot, expression, color=color

and got a blue plot. Make sure to use 0 (zero) and not O (oh), or it won't work correctly.

Haven't tried making individual parts of plots different colors, yet.






2. Writing a procedure.


When it begins to be a lot of commands, or commands become somewhat complicated it is time to leave the interactive mode and write a procedure. Procedures will be created within your favorite editor. Mine is gvim which doesn't even work on yours, so let's say you would use vi.


Step 1. Create a new file with a name which describes the procedure. For example, for the Wein's Displacement Law problem, I named the procedure Planckfn. Therefore the file must be named Planckfn.pro. The name of the file must be the same as the name of the procedure with a .pro added to the end.

Step 2. The first line (not including comments) of the procedure is the word pro plus its name, e.g.,

pro Planckfn

Step 3. Type in the commands, declarations, whatever, just as you would in the interative mode. Remember to save periodically because there is nothing sadder than having to type in an entire procedure all over again when you accidentally delete it.

Step 4. Type the line

end

which tells IDL that the procedure is finished.

Step 5. Save

Step 6. Compile

Step 7. Repeat until you do not get any compiler error messages:

Read any error messages.

Correct any mistakes.

Compile.

Step 8. (Pray it does what you want) Run.


Helpful information:


Go to the online help, and type in

Program Control

You will get the links to these following categories:


IDL contains various constructs for controlling the flow of program execution, such as conditional expressions and looping mechanisms. These constructs include:


Compound Statements


· BEGIN...END


Conditional Statements


· IF...THEN...ELSE

· CASE

· SWITCH


Loop Statements


· FOR...DO

· REPEAT...UNTIL

· WHILE...DO


Jump Statements


· BREAK

· CONTINUE

· GOTO


Below I have printed out the info on the most commonly used control structures.


FOR-loops

What if you want to do exactly the same thing to every member of a data set? If that data set is in an array (most likely) then you can iterate through that array and perform a block of code on every element. And if I forgot to mention it above, the first element in a 1-dimensional array is element 0 (not 1). Likewise in a two dimensional array, the first element would be at location (0,0). Therefore, a FOR-loop will start at 0 (not 1).


Here is the information from the help file on FOR-loops


The FOR statement executes one or more statements repeatedly, incrementing or decrementing a variable with each repetition, until a condition is met.


Note - For more information on using FOR and other IDL program control statements, see Program Control.


Syntax: how a FOR-STATEMENT must be written

You must have the words FOR, DO (with or without BEGIN) and ENDFOR

You must also have a variable = some initial value, final value

It is optional to include how much too increment. The default is +1.


FOR variable = init, limit [, Increment] DO statement


or


FOR variable = init, limit [, Increment] DO BEGIN


statements


ENDFOR


Example


The following example iterates over the elements of an array, printing the value of each element:


array = ['one', 'two', 'three']

n = N_ELEMENTS(array)

FOR i=0,n-1 DO BEGIN

PRINT, array[i]

ENDFOR



WHILE-statements

You will want to use a while-statement when you want to perform the same block of code on a large chunk of data, but you don't know exactly how many data values there are, but you do know an ending condition. The most common example of this is when you want to read in data from a file until you reach the end of the file (you don't know how big the file is, but you know it ends).


The help file has:


The WHILE...DO statement performs its subject statement(s) as long as the expression evaluates to true. The subject is never executed if the condition is initially false.


Note - For information on using WHILE...DO and other IDL program control statements, see Program Control.


Syntax:

You must include the word WHILE and DO (with or without BEGIN)

expression = something which evaluates to true or false AKA a boolean expression


WHILE expression DO statement


or


WHILE expression DO BEGIN


statements


ENDWHILE


Example


i = 0

WHILE (i EQ 1) DO PRINT, i


Because the expression (which is false in this case) is evaluated before the subject statement is executed, this code yields no output.


Be careful not to create infinite loops. This will happen when your expression never changes, i.e., it is always true.


IF-Statements:

Use this when you want a command to be used conditionally.



The IF...THEN...ELSE statement conditionally executes a statement or block of statements.


Note - For information on using IF...THEN...ELSE and other IDL program control statements, see Program Control.


Syntax


IF expression THEN statement [ ELSE statement ]


or


IF expression THEN BEGIN


statements


ENDIF [ ELSE BEGIN


statements


ENDELSE ]


Example


The following example illustrates the use of the IF statement using the ELSE clause. Notice that the IF statement is ended with ENDIF, and the ELSE statement is ended with ENDELSE. Also notice that the IF statement can be used with or without the BEGIN...END block:


A = 2

B = 4

IF (A EQ 2) AND (B EQ 3) THEN BEGIN

PRINT, 'A = ', A

PRINT, 'B = ', B

ENDIF ELSE BEGIN

IF A NE 2 THEN PRINT, 'A <> 2' ELSE PRINT, 'B <> 3'

ENDELSE


IDL Prints:


B <> 3


Under IF - avoiding, I found this, which is very pertinent for programs which will be expected to manipulate large quantities of data.


Note: Programs with array expressions run faster than programs with scalars, loops, and IF statements. Some examples of slow and fast ways to achieve the same results follow.


Example¾Summing Elements


The first example adds all positive elements of array B to array A.


;Using a loop will be slow.

FOR I = 0, (N-1) DO IF B[I] GT 0 THEN A[I] = A[I] + B[I]

;Fast way: Mask out negative elements using array operations.

A = A + (B GT 0) * B

;Faster way: Add B > 0.

A = A + (B > 0)


When an IF statement appears in the middle of a loop with each element of an array in the conditional, the loop can often be eliminated by using logical array expressions.


Example¾Using Array Operators and the WHERE Function


In the example below, each element of C is set to the square-root of A if A[I] is positive; otherwise, C[I] is set to minus the square-root of the absolute value of A[I].


;Using an IF statement is slow.

FOR I=0,(N-1) DO IF A[I] LE 0 THEN $

C[I]=-SQRT(-A[I]) ELSE C[I]=SQRT(A[I])

;Using an array expression is much faster.

C = ((A GT 0) * 2-1) * SQRT(ABS(A))


The expression (A GT 0) has the value 1 if A[I] is positive and has the value 0 if A[I]is not. (A GT 0)* 2 - 1 is equal to +1 if A[I] is positive or -1 if A[I] is negative, accomplishing the desired result without resorting to loops or IF statements.


Another method is to use the WHERE function to determine the subscripts of the negative elements of A and negate the corresponding elements of the result.


;Get subscripts of negative elements.

negs = WHERE(A LT 0)

;Take root of absolute value.

C = SQRT(ABS(A))

;Negate elements in C corresponding to negative elements in A.

C[negs] = -C[negs]


Problem Two


See Wein's displacement Law handout.


Answer:

There are a few new things in here, so I have commented excessively in italics. Also, there are many, many ways to accomplish the same output. Some aspects of coding are just preference. Or knowing of a shortcut. The general rule for beginning programming is that: if it works, it is correct. Later on you can worry about whether it is doing it in the most efficient way, or if you can write it with less lines of code. Your code will become better simply by doing more and trying more - perhaps surfing the web a bit.

;I start first by creating a file in my directory called Planckfn.pro. All procedures must ;start with the word pro and end with the word end. The next word must be the name of ;the file without the suffix .pro. This mandates that the first line of my procedure be


pro Planckfn


;generate array of wavelengths from 10-7 to 10-5

l=findgen(100)

l=l+1

l=l/1E7


;To satisfy my own curiosity, I often ask the computer to show me what I have ; done just so I can see it. The following lines are not necessary, likewise almost

; all print statements in this procedure.

print, "Wavelengths:"

print, l

print, ""


;calc exponential part (exp(hc/klT) -1) for T=300, 350, 400, 450, 500, 550 K

;T = last part of variable name

exppart300=EXP(6.6261*2.998/(1.381*1E3*300*l))-1

exppart450=EXP(6.6261*2.998/(1.381*1E3*450*l))-1

exppart500=EXP(6.6261*2.998/(1.381*1E3*500*l))-1

exppart550=EXP(6.6261*2.998/(1.381*1E3*550*l))-1

exppart400=EXP(6.6261*2.998/(1.381*1E3*400*l))-1

exppart350=EXP(6.6261*2.998/(1.381*1E3*350*l))-1


; to see if results what I expect, look at data for T = 500K

print, exppart500


;calculate intensity

print, ""

print, "Intensities at 500 K"

intensity450=2*(6.6261/1E34)*2.998^2*1E16/(l^5*exppart450)

intensity350=2*(6.6261/1E34)*2.998^2*1E16/(l^5*exppart350)

intensity300=2*(6.6261/1E34)*2.998^2*1E16/(l^5*exppart300)

intensity500=2*(6.6261/1E34)*2.998^2*1E16/(l^5*exppart500)

intensity400=2*(6.6261/1E34)*2.998^2*1E16/(l^5*exppart400)

intensity550=2*(6.6261/1E34)*2.998^2*1E16/(l^5*exppart550)


; print array of intensities

print, intensity500

print, ""


; get max of intensity array

;This is done by using the WHERE , EQ , and MAX functions.

;WHERE will return the subscipt from the array where the condition in the ;following parentheses is true. First I check to see if indeed the correct value is ;chosen by the MAX function which returns the maximum value of the given ;array.

print, 'The maximum intensity is ', MAX(intensity500)


;Now get the subscripts

sub300 = WHERE(intensity300 EQ MAX(intensity300))

sub450 = WHERE(intensity450 EQ MAX(intensity450))

sub350 = WHERE(intensity350 EQ MAX(intensity350))

sub500 = WHERE(intensity500 EQ MAX(intensity500))

sub400 = WHERE(intensity400 EQ MAX(intensity400))

sub550 = WHERE(intensity550 EQ MAX(intensity550))


print, 'The max intensity is at subscript ',sub500


; get the wavelength with max intensity. Here I am assigning the variable ; maxima500 the value at sub500 of the l array.

maxima500 = l(sub500)

print,'The maxima is', maxima500


; calculate maxima - max(lambda)

; max(lambda) = 2898um-K/T

max450 = 2.898/(1E3*450)

max350 = 2.898/(1E3*350)

max300 = 2.898/(1E3*300)

max500 = 2.898/(1E3*500)

max400 = 2.898/(1E3*400)

max550 = 2.898/(1E3*550)


print, 'Lambda max is ',max500

print, 'maxima - lambda max = ', maxima500 - max500


; Now will start making the plots.

; Change font type for everything following

!p.font=1


; The following is how to make an overplot using two different y-axes.

; Set up axes with no data (using /nodata command). By choosing plot with ; largest variation, guarantee that the range will be appropriate,

; although also set ranges with x and yrange commands.

plot, l, intensity550, /nodata, ystyle=4,xstyle=1,$

title='!7Planck''s Function',$

xrange=[1/1E7, 1/1E5], xtitle='lambda (!7m)',$

thick=2, xthick=2, charthick=2, charsize=2


;left axis is yaxis=0. /save causes the scaling established by axis to be saved for

; subsequent overplotting.

axis, yaxis=0,yrange=[0, 2.3E8], /save, ystyle=1, $

ytitle='intensity ()', ythick=2, charthick=2, charsize=2


;overplot Planck's function for given T

;plot intensity as a function of lambda at all temperatures using different symbols

; for each set

oplot, l, intensity550,psym=1

oplot, l, intensity450, psym=2

oplot, l, intensity350, psym=3

oplot, l, intensity300, psym=4

oplot, l, intensity500, psym=5

oplot, l, intensity400, psym=6


;overplot lambda max array onto the Planck fn plot

lmax = [max300, max350, max400, max450, max500, max550]

T = [300, 350, 400, 450, 500, 550]


; Now the right axis (yaxis=1) is created.

axis, yaxis=1, yrange=[300, 600], /save, $

ytitle='temperature(K), Xs)', $

ythick=2, charsize=2, charthick=2

oplot, lmax, T, psym=7


; perform a linear regression on lmax vs T

; define a 6-element vector of Poisson measurement errors. I found the

; information for a linear regression in the help file.

measure_errors=SQRT(ABS(T))

result = LINFIT(lmax, T, MEASURE_ERRORS=measure_errors)


; The result printed here is two values, a and b where a is the slope of the linear ; regression and b is the y-intercept. To actually plot the line, I would need to ; generate a line using a and b which I haven't tried yet.

print, result


; Must always have the following word at the end of a procedure.

end


; Another way to do it using a for loop.


pro Planckfn2


;generate array of wavelengths

l=findgen(100)

l=l+1

l=l/1E7


; generate array of temperatures

T= [300, 350, 400, 450, 500, 550]


; create array for wavelengths at given T with max intensity

lmax = findgen(5)


;change font type

!p.font=1


; for T=550 calc exponential part (exp(hc/klT) -1)

exppart=EXP(6.6261*2.998/(1.381*1E3*550*l))-1


; calc intensity

intensity=2*(6.6261/1E34)*2.998^2*1E16/(l^5*exppart)


; get max of intensity array

sub = WHERE(intensity EQ MAX(intensity))


; get the wavelength with max intensity

maxima = l(sub)


;calc maxima - max(lambda)

; max(lambda) = 2898um-K/550

lmax(5) = 2.898/(1E3*550)

print, 'maxima - max(lambda) = ', maxima-lmax(5)


;set up axes with no data

plot, l, intensity, /nodata, ystyle=4,xstyle=1,$

title='!7Planck''s Function',$

xrange=[1/1E7, 1/1E5], xtitle='lambda (!7m)',$

thick=2, xthick=2, charthick=2, charsize=2


axis, yaxis=0, yrange=[0, 2.3E8], /save, ystyle=1, $

ytitle='intensity', ythick=2, charthick=2, charsize=2


; for each T do calculations and make an overplot

for i = 0, 5 do begin

;calc exponential part (exp(hc/klT) -1)

exppart=EXP(6.6261*2.998/(1.381*1E3*T(i)*l))-1


; calc intensity

intensity=2*(6.6261/1E34)*2.998^2*1E16/(l^5*exppart)


; get max of intensity array

sub = WHERE(intensity EQ MAX(intensity))


; get the wavelength with max intensity

maxima = l(sub)


;calc maxima - max(lambda)

; max(lambda) = 2898um-K/T

lmax[i] = 2.898/(1E3*T(i))

;print, 'maxima - max(lambda) = ', maxima-lmax(i)


;overplot Planck's function for given T

oplot, l, intensity,psym=i+1

endfor


;overplot lambda max array onto the Planck fn plot

axis, yaxis=1, yrange=[300, 600], /save, $

ytitle='temperature(K), Xs)', $

ythick=2, charsize=2, charthick=2

oplot, lmax, T, psym=7


; perform a linear regression on max vs T

; define an 6-element vector of Poisson measurement errors:

measure_errors=SQRT(ABS(T))

result = LINFIT(lmax, T, MEASURE_ERRORS=measure_errors)

print, result

end




3. Reading or Writing Files


The following is a program which opens a file called dump.txt which has an unknown number of lines of 5 numbers. We know the numbers are longitude, latitude, SO2, and we don't remember what the last two columns are. We are going to read the first three columns into 3 separate arrays and the last two columns into two dummy arrays, for now.

NOTE: There is a rule in IDL that says that you cannot read into a subscripted variable, so must first create temporary variables (which are obviously not subscripted) and reassign these to subscipted variables.

; ***********************PROGRAM OPENDATA*************************


pro openData


; I could put dump.txt wherever I use the variable name filename, but

; doing it this way makes the procedure more reusable.

filename = 'dump.txt'


; open the file dump.txt

; When files are opened, they need to be assigned a logical unit number (lun).

; Rather than make up one myself, using the following commands will let

; the computer assign the lun automatically and keep track of it for us.

; The important thing is to release the lun when the procedure is done.

OPENR, lun, filename, /GET_LUN


;read data into 5 arrays

lon = FLTARR(10000) ; longitude

lat = FLTARR(10000) ; latitude

so2Data = FLTARR(10000) ; SO2 data

garbage = FLTARR(10000) ; dummy

stuff = FLTARR(10000) ; dummy

lo = 0.0 ; The following are temporary variables

la = 0.0 ; See note above.

s = 0.0

g = 0.0

st = 0.0

count = 0 ; a variable which will keep track of

how many data rows are read


; This loop will keep reading formatted lines until it hits the end of the file

; which has been assigned the lun we are using. If we had known how many

; rows were in the file, or if that were header information, we could have used

; a for loop. Or even faster, dump the data into an one array with known

; dimensions.

WHILE (NOT EOF(lun)) DO BEGIN

READF, lun, lo, la, s, g, st

lon(count) = lo

lat(count) = la

so2Data(count) = s

garbage(count) = g

stuff(count) = st

count = count + 1

ENDWHILE

lon = lon(0:count-1)

lat = lat(0:count-1)

so2Data = so2Data(0:count-1)

garbage = garbage(0:count-1)

stuff = stuff(0:count-1)


;print, lon

; CAN PUT WHATEVER MANIPULATIONS HERE.

; SEE BELOW FOR POSSIBILITIES


; free lun to be used for another day

Free_Lun, lun

END



4. Contour, Surface, etc. plots.


A. Contour and surface are very straightforward with the commands, Contour and Surface, respectively, followed by a 2-dimensional array of data, then the two vectors the data is plotted with respect to. I think the 2-dimensional array is the data indexed by (x-value index, y-value index). In my data, this means that I create an array, but the only values I am interested in are the diagonals.

Both Contour and Surface have many keywords which can be played around with to make everything look nice.

There is also a Shade_Surf command which creates a surface which is shaded instead of gridded like Surface. Again, it has many keywords.


;contour plot of so2 data

Contour, so2Data, lon, lat, XStyle=1, YStyle=1, NLevels=10


;surface plot of so2 data

Surface, so2Data, lon, lat, Xtitle='longitude', Ytitle='latitude'


B. Maps.


;location of data on plot centered at 30 long, 0 lat, using a lambert

; projection and a scale of 1:50x106

map_set, 30, 0, /lambert, scale=50E6

map_continents ; draws continents

oplot, lon, lat, psym=1 ; overplots data with + symbol







C. Combining plots. Use the !P.Multi = [a, b, c, d, e] command. Here is an example:


;put next 2 plots onto one screen in one row

!P.Multi = [0, 2, 1, 0, 0]


a ==> The number of graphics plots remaining to plot on the display. It is normally set to 0, meaning that as there are no more graphics plots remaining to be output on the display, the next graphic command will erase the display and start with the first of the multiple graphics plots.

b ==> number of graphics columns on the page

c ==> number of graphics rows on the page

d ==> number of graphics plots stacked in the Z direction

e ==> The graphics plots are going to be displyed by filling up the rows first,

if = 0.

If = 1, then columns are to be filled up first.


D. Here is an entire program much like the one above, but I changed the way the data is read. Instead of reading the SO2 data into a 2-dimensional array, I read it into a vector. Then I used a special program, builtin IDL, which plots x-data vs y-data vs z-data. Try it.


; ***********************PROGRAM OPENDATA*************************


pro plotData3d


filename = 'dump.txt'


; open the file dump.txt

OPENR, lun, filename, /GET_LUN


; read data into 5 arrays

lon = FLTARR(10000)

lat = FLTARR(10000)

so2Data = FLTARR(10000)

garbage = FLTARR(10000)

stuff = FLTARR(10000)

lo = 0.0

la = 0.0

s = 0.0

g = 0.0

st = 0.0

count = 0

WHILE (NOT EOF(lun)) DO BEGIN

READF, lun, lo, la, s, g, st

lon(count) = lo

lat(count) = la

so2Data(count) = s

garbage(count) = g

stuff(count) = st

count = count + 1

ENDWHILE

lon = lon(0:count-1)

lat = lat(0:count-1)

so2Data = so2Data(0:count-1)

garbage = garbage(0:count-1)

stuff = stuff(0:count-1)


; put next 2 plots onto one screen

; !P.Multi = [0, 2, 1, 0, 0]


; first plot is the 3D plot of so2 data

XPLOT3D, lon, lat, so2Data


; second plot is location of data on world plot

map_set, 30, 0, /lambert, scale=50e6

map_continents

OPLOT, lon, lat, PSYM=1


; release lun value

FREE_LUN, lun


END


+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


Random data:


; create a set of data, all values random

lon = RANDOMU(-100L, 100)*360.0 - 180.0


PRINT, lon


lat = RANDOMU(-200L, 100)*180.0 - 90.0


PRINT, lat


; so2Data = lon + lat


; print, so2Data


WIDGETS


XMANAGER

Problem Three - Comparing datasets of different spatial resolutions.


This problem could be easily encountered whenever you are working with more than one set of sensor observations. In this case we are dealing with TOMS and MODIS datasets, which have spatial resolutions of 39 km and 1 km, respectively.


We want to compare the SO2 retrievals. Because of the spatial disparity, we would like to find all the MODIS pixels which occur within each TOMS pixel (roughly 1500), get the mean value and other statistical tests, and compare them directly.


There are several components to solve for this problem, and you can basically write your program to solve them sequentially:


(1) read in each of the satellite datasets. We assume that you are working with datasets that are composed of columns of latitude, longitude, and SO2 value.

(2) make a geometric calculation of a TOMS pixel boundary, given it's midpoint

(3) identify subsets based on specific conditions (e.g., all MODIS pixels contained within a TOMS pixel)

(4) perform statistical analyses of a dataset: mean, median, mode, max, min, and percent MODIS coverage of each TOMS pixel.

(5) plot the data: e.g., mean MODIS value versus TOMS


(6) once it works, you will probably find additional problems to solve or variations. Another component to consider then, is to put in "safety" features, which involve things like potential rejection criteria, or dealing with unexpected datasets (e.g., a smaller or larger set, missing data).

1At the end of this section I have included the entire help file on contour, if you are curious.