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.