IDL Tutorial for Astronomy 361 (Optical)
Prof. J. D. Monnier
This tutorial is meant to take approximately 2-3 hours and expose you to the basic functionality of IDL. After this tutorial (and with frequent consultation of idl help '?'), you should be able to effectively use IDL to analyze and present the data taken during subsequent experiments. There is no separate Worksheet to fill out, but rather the tutorial instructs you to create a 'JOURNAL' file of your IDL session and also asks you to create postscript plots at a number of junctures. Please name these files as instructed so I can review them later as part of the grading.
I have included links to other online IDL tutorials if you want to learn more; also, there are lessons in the IDL reference books in the 8th floor Dennison Computer Room.
For those finding this tutorial on the Web, you can try to follow along by installing the IDL Astronomy Library and downloading this tar file: FILES
|?||Invoke IDL Online Help|
|<up arrow>||Scroll to previous commands|
|&||Separate multiple commands on same line|
|$||Used to continue a command on the next line|
|Printing||set_plot,'ps' & <plot commands> & device,/close & set_plot,'x' ; [or set_plot,'win']|
|.run filename||Execute commands in the filename (run script!)|
|.continue||(or just '.c') allows a program/script to continue after a 'STOP' or break|
|save,/all,file='savesession.sav'||Save your IDL session to a file. Can start again next login using restore.|
|restore,'savession.sav'||Restores all files and compiled routines from previous session|
I. IDL environment (idlde)
II. Command line exercises
III. Reading in data
IV. Displaying 1-D data
V. Displaying 2-D data (images)
I. IDL environment
After logging in as yourself on the Angell Hall Mac cluster, you should copy the example startup file (~/../as361/Idlstuff/idl_startup.pro) into your own space (~/idlwork/). To start IDL, you should open an X11 terminal and type:
[carina:~] monnier% idlde
When you begin idlde, you will see the idl development environment (below picture for Windows):
The idlde environment allows you to look at local variables, edit procedures/scripts, inspect previous commands, and issue new commands. You will need to edit this file to customize it. Please open your personal copy of the idl_startup.pro using the File-> Open tool inside idlde. You should be able to edit this file in the top-right window using familiar text-editing methods. Please find the line where home_dir is dfined. Please make this say home_dir ='~/idlwork/' Once you are done. Please save and continue..
To setup the environment properly, you need to choose a startup file. To do this, use the FILE pulldown menu, choose PREFERENCES, and click on the STARTUP tab. Now find the idl_startup.pro file you just modified and select this as your startup file. Next time IDL is run, this file will be loaded and will set up your path and a few other things (I recommend you restart idlde and then take a look at the startup file). For instance, from the startup file you can see the location of the important libraries, such as the Astronomy library in as361/Idlstuff/Libraries/Astrolib/pro).
Now, moving on to some basic IDL programming. Unlike in more structured programming languages (e.g., C++, Fortran), one does not have to explicitly define variables in IDL; rather, variables are defined implicitly through their use. For example,
A FLOAT = 1.00000
I INT = 1
You see that 'a' was declared as a float, while 'i' was declared as an integer. This information is also shown in the 'Variable Watch Window.' For more information on IDL data types, you can invoke IDL online help by typing a '?': try a search for 'data types.'
IDL can also save all your interactive session by using the JOURNAL command, and this will part of your tutorial today instead of filling out a separate worksheet. Start the recording process by typing the following command.
IDL> journal,'YourName.Date.journal' ; e.g., journal, 'monnier.2005Mar23.journal'
When you exit IDL (or type 'journal' again with no filename), the recording session will end. Please use JOURNAL to record your sessions for future reference. After typing a few commands, check that the file has been created in the proper place (~/idlwork/). You may need to copy this file into as361/Deposit/ for grading -- ask your professor.
II. Command line exercises
One can create arrays, 1-D or 2-D, easily as well.
IDL> y=fltarr(50) ; semicolon marks a comment line. Use IDL> ? to find fltarr and see what this line does.
Try to understand what the following does using online help.
IDL> x= findgen(100)
IDL> plot,x,y,title='First plot',xtitle='X axis',ytitle='y axis label' ; IDL uses keywords to pass (often optional) information. Here we use the optional keyword XTITLE to label the x axis -- see online help for complete list of keywords
Note that a window (IDL 0) was created. You can create multiple plotting windows using IDL, and move between them easily. See online help for details on how to modify their default sizes, title bars, positions, etc.
IDL> window,0 & plot,x,y ; ampserand used to separate commands on the same line
IDL> window,1 & plot,x,y*y ; create new window and make a plot
IDL> wset, 0 & plot,x,x^2 ; go back to the first window and make a new plot
There are many ways to make a plot in IDL. Simple plots can be interactively created using LIVE_PLOT, or some version of LIVE_*. More generally, one makes a plot by opening a postscript device, issuing the plot commands, then closing the device -- thereby, creating a nice postscript file. See below:
IDL> set_plot,'ps' ;set postscript device
IDL> device, /landscape, filename = 'yourname.plot1.ps' ;e.g. monnier.plot1.ps
'Your Title Here'
IDL> device,/close ; close the postscript device... Did it create a good file?? Check before continuing using unix application ghostview gv (e.g.,type IDL> $ gv filename & ) and don't worry if its upside down).
IDL> set_plot,'x' ; This is important -- otherwise your next graphing commands will continue to be sent to the PS (postscript) device. "X" is short for X-windows, the graphical interface used for idl on mac and lnius. (for Windows, you would use use set_plot,'win').
What about Images? IDL is probably best used for doing image analysis. A quick way to create a dummy image in IDL is using the dist() function (IDL> ? dist):
and there are many built-in ways to display this data (use online help to see what each of these commands do):
IDL> contour,im,findgen(100)-50,findgen(100)-50,/c_ann,nlevel=20,xstyle=1,ystyle=1 ; contour has lots of options...
IDL> image_cont,im,/aspect ; using a / in front of a keyword is equivalent to setting it equal to 1, a useful shortcut for 'turning on' or enabling extra features in procedures/functions
IDL> show3, im
IDL> live_surface,im ;This program allows you to even rotate the surface around..
You can choose many color tables when viewing images. You can choose a palatte using 'xloadct' or:
IDL> loadct,3 ; you can see a list by just typing loadct w/o a number
Now replot images to see the new color table take effect:
Try other color tables (I like 3 and 10 for general use) (if color does not work for you, you might try typing: IDL> device,decompose=0)
You can plot slices through the array by selecting certain rows/columns of the image:
IDL> plot,im[50,*] ; Note that  brackets are prefered syntax for addressing array elements, while () paranthesis are meant for function calls
IDL> oplot,im[45,*],linestyle=2 ; the * wildcard lets you plot an entire row or column of an array
And it is easy to add a simple legend using the LEGEND command from the Astrolibrary:
IDL> legend,['row 50','row 45'],line=[0,2]
Since LEGEND is not a built-in IDL routine, you can not find help on it using the online help. However, you can use the following command to print out the routine's header.
You can access subarrays (a range of indices) easily using a colon
Study the system variable !p.multi using the online help. This variable is used to make multiple plots on one page. Please use !p.multi to create one page of plots 3 across and 2 up/down. You should include plot, contour, and legend commands, and please label the axes.
Save the multipanel plot as a postscript file: 'YourName.plot2.ps'
You can return to one plot per page by typing
III. Reading in image data
One of the most important functions of a data analysis package is reading in data in variety of formats. IDL has some data formats built-in or programmed in the Astronomy library such as FITS. For instance, one can read in data from the Hubble Space Telescope using available libraries. Here, you can look at short exposure image of a binary star taken with the Keck Telescope:
IDL> image = readfits(dialog_pickfile()) ; Choose the file 'phiand.fits' in the '~/../as361/Idlstuff/AY361/' directory
Now, look at this image...
IDL> shade_surf,image ;loadct,3 produces a nice color table for shade_surf
Try the following and see if you determine what the result means.
You can see other data formats supported from the online help ('? read'), including GIF, JPEG, TIFF,. PNG, etc. etc. [Note that newer versions of IDL do not continue full support of GIF because of a licensing issues -- one can use the open format PNG.]
The smooth() command is helpful for visualizing images in some circumstances. Experiment with this function. What does the EDGE keyword do?
Compare the smoothed and unsmoothed images by creating a plot (put them on the same page using !p.multi). Note that tv/tvscl do not work easily with postscript device, and you could use imcont, contour, image_cont, shade_surf, etc. Please save a copy of this plot as 'YourName.plot3.ps'.
IV. Displaying 1-D Data
In this section, we will use variables X & Y to hold columns 1 and 2 respectively of data1.txt. Here we will just explore some of the features of 1-D plotting, including curve fitting, linestyles, plot symbols, legends, and error bars.
First, load some 1-d data using the script custom_data1.pro (you can load this file into the IDL editor if you are curious to see the code -- it is also in ~/../as361/Idlstuff/AY361/). This procedure takes 4 arguments: the first one is the filename (INPUT), while the last three (header,x,y) are output variables.
IDL> filename = dialog_pickfile() ; choose data1.txt in ~monnier/AY361/General_Files
IDL> plot,x,y,xtitle="X",ytitle="Y",title="Test Plot", psym=4
Experiment with different plot symbols and linestyles-- see online help for the PLOT command.
You can easily overplot data onto an existing plot
IDL> oplot,x,y + 0.15, psym=2 ;use a different plot symbol or linestyle to tell the two datasets apart
The Astronomy Library has extra plot symbols you can use.
IDL> plotsym ; this will show you a list
IDL> plotsym,0,/fill ; this will set up plotsymbol 8 as a filled circle
IDL> plot,x,y,xtitle="X",ytitle="Y",title="Test Plot", psym=8
Usually, one has error bars associated with data. The Astronomy Library has a nice routine called ploterror (and oploterror) which allows you to plot error bars. See the following example:
IDL> err = replicate(.05,100) ; Create an array of error bars with value 0.05
IDL> ploterror, x,y,err,psym=8, xrange = [-10,10] ; XRANGE plots only a limited range along the x axis (similar for YRANGE)
IDL> plotsym,3,/fill ; setup a new plot (user defined) symbol
Remember, that you can get the help on some IDL routines using
IDL> doc_library,'ploterror' ; note that ploterror take the normal labeling keywords (xtitle,title, etc) as PLOT
For the rest of this part we will work with a reduced set of the data, and this will also introduce you to the useful IDL function called 'WHERE.'
IDL> index_array = where(x ge -10 and x le 10)
This will return the indices of the x-vector which satisfy the conditions inside the parentheses. 'ge' means 'greater than or equal to', 'le' means 'less than or equal to.' See online help for other conditional expressions ('? conditional'), such 'eq' = equals.
IDL> x0 = x[index_array]
IDL> y0 = y[index_array]
Inspect these variable arrays and plot the results. We now are only dealing with the part of the dataset near the origin. Lets try using a few of the built-in fitting routines in IDL.
IDL> result1= gaussfit(x0,y0,gauss_params)
Read the online help enough to understand what the parameters mean -- note how they have pretty different calls.
Now lets plot all of this and compare the fits:
IDL> plot,x0,y0,psym=4,title='Compare Fits',xtitle="X Axis",ytitle="Yaxis"
Now overplot (using OPLOT) result1 and yfit using different linestyles. Use the LEGEND procedure to label each curve. You see that this is quite a respectable plot, with labeled axes, a title, legend, etc. Please save a copy of this plot as 'YourName.plot4.ps'.
V. Displaying 2-D Data
Analyzing 2-dimensional data is not that different in principle from analyzing 1-D data. Different routines are used to show contour plots and intensity images. When we take CCD frames from the Angell Hall Observatory or analyze GRB fields, we will want to use IDL to graph and analyze the image. You saw some commands in the top section (surface, shade_surf, tvscl, etc); we work with these a little more here.
First let us load in an image of Andromeda taken with the Angell Hall CCD camera:
IDL> im=readfits('../../as361/Idlstuff/AY361/andromeda.fits',header1) ; or find it using restore,dialog_pickfile() if your path is not right.
We can get a feeling for this image by using help. The dimensions and datatype of 'im' can be found easily ..
IM UINT = Array[510, 340]
These data are stored as Unsigned Integers. Since we may want to perform floating operations (e.g., division), lets convert to FLOAT:
IM FLOAT = Array[510, 340]
You can also use the hprint command to look at the header fields in a FITS file.
Now we can display this image pixel for pixel onto the screen.
IDL> tvscl,im ; remember to use loadct,3 and device,decompose=0 to see colors
You can also inspect using image_cont or imcont routines.
Unfortunately, 'image_cont' is not a very flexible routine making axis-labelling problematic. One often has to write scripts to make publication quality plots for professional journals. I've included here a customized version called "imcont.pro" which has more options.
IDL> imcont,im,title="Andromeda at Angell Hall",xtitle="Pixels (x-axis)",ytitle="Pixels (y-axis)",/aspect,/nocontour
For this observation, data original CCD pixels were binned 3x3, corresponding to 1.71 arcseconds. We can relabel the the axes in terms
of arseconds by creating x and y-axis arrays as we did above for the contour commands. Using help (see above), one can see the image dimensions and
create corresponding vectors containing label information.
IDL> xvector = findgen(510) * 1.71 ; will convert "pixel number" to "relative arcseconds"
IDL> yvector = findgen(340) *1.71
IDL> imcont,im,title="Andromeda at Angell Hall",xv=xvector,yv=yvector, $
Now make a postscript version of this plot, but label axes in terms of arcminutes instead. Save as a single plot in the file 'YourName.plot5.ps'
You were introduced to a number of routines above during the command-line exercises. Try making a contour plot of 'im' array, and label the axes with the values found in the arrays 'xvector','yvector'. Other routines to practice with: shade_surf, others ?
Before continuing, make sure to view the histogram equalized version of this frame.
You can clearly see the extended structure of Andromeda, but also note the many foreground (faint) foreground stars also in the frame.
Let us concentrate on one of these by creating a subarray:
Look at this image -- note there is a faint star in the center. IDL also allows 2-dimensional fits to data using built-in fitting function.
IDL> imfit=gauss2dfit(imsub,params) ;see online help for details of the call
You can compare the image with the fit:
IDL> imcont,imsub,title="Field Star Near Andromeda",/nocont,/asp
IDL> contour,imsub,/overplot, nlev=10
IDL> contour,imfit,/overplot, nlev=10,c_linestyle=2,thick=3
IDL> legend,["Data","Gaussian Fit"],line=[0,2]
You can see the fit is good. By the way, what were the parameters of the fit?
Given the high noise level in this image, you might want to repeat the above comparison but first smoothing the image.
IDL> imcont,smooth(imsub,2),/asp,title="Field Star (smoothed)"
IDL> contour,imfit,/overplot, nlev=10,c_linestyle=2,thick=3
Another way to compare data is to make multiple plots on the screen using !p.multi. For instance, to make a 3x1 array of plots:
IDL> imcont,imsub,/asp,title="Field Star"
IDL> legend,['Gaussian Fit']
IDL> imcont,imsub-imfit,/asp ; IDL adds/subtract arrays as easy as if scalars, if arrays are compatible
IDL> legend,['Residual Image']
The exercise (and postscript plot) of this section will be to draw 1-D scans through this image, one along the peak in x and one along the peak in y (curves should be labeled in the legend). Save as a single plot in the file 'YourName.plot6.ps' For extra credit, also overplot Gaussian fits to each scan and include the fitted parameters in the legend.
Hint: You should be able to use contour and print statements to determine which pixel is the maximum. You could also experiment with the CURSOR command. Once you have found the peak pixel (x0,y0), you can plot the scans as shown in section II above.
Scripts allow you to automatically run a series of commands without having to type them in at the command line. This is particularly useful when executing a repetitive task, such as analyzing a number of data sets with the same procedure. Also, scripts are useful for making complicated plots since one generally wants to see the figure in the window before making a postscript figure. Here, we walk through a simple example of this.
For this exercise, you will be creating a file ("script") to run. As for the figures, you should name your scripts in a unique way to help us figure out which files belong to whom when grading. All your scripts should have names which obey the following style "YourName.script#.pro" so there are no conflicts (e.g., 'monnier.script1.pro'). While you work on your files in the ~/idlwork directory, you should copy them over to the as361/Deposit directory when you want them graded.
So, open a new file using ctrl-N or using pulldown menu (File->New->Editor). Now you can type commands you would on the command line and these lines will be executed when the script is run. The only different rule is that this file must end with an 'end' statement.
As an example, lets put some of the commands from the previous exercises in a script:
in your new file:
; Test script written by -name- on -date- to do -what
; (header comments can be useful for future reference)
plot, x, y, xtitle="X Axis", ytitle="Y Axis",title="Plot made on:"+systime()
Now save this file as 'YourName.script1.pro' and then compile/run the program (remember to choose the datafile "data1.txt"). The easiest way to do this with the idlde gui is to right click in the editor window and choose 'Compile...' followed by 'Run ...'. After doing this, look at the plot that it made.
What about printing this? Using scripts, this is easy:
IDL> device,/landscape, filename = 'YourName.plot7.ps'
IDL> .r filename.script ; or if you don't know the full path, etc., use the GUI to run script
IDL> set_plot,'x' ; New plots will go to the screen as usual (with linux/solaris, type: set_plot,'x' instead)
You can now print your file using a postscript viewer like ghostview (gv). You can call gv from idl, like so:
[UNIX only:] IDL> spawn,'gv filename.ps'
On a unix machine, one can also easily print this by spawning an lpr command:
IDL> spawn, 'lpr filename.ps' ; this will execute a unix shell command [not for windows]