PRO iris_burst_check, filename, threshold=threshold, ixrange=ixrange, $ no_median=no_median, output=output, quiet=quiet ;+ ; NAME: ; IRIS_BURST_CHECK ; ; PURPOSE: ; This routine aims to quickly identify bursts in a data-set. I ; compared the results against my method for the 22-Oct-2013 ; data-set (where I used a Gaussian fit) and most of the events ; are found, although some are missing. Only one event was found ; with this method and not with the Gaussian method. ; ; The IRIS calibration has decreased significantly with time, and ; so the threshold for finding bursts is scaled based on the ; 1402.77 response on the input file date relative to the ; reference date of 22-Oct-2013 21:00 UT. ; ; INPUTS: ; Filename: Name of an IRIS spectrum file. Can be an array of ; filenames. ; ; OPTIONAL INPUTS: ; Threshold: The threshold intensity for flagging burst ; pixels. Default is 500 DN/s. ; Ixrange: Used to specify a sub-range in X to search. This is ; intended for very long sit-and-stare sequences for ; which you only want to search a part of the ; sequence. The range is specified in pixel numbers. ; ; KEYWORD PARAMETERS: ; NO_MEDIAN: The routine checks ; QUIET: If set, then no text is printed to the IDL ; screen. Use OUTPUT= to send the results to a ; structure. ; ; OPTIONAL OUTPUTS: ; OUTPUT: A structure containing the results. The tags are: ; .time Time for pixel in CCSDS format. ; .file File containing pixel. ; .xpix Index of X-pixel. ; .ypix Index of X-pixel. ; .x X-position in heliocentric coords (arcsec). ; .y Y-position in heliocentric coords (arcsec). ; .thresh Intensity threshold (DN/s). ; .ybin No. of Y-pixels that are binned. ; .median Median intensity for pixel (DN/s). ; .intensity Intensity of pixel (DN/s). ; .group Index of pixel group to which event belongs. ; .file_index Index of file sequence. ; ; OUTPUTS: ; The list of burst pixels is printed to the IDL screen. ; ; PROGRAMMING NOTES: ; For each spatial pixel, the routine averages the intensity bins ; between -50 and +50 km/s from the line center of Si IV ; 1402.77. The average intensity is required to be greater than ; the THRESHOLD value in order for the spatial pixel to be ; flagged. ; ; To prevent very strong particle hits from being flagged, the ; routine also takes the median intensity value in the -50 to +50 ; km/s range and requires to the average intensity to be less than ; a factor MED_FACTOR larger than the median value. To prevent ; this check, use the /NO_MEDIAN keyword. ; ; The routine checks for Y-binning, and increases the threshold ; value accordingly. ; ; MODIFICATION HISTORY: ; Ver.1, 11-Jul-2016, Peter Young ; Ver.2, 27-Jul-2016, Peter Young ; I now check for binning in the Y-direction. ; Ver.3, 15-Dec-2016, Peter Young ; Now prints out time (useful for sit-and-stare data); Y-binning ; was not correctly treated, so I've fixed this now; I ; now check to make sure the intensity is not hugely different ; from the median intensity (sign of a particle hit); to switch ; this off use the /no_median keyword. ; Ver.4, 21-Dec-2016, Peter Young ; Expanded header; added OUTPUT optional output. ; Ver.5, 8-Mar-2017, Peter Young ; Now adjusts threshold based on the change in IRIS response ; between the input file date and my reference file date. ; Ver.6, 14-Mar-2017, Peter Young ; Now assign a group to each pixel. The group is determined by ; neighboring pixels. Updated header, expanded output structure, ; and modified print-out to screen. ; Ver.7, 20-Apr-2017, Peter Young ; Added X and Y position to output structure. ;- IF n_params() EQ 0 THEN BEGIN print,'Use: IDL> iris_burst_check, filename [, output=, threshold=, /quiet, ixrange=, /no_median]' return ENDIF ; ; Delete output structure if it was previously defined. ; IF n_tags(output) NE 0 THEN junk=temporary(output) n=n_elements(filename) lref=1402.77 ; ; The average intensity in a line profile (at a single Y-pixel) is ; computed over the velocity range -vmax to +vmax from the rest ; wavelength. The intensity needs to be greater than THRESHOLD DN/s. ; vmax=50.0 IF n_elements(threshold) EQ 0 THEN BEGIN threshold=500.0 ; ; Adjust threshold based on observation date ; d=iris_obj(filename[0]) date=d->getinfo('DATE_OBS') obj_destroy,d iresp=iris_get_response(date,version=version) getmin=min(abs(iresp.lambda-1402.77/10.),imin) area_sg=iresp.area_sg[imin,0] ; iresp0=iris_get_response('22-oct-2013 21:00',version=version) getmin=min(abs(iresp0.lambda-1402.77/10.),imin) area_sg0=iresp0.area_sg[imin,0] threshold=threshold*area_sg/area_sg0 ENDIF ; ; If /no_median set, then set med_factor very high. ; IF keyword_set(no_median) THEN med_factor=1e10 ELSE med_factor=10.0 str={time: '', file: '', xpix: 0, ypix: 0, thresh: 0., ybin: 0, median: 0., $ intensity: 0., group: 0, file_index: 0, x: 0., y: 0.} ; ; Comment on ybin ; I assume that when two pixels are binned in the Y-direction the DN ; values are summed. Therefore to implement the threshold correctly, ; I need to double the threshold (in the case of x2 binning). ; Comment on lambda-binning ; Because I compute the line width over a fixed velocity width, then ; I don't need to worry about lambda-binning. ; group_num=1 FOR i=0,n-1 DO BEGIN wd=iris_getwindata(filename[i],lref,/normalize,/keep_sat,ixrange=ixrange) ; ybin=round(wd.scale[1]/0.166) ; v=lamb2v(wd.wvl-lref,lref) k=where(v GE -vmax AND v LE vmax,nk) int=average(wd.int[k,*,*],1,missing=wd.missing) med_int=median(wd.int[k,*,*],dimension=1) k=where(int GE threshold*float(ybin) AND int LT med_factor*med_int,nk) ; IF nk NE 0 THEN BEGIN ; ; The following code groups neighboring pixels into groups, with ; each group assigned an index (group). ; I have to define imgx as region_grow does not work on edge ; pixels. ; s=size(int,/dim) img=bytarr(s[0],s[1]) img[k]=1b imgx=bytarr(s[0]+2,s[1]+2) imgx[1:-2,1:-2]=img group=bytarr(s[0]+2,s[1]+2) ; WHILE total(imgx) NE 0 DO BEGIN j=where(imgx EQ 1b,nj) IF nj EQ 1 THEN BEGIN roi_pix=j[0] ENDIF ELSE BEGIN roi_pix=region_grow(imgx,j[0],/all_neigh) ENDELSE group[roi_pix]=group_num group_num=group_num+1b imgx[roi_pix]=0b ENDWHILE group=group[1:-2,1:-2] ; ; Now go through each pixel and store information in the output ; structure. ; FOR j=0,nk-1 DO BEGIN ij=array_indices(int,k[j]) tstr=anytim2utc(/ccsds,wd.time_ccsds[ij[0]],/time,/trunc) ; str.time=wd.time_ccsds[ij[0]] str.file=filename[i] str.file_index=i str.xpix=ij[0] str.ypix=ij[1] str.x=wd.solar_x[str.xpix] str.y=wd.solar_y[str.ypix] str.thresh=threshold*float(ybin) str.ybin=ybin str.median=med_int[k[j]] str.intensity=int[k[j]] str.group=group[k[j]] IF n_tags(output) EQ 0 THEN output=str ELSE output=[output,str] ENDFOR ENDIF wd=0 ENDFOR IF NOT keyword_set(quiet) AND n_tags(output) NE 0 THEN BEGIN ng=max(output.group) print,'File index Event Time X-ind Y-ind Int (DN/s) npix' FOR i=1,ng DO BEGIN k=where(output.group EQ i,nk) getmax=max(output[k].intensity,imax) ii=k[imax] ; tstr=anytim2utc(/ccsds,output[ii].time,/time,/trunc) print,format='(5x,i5,3x,i5,3x,a8,2x,2i6,f12.1,i6)', $ output[ii].file_index,i,tstr,output[ii].xpix, $ output[ii].ypix,output[ii].intensity, $ nk ENDFOR ENDIF IF NOT keyword_set(quiet) THEN BEGIN print,'' print,format='("**NOTE: threshold adjusted for calibration and pixel-binning: ",f8.1," DN/s")',threshold*float(ybin) ENDIF END