;+ ; NAME ; ; EIS_REPLACE_MISSING ; ; PROJECT ; ; Hinode/EIS ; ; EXPLANATION ; ; This routine takes the output structure from EIS_GETWINDATA, and ; calculates replacement values for the intensity and error of ; 'missing' pixels. The replacement values can either be placed ; directly into the WINDATA array (using /REPLACE), or the new ; intensity and error arrays can be output separately. ; ; The intensity value of a missing pixel is derived from the ; intensities of non-missing neigbouring pixels in the detector-Y ; direction. A hierarchy system is used that accounts for various cases ; of neighbouring pixels being also flagged missing. E.g., if ; method 1 is not possible due to a neighbour pixel being missed, ; then method 2 is tried. If this also is not possible due to a ; missing pixel, then method 3 is applied, etc. If none of the 5 ; methods is possible then the pixel remains flagged as being ; missing. The methods are given below, where we denote the ; missing pixel by index i and the intensity of a pixel as P(i) ; ; 1. Use the average of the two neighbouring pixels (P(i-1), ; P(i+1)). ; ; 2. Use the weighted average of one of the neighbouring pixels, and the ; next-neighbouring pixel that is on the opposite side of the ; missing pixel (P(i-1) and P(i+2), or P(i+1) and P(i-2)). ; ; 3. Use the weighted average of one of the neighbouring pixels, and the ; next-next-neighbouring pixel that is on the opposite side of ; the missing pixel (P(i-1) and P(i+3), or P(i+1) and P(i-3)). ; ; 4. Use the average of the two next-neighbouring pixels (P(i-2) ; and P(i+2)). ; ; 5. Use the value from one of the missing neighbours (P(i-1) or ; P(i+1)). ; ; The error values assigned to the replaced pixels are derived by ; first performing a linear fit to the error^2 vs. intensity points ; for all the non-missing data points. For fit parameters a and b, ; where ; ; error^2 = a + b*intensity ; ; then an error estimate is obtained from this relationship by ; inserting the newly-replaced intensity value. ; ; Since the five methods described above represent increasingly ; poor approximations to the intensity of the replaced pixels, an ; additional factor is included in the error value such that: ; ; error_new = error * frac(k) ; ; where k denotes which of the interpolation methods was used. The ; default values of frac lie between 1.0 and 1.3. ; ; INPUTS ; ; WINDATA The output structure from EIS_GETWINDATA. If the ; /REPLACE keyword is set, then the intensity and error ; array tags will be modified. ; ; KEYWORDS ; ; QUIET If set, then no information messages are printed to the ; screen. ; ; REPLACE If set, then replace the .INT and .ERR array tags of ; WINDATA with the newly computed intensity and error ; arrays. If not set, then WINDATA remains unaltered by ; EIS_REPLACE_MISSING. ; ; OPTIONAL INPUTS ; ; FRAC A five element array whose significance is explained ; above. ; ; OPTIONAL OUTPUTS ; ; WDOUT The intensity array from WINDATA, with missing pixel ; values replaced. ; ; ERROUT The error array from WINDATA, with missing pixel values ; replaced. ; ; OUTPUTS ; ; If /REPLACE is set, then WINDATA will be returned with modified ; INT and ERR array tags. Otherwise WINDATA will remain the same. ; ; CALLS ; ; EIS_INT_ERR_FIT, EIS_REPLACE_ERRORS ; ; HISTORY ; ; Ver.1, 30-Mar-2009, Peter Young ; Ver.2, 14-Apr-2009, Peter Young ; corrected bug related to pixels with negative values; tidied ; up routine, introducing new functions erm_calc_average_3 and ; erm_calc_average_4 for working out the averages. ; Ver.3, 16-Apr-2009, Peter Young ; error interpolation method was not correct for broad ; wavelength windows so introduced calls to eis_int_err_fit and ; eis_replace_errors to do this correctly. ; Ver.4, 29-Apr-2009, Peter Young ; updated header and added additional information printing ; Ver.5, 15-Apr-2010, Peter Young ; the routine now checks for any exposures with a zero exposure ; time. If found, all the pixels in the exposures will be set to ; missing. ; Ver.6, 25-Aug-2010, Peter Young ; corrected bug when there are multiple exposures per raster ; position (NEXP_PRP>1). ; Ver.7, 5-Nov-2010, Peter Young ; the routine now checks for data columns that are anomalously ; bright and flags them as missing. ; Ver.8, 20-Dec-2010, Peter Young ; removed change introduced in Ver. 7 as it removed good data in ; one case. ; Ver.9, 13-Apr-2011, Peter Young ; added check to make sure eis_int_err_fit ran correctly. If not ; then no data are replaced. ; Ver.10, 10-Jan-2013, Peter Young ; if the /noabs keyword was used in eis_prep (so intensity given ; in DN units), then the error array will not be ; interpolated. Note that windata.hdr.cal_abs flags status of ; absolute calibration. ; Ver.11, 13-May-2025, Peter Young ; eis_prep now produces a correct error array for the /noabs ; option, so I've modified the routine to correctly handle these ; data (previously the err array was not modified for /noabs data). ; Ver.12, 07-Aug-2025, Peter Young ; Updated calls to eis_int_err_fit and eis_replace_errors. ;- FUNCTION erm_calc_average_3, wd, method=method, miss_map=miss_map, chck=chck ; ; WD Input intensity array. ; METHOD Interpolation method ; MISS_MAP The map of missing pixels ; CHCK The check array for whether method has worked. ; siz=size(wd) ny=siz[3] wd_avg=fltarr(siz[1],siz[2],siz[3]) chck=bytarr(siz[1],siz[2],siz[3])+1b ; CASE method OF 1: BEGIN ; -- average neighbours wd_avg[*,*,1:ny-2]=0.5*(wd[*,*,0:ny-3]+wd[*,*,2:ny-1]) chck[*,*,1:ny-2]=miss_map[*,*,0:ny-3]+miss_map[*,*,2:ny-1] END ; 2: BEGIN ; -- weighted average of neighbour & n-neighbour wd_avg[*,*,1:ny-3]=(2.*wd[*,*,0:ny-4]+wd[*,*,3:ny-1])/3. chck[*,*,1:ny-3]=miss_map[*,*,0:ny-4]+miss_map[*,*,3:ny-1] END ; 3: BEGIN ; -- weighted average of neighbour & n-neighbour wd_avg[*,*,2:ny-2]=(wd[*,*,0:ny-4]+2.*wd[*,*,3:ny-1])/3. chck[*,*,2:ny-2]=miss_map[*,*,0:ny-4]+miss_map[*,*,3:ny-1] END ; 4: BEGIN ; -- weighted average of neighbour & n-n-neighbour wd_avg[*,*,1:ny-4]=(3.5*wd[*,*,0:ny-5]+wd[*,*,4:ny-1])/4.5 chck[*,*,1:ny-4]=miss_map[*,*,0:ny-5]+miss_map[*,*,4:ny-1] END ; 5: BEGIN ; -- weighted average of neighbour & n-n-neighbour wd_avg[*,*,3:ny-2]=(wd[*,*,0:ny-5]+3.5*wd[*,*,4:ny-1])/4.5 chck[*,*,3:ny-2]=miss_map[*,*,0:ny-5]+miss_map[*,*,4:ny-1] END ; 6: BEGIN ; -- average of n-neighbours wd_avg[*,*,2:ny-3]=0.5*(wd[*,*,0:ny-5]+wd[*,*,4:ny-1]) chck[*,*,2:ny-3]=miss_map[*,*,0:ny-5]+miss_map[*,*,4:ny-1] END ; 7: BEGIN ; -- upper neighbour wd_avg[*,*,0:ny-2]=wd[*,*,1:ny-1] chck[*,*,0:ny-2]=miss_map[*,*,1:ny-1] END ; 8: BEGIN ; -- lower neighbour wd_avg[*,*,1:ny-1]=wd[*,*,0:ny-2] chck[*,*,1:ny-1]=miss_map[*,*,0:ny-2] END ; ENDCASE return,wd_avg END ;-------------------------- FUNCTION erm_calc_average_4, wd, method=method, miss_map=miss_map, chck=chck siz=size(wd) ny=siz[3] wd_avg=fltarr(siz[1],siz[2],siz[3],siz[4]) chck=bytarr(siz[1],siz[2],siz[3],siz[4])+1b ; CASE method OF 1: BEGIN ; -- average neighbours wd_avg[*,*,1:ny-2,*]=0.5*(wd[*,*,0:ny-3,*]+wd[*,*,2:ny-1,*]) chck[*,*,1:ny-2,*]=miss_map[*,*,0:ny-3,*]+miss_map[*,*,2:ny-1,*] END ; 2: BEGIN ; -- weighted average of neighbour & n-neighbour wd_avg[*,*,1:ny-3,*]=(2.*wd[*,*,0:ny-4,*]+wd[*,*,3:ny-1,*])/3. chck[*,*,1:ny-3,*]=miss_map[*,*,0:ny-4,*]+miss_map[*,*,3:ny-1,*] END ; 3: BEGIN ; -- weighted average of neighbour & n-neighbour wd_avg[*,*,2:ny-2,*]=(wd[*,*,0:ny-4,*]+2.*wd[*,*,3:ny-1,*])/3. chck[*,*,2:ny-2,*]=miss_map[*,*,0:ny-4,*]+miss_map[*,*,3:ny-1,*] END ; 4: BEGIN ; -- weighted average of neighbour & n-n-neighbour wd_avg[*,*,1:ny-4,*]=(3.5*wd[*,*,0:ny-5,*]+wd[*,*,4:ny-1,*])/4.5 chck[*,*,1:ny-4,*]=miss_map[*,*,0:ny-5,*]+miss_map[*,*,4:ny-1,*] END ; 5: BEGIN ; -- weighted average of neighbour & n-n-neighbour wd_avg[*,*,3:ny-2,*]=(wd[*,*,0:ny-5,*]+3.5*wd[*,*,4:ny-1,*])/4.5 chck[*,*,3:ny-2,*]=miss_map[*,*,0:ny-5,*]+miss_map[*,*,4:ny-1,*] END ; 6: BEGIN ; -- average of n-neighbours wd_avg[*,*,2:ny-3,*]=0.5*(wd[*,*,0:ny-5,*]+wd[*,*,4:ny-1,*]) chck[*,*,2:ny-3,*]=miss_map[*,*,0:ny-5,*]+miss_map[*,*,4:ny-1,*] END ; 7: BEGIN ; -- upper neighbour wd_avg[*,*,0:ny-2,*]=wd[*,*,1:ny-1,*] chck[*,*,0:ny-2,*]=miss_map[*,*,1:ny-1,*] END ; 8: BEGIN ; -- lower neighbour wd_avg[*,*,1:ny-1,*]=wd[*,*,0:ny-2,*] chck[*,*,1:ny-1,*]=miss_map[*,*,0:ny-2,*] END ; ENDCASE return,wd_avg END ;---------------------- PRO eis_replace_missing, windata, wdout=wdout, errout=errout, $ quiet=quiet, replace=replace, $ frac=frac, info_array=info_array IF n_params() LT 1 THEN BEGIN print,'Use IDL> eis_replace_missing, windata [, /replace, wdout=, errout=, /quiet,' print,' frac=, info_array= ]' print,'Examples:' print,' IDL> eis_replace_missing, windata, /replace' print,' IDL> eis_replace_missing, windata, wdout=wdout,errout=errout' return ENDIF IF (windata.hdr.slit_ind EQ 1) OR (windata.hdr.slit_ind EQ 3) THEN BEGIN message,/info,/cont,'This routine only works for narrow slit data.' return ENDIF ; ; If the /noabs or /phot keyword were given to eis_prep, then the effective ; area curve was not applied to the intensities. This modifies how the data ; replacement is done, and I flag these cases with noabs=1. ; ;; noabs=0b ;; photon=0b ;; chck=strpos(windata.hdr.calstat,'ABS') ;; IF chck LT 0 THEN BEGIN ;; noabs=1b ;; ENDIF ELSE BEGIN ;; chck=strpos(windata.hdr.calstat,'PHOT') ;; IF chck GE 0 THEN photon= ;; ENDELSE IF n_elements(frac) EQ 0 THEN frac=[1.0,1.2,1.2,1.3,1.3] nexp_prp=windata.hdr.nexp_prp ; ; When an exposure time is zero eis_prep currently does not set the ; data for that exposure to be missing. The following piece of code ; sets the zero exposures to be missing. Note that it's ; necessary to check if there are multiple exposures per raster ; position. ; IF nexp_prp EQ 1 THEN BEGIN i=where(windata.exposure_time EQ 0.,ni) IF ni NE 0 THEN BEGIN FOR j=0,ni-1 DO BEGIN k=i[j] windata.err[*,k,*]=windata.missing ENDFOR ENDIF ENDIF ELSE BEGIN FOR k=0,nexp_prp-1 DO BEGIN exposure_time=windata.exposure_time[*,k] i=where(exposure_time EQ 0.,ni) IF ni NE 0 THEN BEGIN FOR j=0,ni-1 DO BEGIN ij=i[j] windata.err[*,ij,*,k]=windata.missing ENDFOR ENDIF ENDFOR ENDELSE wdin=windata.int errin=windata.err wvl=windata.wvl missing=windata.missing siz=size(wdin) ny=siz[3] nl=siz[1] wdout=wdin errout=errin i=where(errin EQ missing,ni) ; ; MISS_MAP stores the locations of the missing pixels. ; CASE siz[0] OF 3: BEGIN miss_map=bytarr(siz[1],siz[2],siz[3]) info_array=bytarr(siz[1],siz[2],siz[3]) END 4: BEGIN miss_map=bytarr(siz[1],siz[2],siz[3],siz[4]) info_array=bytarr(siz[1],siz[2],siz[3],siz[4]) END ELSE: BEGIN message,/info,/cont,'Data arrays have non-standard size. They should have 3 or 4 dimensions.' return END ENDCASE ; IF ni NE 0 THEN miss_map[i]=1b ; ; Derive line fit parameters that describe the relation between error ; and intensity. ; f=eis_int_err_fit(windata,/quiet,noabs=noabs,photon=photon) ; ; The value of 255 is used to identify pixels that are missing, but ; which do not get a replacement intensity. I initially set each ; missing pixel to 255 and then change the values depending on which ; of the replacement methods works for that pixel. ; info_array[i]=255b ; ; The loop below takes the array through the eight different averaging ; methods, with the least effective first and the most effective ; last. The average values are stored in WD_AVG, while CHCK indicates ; whether the averaging method was successful. E.g., if CHCK[k]>0 then ; the method was not successful in the sense that one of the averaging ; pixels was flagged as missing, or the method would not work (e.g., ; for edge pixels). ; FOR j=8,1,-1 DO BEGIN CASE siz[0] OF 3: wd_avg=erm_calc_average_3(wdout,method=j,miss_map=miss_map,chck=chck) 4: wd_avg=erm_calc_average_4(wdout,method=j,miss_map=miss_map,chck=chck) ENDCASE ; CASE 1 OF j EQ 7 OR j EQ 8: frac_val=frac[4] j EQ 6: frac_val=frac[3] j EQ 5 OR j EQ 4: frac_val=frac[2] j EQ 3 OR j EQ 2: frac_val=frac[1] j EQ 1: frac_val=frac[0] ENDCASE ; err_interp=eis_replace_errors(wvl,wd_avg,f,noabs=noabs,photon=photon) ; k=where(chck[i] EQ 0b,nk) IF nk NE 0 THEN BEGIN wdout[i[k]]=wd_avg[i[k]] errout[i[k]]=err_interp[i[k]] * frac_val info_array[i[k]]=byte(j) ENDIF ENDFOR ; ; For all remaining pixels, set intensity and error to be missing. ; k=where(info_array EQ 255b,nk) IF nk NE 0 THEN BEGIN wdout[k]=missing errout[k]=missing ENDIF IF NOT keyword_set(quiet) THEN BEGIN i=where(info_array EQ 255b,ni) print,format='("%age of pixels not replaced: ",f5.1)',float(ni)/float(n_elements(wdin))*100. ; i=where(info_array EQ 8b OR info_array EQ 7b,ni) print,format='("%age of pixels replaced, method 5: ",f5.1)',float(ni)/float(n_elements(wdin))*100. ; i=where(info_array EQ 6b,ni) print,format='("%age of pixels replaced, method 4: ",f5.1)',float(ni)/float(n_elements(wdin))*100. ; i=where(info_array EQ 4b OR info_array EQ 5b,ni) print,format='("%age of pixels replaced, method 3: ",f5.1)',float(ni)/float(n_elements(wdin))*100. ; i=where(info_array EQ 3b OR info_array EQ 2b,ni) print,format='("%age of pixels replaced, method 2: ",f5.1)',float(ni)/float(n_elements(wdin))*100. ; i=where(info_array EQ 1b,ni) print,format='("%age of pixels replaced, method 1: ",f5.1)',float(ni)/float(n_elements(wdin))*100. ENDIF ; ; Replace intensity and error arrays in WINDATA ; IF keyword_set(replace) THEN BEGIN windata.int=wdout windata.err=errout IF NOT keyword_set(quiet) THEN BEGIN message,/info,/cont,'The INT and ERR tags of WINDATA have been modified.' ENDIF ENDIF ELSE BEGIN IF NOT keyword_set(quiet) THEN BEGIN print,'% EIS_REPLACE_MISSING: WINDATA has not been modified.' print,' Use the keyword /REPLACE to update the intensity and error arrays' print,' in WINDATA.' ENDIF ENDELSE ; ; Tidy up ; junk=temporary(wd_avg) junk=temporary(miss_map) END