Документ взят из кэша поисковой машины. Адрес оригинального документа : http://hea-www.harvard.edu/AstroStat/SolStat2012/SolarTutorial/solarstats_aia_idl8.pro Дата изменения: Thu Feb 16 00:43:05 2012 Дата индексирования: Tue Oct 2 10:41:13 2012 Кодировка:

pro solarstats_aia_idl8

; SolarStats Workshop: 2012/02
; A sample script to process AIA data, meant primarily to be
; entered line-by-line at the command prompt for instruction
; on basic AIA processing and introdcution of useful IDL concepts.

; NOTE: As the name implies, this script is intended for IDL version
; 8.1. Most people in solar physics are currently using 7.1, as
; there a number of known issues with 8.1 and some of the SolarSoft
; routines (SSW) are not yet compatible. These issues are simple to
; resolve for our purposes here and are noted when they come up with
; **IDL 8 ISSUE** preceeding the relevant text.

; At this point you should have some data downloaded. If not,
; refer to the workshop website. You'll also need to have IDL
; installed along with the SSWIDL suite of routines. Again, refer
; to the workshop website if you're not at this point yet.

; Set dir equal to wherever the AIA data lives. Note that strings
; can be contained in either singe (') or double (") quotes.
dir = '/Users/pmccaule/Desktop/solartutorial/aia_data/'

; We'll also set up a directory to be used for output:
dir_out = '/Users/pmccaule/Desktop/solartutorial/aia_data/output/'

; We can use IDL to test if this directory exists and create it
; if it doesn't. The FILE_TEST() function returns 1 if the
; file is found and 0 if not.
if file_test(dir_out, /directory) eq 0 then file_mkdir, dir_out

; Here we create a string array of all the files in dir that
; end in '.fits', where * can be any set of characters
list = file_search(dir, '*.fits')

; Once you've created a variable, it's often useful to use the
; HELP and PRINT procedures to check that things worked as
; expected. We should have an array of 140 file paths.
help, list
print, list[0]

; Now we need to sort through our list to pull out
; observations of the same wavelength. To do this, we
; can use the WHERE() function, which returns array indices
; that match a specified citeria. The criteria we'll use is
; whether or not an element within list contains the phrase
; '171.fits'. For this we can use STRPOS(), which will search
; for this phrase and return an array with the same number of
; elements as list but populated with either (a) -1 if the
; phrase isn't found or (b) the location within the string
; where the phrase begins, with 0 indicating the first character.
search171 = strpos(list, '171.fits')

; By asking WHERE() to give us every index that does NOT
; contain -1, we can pull out just the indices we want:
indices171 = where(search171 ne -1)

; The resultant array can then be used to directly index list
; to create a new list with just the 171 observation file paths.
list171 = list[indices171]

; Or, more simply, we can wrap that all up in a single line...
list304 = list[where(strpos(list, '304.fits') ne -1)]

; To check that the arrays were populated, we can use help.
; There should be 20 elements in each array.
help, list171, list304

; Now that we have the files that we want, we can
; read in the data for manipulation. For small data
; sets (<~ 10 images) this can be done all at once. Often
; though, you'll want to process larger sets, in which
; case reading and processing is best done within a for
; loop, one image at a time, to manage memory. Here, we'll
; do one image and then later show how it might be looped
; at the very end of this script.

read_sdo, list171[0], index171, data171, /uncomp_delete
read_sdo, list304[0], index304, data304, /uncomp_delete

; Here we have two 4096x4096 data arrays and two corresponding
; 1-element structure arrays.
help, data171, index171, data304, index304

; To have a look at what info is contained in the structure...
help, index171, /str

; The next step is to 'prep' the data. This corrects for
; translational and rotational offets between the different
; telescopes and isn't necessary for every application, but
; is imporant for aligning/overlaying data from the different
; passbands.
aia_prep, index171, data171, index171p, data171p
aia_prep, index304, data304, index304p, data304p

; To see what AIA_PREP did...
print, index171.xcen, index171.ycen, index171.crota2
print, index171p.xcen, index171p.ycen, index171p.crota2

; There are many ways to plot data in IDL. One of easiest ways
; to get a quick look at the data is to use PLOT_IMAGE, which wraps
; a few of the proceeding scaling steps into one procedure.
plot_image, data171p

; Now we'll scale the data for optimum display. A good first
; step is generally to normalize for exposure time. Note the
; use of the TEMPORARY() function, which is useful for conserving
; memory when performing operations on large arrays as it avoids
; making a new copy of the variable being modified.
data171p = temporary(data171p) / index171.exptime
data304p = temporary(data304p) / index304.exptime

; Usually a logarithmic scaling is best. Note that the '>1e-6' tells IDL
; to only take the log of values within data171p that are
; positive. The calibration already applied before downloading the data
; will have set some pixels to negative values or to 0, so this is
; one way to account for that.
log171 = alog(data171p >1e-6)
log304 = alog(data304p >1e-6)

; To see the difference that the log scaling made...
plot_image, log171

; Now we'll take two steps that PLOT_IMAGE did for us. The first
; is byte scaling, which scales all the values in our array to
; a specified range.
log171 = bytscl(temporary(log171), min=1.0, max=8.0)
log304 = bytscl(temporary(log304), min=1.0, max=8.0)

; Often we want to output images that aren't quite as large as
; the full resolution AIA data. For this we can use the following,
; which will shrink our array down to 1024x1024px. Note that REBIN()
; only works for resizing images to integer multiples of the
; original dimensions. CONGRID() can be used to resize to any
; arbitrary set of dimensions.
rebin171 = rebin(log171, 1024, 1024)
rebin304 = rebin(log304, 1024, 1024)

; Now we'll set our plot window to be exactly the size of the
; image we'd like to output.
window, xsize=1024, ysize=1024

; To display, we'll use the TV procedure. This one of the
; more basic ways to plot images because it doesn't attempt to do
; any scaling. TVSCL and PLOT_IMAGE can also be used and will
; automatically scale data with generally good, but mixed, results.
tv, rebin171
tv, rebin304

; To add color we'll use AIA_LCT, which is the default way to
; colorize AIA images. This routine is adapted from the basic IDL
; routine LOADCT. The vectors rr, gg, and bb store the color
; table values, which will be needed to write color images to files.
aia_lct, rr, gg, bb, wavelnth = 171, /load

tv, rebin171

; Often you'll want to annotate the image with the wavelenth and
; observation time. To do that, use the XYOUTS procedure. Note that
; the /NORMAL keyword specifies a normalized coordinate system from
; 0.0 to 1.0. The function STRCOMPRESS here is used to remove the blank
; spaces that are introduced after converting the wavelength number
; to a string.
xyouts, 0.01, 0.01, index171p.date_obs, /normal, size=2
xyouts, 0.01, 0.03, 'SDO/AIA '+strcompress(string(index171p.wavelnth), $
/remove_all), /normal, size=2

; Now we're ready to output the image to a file. WRITE_PNG asks for a
; filename, image array, and color vectors. For the image array, here
; we can use the TVRD() function, which returns whatever is on the current
; display window; use the TRUE=1 keyword to preserve color. Analogous
; procedures exist for other file types: WRITE_GIF, WRITE_JPEG, WRITE_TIFF

; **IDL 8 ISSUE** There is a problem in 8.1 when using the X-window
; display and TVRD() to output image files. This is the default graphics
; device and what we've been using so far. As I understand it, this issue
; will be resolved in 8.2 and the 'Z' buffer can be used instead for now.
; The Z buffer is a 'pseudo-device' that 'displays' graphics only in
; memory (as opposed to on the screen) and can then be used to output
; images files.

set_plot, 'Z' ;to switch to Z buffer

; Instead of the WINDOW procedure, we use the DEVICE procedure and the
; set_r variable to specify our window size. We also need to specify the
; pixel depth as 24-bit if we want true color output (the default 8-bit
; won't allow us to silmutaneously use the different AIA color tables)
; as well as setting the DECOMPOSED keyword to 0. Refer to the help files
; for additional information. Getting color to output properly can be
; tricky in IDL but there is a lot of information out there on it.
device, set_r=[1024, 1024], set_pixel_depth=24, decomposed=0
tv, rebin171
xyouts, 0.01, 0.01, index171p.date_obs, /normal, size=2
xyouts, 0.01, 0.03, 'SDO/AIA '+strcompress(string(index171p.wavelnth), $
/remove_all), /normal, size=2

write_png, dir_out+'solstats171.png', tvrd(true=1), rr, gg, bb

; Now for the 304 image:
aia_lct, rr, gg, bb, wavelnth = 304, /load
tv, rebin304
xyouts, 0.01, 0.01, index304p.date_obs, /normal, size=2
xyouts, 0.01, 0.03, 'SDO/AIA '+strcompress(string(index304p.wavelnth), $
/remove_all), /normal, size=2
write_png, dir_out+'solstats304.png', tvrd(true=1), rr, gg, bb

set_plot, 'x' ;to restore x-window graphics display

; Now we'll consider plotting both images in the same plot window.
; To do this, we change the !P.MULTI system variable. (Note that IDL
; system variables all begin with '!'.) !P.MULTI is a 5-element array
; that begins as all 0's. To change the number of plot columns, modify
; !P.MULTI[1] and !P.MULTI[0] for rows. Refer to help files for more info.
window, xsize=1024, ysize=512
!p.multi = [0,2,0,0,0]
aia_lct, rr, gg, bb, wavelnth = 171, /load

; PLOT_IMAGE is great to use for multiple plots in the same window
; because it does a good job resizing for you. Note that here we can use
; the full size image array, as opposed to the binned down version, because
; of the automatic resizing. Also notice the SCALE keyword, which tells
; PLOT_IMAGE how many units per pixel there are in our image. If we had used
; the 1024px binned down version, this number would have been 2.4.
plot_image, log171, title='AIA 171 '+STRING(197B), xtitle=index171p.date_obs, $
ytitle='arcsec', scale=0.6
aia_lct, rr, gg, bb, wavelnth = 304, /load
plot_image, log304, title='AIA 304 '+STRING(197B), xtitle=index304p.date_obs, $
ytitle='arcsec', scale=0.6

;**IDL 8 ISSUE** Same as previous
set_plot, 'z'
device, set_r=[1024,512], set_pixel_depth=24, decomposed=0
aia_lct, rr, gg, bb, wavelnth = 171, /load
plot_image, log171, title='AIA 171 '+STRING(197B), xtitle=index171p.date_obs, $
ytitle='arcsec', scale=0.6
aia_lct, rr, gg, bb, wavelnth = 304, /load
plot_image, log304, title='AIA 304 '+STRING(197B), xtitle=index304p.date_obs, $
ytitle='arcsec', scale=0.6
write_png, dir_out+'solstats_combo.png', tvrd(true=1), rr, gg, bb

; Now if you're only interested in a specific region, we can index the array to
; only display those pixels. The syntax works like: Array[x1:x2,y1:y2]. This also
; illustrates one of IDL's quirks, which is to use the opposite of standard
; matrix notation, allowing for intuitive x/y mapping of images and the like.
aia_lct, rr, gg, bb, wavelnth = 171, /load
plot_image, log171[1624:2647,1124:2147], title='AIA 171 '+STRING(197B), $
xtitle=index171p.date_obs, ytitle='arcsec', scale=0.6
aia_lct, rr, gg, bb, wavelnth = 304, /load
plot_image, log304[1624:2647,1124:2147], title='AIA 304 '+STRING(197B), $
xtitle=index304p.date_obs, ytitle='arcsec', scale=0.6
write_png, dir_out+'solstats_combo_cutout.png', tvrd(true=1), rr, gg, bb

set_plot, 'x' ;to restore x-window graphics display

; Next we'll consider one of most common quantitative means of analysis,
; the light curve. This is as simple as adding up all the pixels in
; some box and plotting that over time. Even more simply, AIA structure indexes
; already come with a keyword containing the average pixel value in the
; corresponding data array (index.DATAMEAN). So, using the /NODATA keyword,
; we can quickly read in just the indexes and create our curve:
read_sdo, list171, in171, /nodata
read_sdo, list304, in304, /nodata

; To make a more interesting plot, lets put in all of the different EUV bands.
; Notice that I can wrap up several previous steps at once here.
read_sdo, list[where(strpos(list, '94.fits') ne -1)], in94, /nodata
read_sdo, list[where(strpos(list, '131.fits') ne -1)], in131, /nodata
read_sdo, list[where(strpos(list, '193.fits') ne -1)], in193, /nodata
read_sdo, list[where(strpos(list, '211.fits') ne -1)], in211, /nodata
read_sdo, list[where(strpos(list, '335.fits') ne -1)], in335, /nodata

!p.multi[1] = 0 ;set back to one plot per window

; We can use the procedure UTPLOT to easily create a time series plot and
; OUTPLOT to overlay different series. These two procedures are derived from
; the native IDL procedures PLOT and OPLOT and inherit all the keyword
; functionality that those two procedures have, though inherited keywords
; tend not to be listed in the header documents.
utplot, in171.date_obs, in171.datamean
outplot, in304.date_obs, in304.datamean
outplot, in94.date_obs, in94.datamean
outplot, in131.date_obs, in131.datamean
outplot, in193.date_obs, in193.datamean
outplot, in211.date_obs, in211.datamean
outplot, in335.date_obs, in335.datamean

; Now we'll introduce a number of features to increase the awesomeness of our
; plot. The first is to set up our y-axis on a log-scale since there is such
; a large range of values. To do this we use the /YLOG keyword. Using the
; YRANGE keyword, we can specify the range over which to plot. IDL treats this
; range as more of a suggestion unless you also set the YSTYLE keyword equal
; to 1, which forces IDL to give you the range you asked for.
utplot, in171.date_obs, in171.datamean, ystyle=1, /ylog, yrange=[1,500]
outplot, in304.date_obs, in304.datamean
outplot, in94.date_obs, in94.datamean
outplot, in131.date_obs, in131.datamean
outplot, in193.date_obs, in193.datamean
outplot, in211.date_obs, in211.datamean
outplot, in335.date_obs, in335.datamean

; Now lets add some color and a legend. Right now we have one of the AIA
; color tables loaded, which means we can only use the range of colors
; contained in that table. It'd be nicer to have a variety of colors
; to work with and one way to do that is to create our own color table
; (or to use LOADCT). First, I'll create a structure with all the colors
; I want. This way we can access the color table with something like
; ct.cyan instead of remember that cyan is the 7th color in the table.
ct = {black : 0, $
white : 1, $
red : 2, $
green : 3, $
blue : 4, $
yellow : 5, $
magenta : 6, $
cyan : 7, $
orange : 8}

; Next we create red/green/blue vectors with value cominations that
; give us the colors we want
; 0 1 2 3 4 5 6 7 8
rr = [ 0,255, 255, 0, 0,255,255, 0,255]
gg = [ 0,255, 0,255, 0,255, 0,255,125]
bb = [ 0,255, 0, 0,255, 0,255,255, 0]

; Lastly, we load the new color table
tvlct, rr, gg, bb

; Color and titles can be added to the plot using the COLOR, YTITLE,
; and TITLE keywords. Note that the color values pertain only to the new
; color table we just loaded. A different table would produce different colors.
utplot, in171.date_obs, in171.datamean, ystyle=1, /ylog, yrange=[1,500], $
color=ct.yellow, ytitle='DN', title='AIA Light Curves for X2.2 Flare on 2011/02/15'
outplot, in304.date_obs, in304.datamean, color=ct.red
outplot, in94.date_obs, in94.datamean, color=ct.green
outplot, in131.date_obs, in131.datamean, color=ct.cyan
outplot, in193.date_obs, in193.datamean, color=ct.orange
outplot, in211.date_obs, in211.datamean, color=ct.magenta
outplot, in335.date_obs, in335.datamean, color=ct.blue

; **IDL 8 ISSUE** Below we add a legend to our plot using the Astrolib
; LEGEND procedure that should have come bundled with your installation
; of SolarSoft. It seems that IDL 8 now includes its own native LEGEND()
; function, so a namespace conflict has been introduced. We could simply
; use the built-in function, but this provides a useful diversion into
; how IDL finds the routines you ask for and what to do if you run into
; such a problem.

; The directories holding all the routines available to you are contained
; in the environment variable !PATH.
print, !PATH

; When IDL is asked to bring up 'legend', it finds the first instance of
; 'legend.pro' in the directories contained within the !PATH variable. To
; check if there are multiple, we can use the WHICH procedure.
which, 'legend', /all

; Whichever version is listed first is the one IDL grabs. If we want
; IDL to be grabbing a different version, then we can append that path
; to the BEGINNING of our !PATH variable so IDL finds that one first.
; I want the legend contained within /ssw/gen/idl_libs/astron/plot/, so if
; it's not listed first I add it to the beginning. Note that this file
; structure may be a bit different on your machine.
!PATH = expand_path('+/usr/local/ssw/gen/idl_libs/astron/plot/') + ':' + !PATH

; Now we should see the right 'legend.pro' at the top of the list and
; we're good to go:
which, 'legend', /all

; And here we add the legend:
legend, ['193', '171', '211', '304', '131', '335', '094'], $
colors=[ct.orange, ct.yellow, ct.magenta, ct.red, ct.cyan, ct.blue, ct.green], $
linestyle=[0,0,0,0,0,0,0], /normal, position=[.86,.65], charsize=1.2

; That's it! After the end statement you'll find a sample setup for looping
; some of the processing steps introduced here and some suggestions for additional
; excercises using the data you've downloaded.
end

; Although procedures like read_sdo and aia_prep allow for input of image
; arrays, rather than one image at a time, it is perhaps more often than not
; best to read and process images one at a time within a loop. Due to the
; large size of AIA data, attempting to read in any more than say 10 images at
; once can be overly taxing on computational resources. What follows is the
; syntax of a for loop that will read, process, and output one image at a time
; using steps that are explained above.

;window, xsize=1024, ysize=1024
;for i=0, n_element(list)-1 do begin
; read_sdo, list[i], index_temp, data_temp, /uncomp_delete
; aia_prep, index_temp, data_temp, index_prep, data_prep
; data_prep = temporary(data_prep) / index_prep.exptime
; data_log = bytscl(alog10(temporary(log171) >1e-6), min=1.0, max=8.0)
; aia_lct, rr, gg, bb, wavelnth=index_prep.wavelnth, /load
; tv, rebin(data_log, 1024, 1024)
; xyouts, 0.01, 0.01, index171p.date_obs, /normal, size=2
; xyouts, 0.01, 0.03, 'SDO/AIA '+strcompress(string(index171p.wavelnth), $
; /remove_all), /normal, size=2
; fname = 'fname_'+trim(i,'(i3.3)')+'.png'
; write_png, dir_out+fname, tvrd(true=1), rr, gg, bb
;endfor

; A useful excercise might be to use the loop code above and the information
; provided in this sample script to create a series of images showing each EUV
; channel in the same window for each of the 20 times in this data set. The same
; thing could also be done for a cutout region. Since the sun is rotating, this
; would require moving your cutout box by some offset for each new image. (The
; BLINK() routine would be useful for this.) A light curve like the one we
; already created could then be generated for just the cutout region by summing
; the pixels within the range.