;+ ; $Id: scc_mk_daily_med.pro,v 1.73 2022/06/10 22:17:28 secchia Exp $ ; ; Project : STEREO SECCHI ; ; Name : SCC_MK_DAILY_MED ; ; Purpose : This procedure generates an image by taking all the files of a given ; type (up to 40) for one day and finding the median value for each pixel. ; ; Explanation: ; ; Use : IDL> scc_mk_daily_med,'CAM','SC','YYYYMMDD',polar='POL' ; ; Inputs :CAM = cor1, cor2, hi1, hi2 ; SC = 'A', 'B' ; YYYYMMDD = date to be processed (nrlanytim) ; polar= '0', '120', '240', 'tbr' (tbr = total brightness), 'DBL' (type DOUBLE) ; ; Outputs : FITS file in $SECCHI_BKG/[ab]/daily_med/YYYYMM/ ; ; Keywords : SAVEDIR: specify your own directory to save file to ; FILES: specify your own file list (this may not work properly yet!!) ; /DOREBIN output 512x512 ; OUTSIZE= desired size of output (default is 1024x1024); native COR2 is 2048x2048 for prime mission ; /DOEXPTIMECORRECTION Do exptime correction in secchi_prep; otherwise just divide by exptime ; /NOSHIFT Do not shift data to match pointing of one image (applies to HI only) ; /NOSMOOTHPEAKS Do not do detection of dark pixels (implies /NOSHIFT and SMOOFILT=0). By default ; for COR2, smoothpeaks() is called for dark pixel detection and recording (only). ; SMOOFILT= Set equal to box size (default=17) for doing peak/valley removal; set by default for HI. ; ROLL= Use images with this sc_roll value. ; /DARKPIX Save dark pixel mask ; SOURCE = Telemetry source, either "rt" (realtime), "pb" ; (playback), or "lz" (level-zero, default). ; ; ; ; Calls from LASCO : ; ; Common : ; ; Restrictions: Need $SECCHI_BKG and appropriate permissions if not using savedir option ; ; Side effects: ; ; Category : DAILY display background stray light model ; ; Prev. Hist. : None. ; ; Written : Karl Battams, NRL/I2, Jan 2007 ; ; $Log: scc_mk_daily_med.pro,v $ ; Revision 1.73 2022/06/10 22:17:28 secchia ; removed cronus from mail statements ; ; Revision 1.72 2017/10/17 12:54:05 secchia ; corrected darkpix dir for new earth ; ; Revision 1.71 2017/09/08 15:51:52 nathan ; nr - try to fix scrolls ; ; Revision 1.70 2017/08/10 21:53:04 nathan ; fix incorrect filters for cor2 ; ; Revision 1.69 2015/12/03 14:14:35 secchia ; handles yearly directory for COR2 P0 images ; ; Revision 1.68 2015/07/23 15:32:09 secchia ; moved call to hi_mread_pointing after if file exists statement ; ; Revision 1.67 2015/07/20 18:19:58 secchia ; compare scroll post_conjunction to crota ; ; Revision 1.66 2015/07/20 17:06:36 secchia ; removed typo ; ; Revision 1.65 2015/07/20 16:59:29 hutting ; added post_conjunction keyword to get_stereo_roll ; ; Revision 1.64 2015/07/16 11:32:09 secchia ; added check for nmissing ; ; Revision 1.63 2014/09/03 16:28:56 secchia ; will write a daily median file if only one image is available for sidelobe operations ; ; Revision 1.62 2014/09/02 20:51:27 hutting ; reduced the number of images required if using beacon data ; ; Revision 1.61 2014/09/02 20:28:22 hutting ; useing beacon images if no others are available for sidelobe images ; ; Revision 1.60 2014/08/22 14:07:01 secchia ; added check for directory before changing ; ; Revision 1.59 2013/11/29 22:45:30 secchia ; ensure odd sizes not returned ; ; Revision 1.58 2013/11/29 13:22:16 secchia ; added source keyword ; ; Revision 1.57 2013/11/24 16:42:28 secchia ; nr - set /NOSMOOTHPEAKS if not /DARKPIX and not ishi; move darkpix criteria (perc) out of for loop ; ; Revision 1.56 2013/11/20 22:41:59 nathan ; Require /DARKPIX to save dark pixel mask (requires change in secchi_reduce) ; ; Revision 1.55 2013/11/20 17:50:39 nathan ; re-enable /SMOOFILT with updated smoothpeaks implementation ; ; Revision 1.54 2013/11/20 16:12:20 nathan ; disable smoofilt until it works ; ; Revision 1.53 2013/11/19 22:26:13 nathan ; add ROLL= keyword ; ; Revision 1.52 2012/08/23 20:59:00 secchib ; bug ; ; Revision 1.51 2012/08/23 15:38:07 secchib ; nr - return if SECCHI_P0 not available ; ; Revision 1.50 2012/05/16 21:30:15 secchib ; nr - change valid sc_roll range from 1 deg to 2 ; ; Revision 1.49 2012/05/16 20:42:10 nathan ; update roll for postd to account for B early 2007 ; ; Revision 1.48 2012/05/09 14:37:07 nathan ; fix SMOOFILT logic ; ; Revision 1.47 2012/05/08 18:55:00 nathan ; Add /SMOOFILT option; fix conversion to type integer for COR2 (bscale, bzero) ; ; Revision 1.46 2012/05/02 15:51:11 secchia ; nr - do file_exist check before mreadfits; update n_images in dpfi ; ; Revision 1.45 2012-04-20 22:46:05 nathan ; change change output type to UINT ; ; Revision 1.44 2012/04/20 22:40:38 nathan ; change criteria for using bscale because of differenc COR2 DBL sizes ; ; Revision 1.43 2012/03/22 15:43:57 secchia ; nr - smoothpeaks(debug=verbose) ; ; Revision 1.42 2012/03/19 20:07:07 secchia ; nr - fix typo in smoothpeaks section for cor2 ; ; Revision 1.41 2012/03/08 17:45:52 nathan ; fix bugs; add /NOSMOOTHPEAKS; refine setting of WCS values in hdr0 ; ; Revision 1.40 2012/03/07 21:21:35 nathan ; For matchpos, always choose image pointed most toward sun to match to; correct ; dark pixels before hi_align; record dark pixel stats for all images in ; /net/earth/data4/secchi; put keywords used in comment in hdr ; ; Revision 1.39 2011/12/27 19:37:50 nathan ; streamline summary.prog query ; ; Revision 1.38 2011/11/17 19:30:25 nathan ; fix conflict between files var and FILES= keyword ; ; Revision 1.37 2011/11/17 18:31:03 nathan ; moved set matchpos ; ; Revision 1.36 2011/11/16 19:43:11 nathan ; implement matchpos for HI ; ; Revision 1.35 2011/08/11 21:48:32 nathan ; fix imtyp for cor2 tbr ; ; Revision 1.34 2011/08/09 18:57:39 nathan ; add /DOEXPTIMECORRECTION option; turn on update_hdr option to secchi_prep ; ; Revision 1.33 2011/08/05 17:54:48 nathan ; update SC_ROLL in header; do not check file sc_roll if FILES= set ; ; Revision 1.32 2011/08/04 22:35:02 nathan ; Implement FILES= keyword so no other input is required, and generate filenames ; with rDEG indicating S/C roll (if GT 8 deg). ; ; Revision 1.31 2011/07/28 21:00:29 secchib ; removed stop ; ; Revision 1.30 2011/07/28 20:54:33 nathan ; fix /totalb use; fix s/c roll check ; ; Revision 1.29 2011/07/21 14:46:42 secchia ; nr - distribute images over whole day and change how screened for roll ; ; Revision 1.28 2011/07/20 19:10:38 secchia ; last changes not quite done ; ; Revision 1.27 2011/07/20 18:19:24 secchia ; nr - call COR2 dbl and polar with /nowarp ; ; Revision 1.26 2011/07/19 20:05:30 secchia ; nr - finish previous fix ; ; Revision 1.25 2011/07/19 18:29:49 secchia ; accept 1024 images for cor2 ; ; Revision 1.24 2010/11/18 18:41:52 secchib ; set imgtyp to img for cor2 dbls ; ; Revision 1.23 2010/05/12 14:47:13 nathan ; added OUTSIZE keyword ; ; Revision 1.22 2009/05/15 16:10:58 secchia ; nr - set TYPE= for scc_read_summary ; ; Revision 1.21 2009/04/09 12:03:40 mcnutt ; added cor2 pol=dbl backgrounds ; ; Revision 1.20 2009/03/19 16:26:26 mcnutt ; gets average roll for his roll= ( roll@date_cmd+roll@date_end ) /2 ; ; Revision 1.19 2008/12/11 19:18:51 nathan ; do not include zero values in median ; ; Revision 1.18 2008/08/06 19:16:58 nathan ; correct cor1 inp directory; fix bug if savedir=. ; ; Revision 1.17 2008/02/12 14:55:18 secchia ; nr - added /new to scc_putin_array call ; ; Revision 1.16 2008/02/11 15:31:06 secchia ; nr - fixed case where ng=0 ; ; Revision 1.15 2007/11/19 19:22:05 nathan ; use normal SECCHI_BKG directories ; ; Revision 1.14 2007/11/01 21:49:14 nathan ; added image roll check and dstart/stop correction ; ; Revision 1.13 2007/10/26 22:24:12 nathan ; Added /DOREBIN; added cadence, readtime, cleartim, crval, nmissing, cosmics ; to output header; tried to get dsatval right; use start of first image for ; DATE-OBS and midpoint is DATE-AVG; scale COR2 images to type INT before ; saving; save HI as type FLOAT; saves results in $SECCHI_BKG/../newbkg ; ; Revision 1.12 2007/09/25 22:55:52 nathan ; mostly done--still need to update scc_monthly_min.pro ; ; Revision 1.11 2007/09/25 16:00:11 nathan ; This has been re-written to utilize secchi_prep and the full SECCHI FITS header ; and to resolve Bug 224 ; ; Revision 1.10 2007/08/16 15:08:12 nathan ; made sure default output dir is group writable ; ; Revision 1.9 2007/07/25 14:15:05 reduce ; needed uppercase S/C identifier. Karl B ; ; Revision 1.8 2007/07/05 19:03:36 reduce ; Fixed 0-deg polarizer filename error. Karl B. ; ; Revision 1.7 2007/07/02 19:28:26 reduce ; Added a 'cd,orig' cmd. Karl B. ; ; Revision 1.6 2007/06/27 14:01:15 reduce ; Exits gracefully if too few files found. Karl B. ; ; Revision 1.5 2007/06/22 20:46:25 nathan ; cd,orig before return in Error case ; ; Revision 1.4 2007/06/21 19:20:00 reduce ; No longer crashes on missing images. Karl B. ; ; Revision 1.3 2007/06/18 20:48:57 reduce ; Couple of bug fixes. Should be ready for release now. Karl B. ; ; Revision 1.2 2007/02/01 20:43:52 reduce ; More mods and bug fixes. KB ; ; Revision 1.1 2007/01/30 20:29:21 reduce ; Initial Release by KB ; ;- ; ; PROCEDURE: ; For each image that satifies the selection conditions, (default naxis1=1024, ; filter and polarizer as requested), the median image is computed of ; the median value of all the images for a single day after being ; normalized to the median exposure time. ; ; If the number of images is less than 7, there is a second pass. ; In the second pass, images within +/- 2 days of the given day are used, ; up to 15 per day. ; ; ; Convert the telescope number into lower case ; And select the standard size parameters ; pro scc_mk_daily_med,tel,sc,Date0,polar=polar, DOREBIN=dorebin, SAVEDIR=savedir,FILES=files, $ VERBOSE=verbose, OUTSIZE=outsize, DOEXPTIMECORRECTION=doexptimecorrection, NOSHIFT=noshift, $ NOSMOOTHPEAKS=NOSMOOTHPEAKS, SMOOFILT=smoofilt, ROLL=roll, DARKPIX=darkpix, SOURCE=source IF keyword_set(FILES) THEN BEGIN jk=sccreadfits(files[0],h0,/nodata) tel=h0.detector sc=rstrmid(h0.obsrvtry,0,1) date0=h0.date_obs polar=1 IF tel EQ 'EUVI' THEN polar=trim(h0.wavelnth) ELSE $ IF tel EQ 'HI1' or tel EQ 'HI2' THEN polar='tbr' ELSE $ IF h0.seb_prog EQ 'DOUBLE' THEN polar='dbl' ELSE $ IF h0.polar GT 999 THEN polar='tbr' ELSE polar=trim(h0.polar) ENDIF ELSE IF (n_params() LT 3) then BEGIN PRINT,'Incorrect number of input parameters.' PRINT,' PRINT,'Usage: PRINT,' scc_mk_daily_med,CAM,S/C,DATE,POLAR PRINT,'EXAMPLE: PRINT,' IDL> scc_mk_daily_med,"hi1","a","20070115",polar="tbr"' PRINT,' --- OR --- PRINT,' IDL> scc_mk_daily_med,"cor2","b","20070222",polar="120"' PRINT,' return ENDIF IF (STRLOWCASE(tel) EQ 'cor1') THEN BEGIN PRINT,' PRINT,'**********************************************************' PRINT,' Please use cor1_mk_daily_med.pro PRINT,'**********************************************************' PRINT,' return ENDIF version='$Id: scc_mk_daily_med.pro,v 1.73 2022/06/10 22:17:28 secchia Exp $' len=strlen(version) version=strmid(version,1,len-2) src='lz' IF keyword_set(SOURCE) THEN src=source IF keyword_set(ROLL) THEN n_median_min = 2 ELSE $ n_median_min = 5 ; number of images needed to generate median cam = STRLOWCASE(tel) xaxis=1 yaxis=1 if not keyword_set(polar) then BEGIN PRINT,'Keyword "polar" needs to be defined. PRINT,'Options are as follows:' PRINT,' "tbr" -- total brightness (COR or HI) PRINT,' "0", "120" or "240" -- polarized (COR only) PRINT,' "dbl" -- total brightness (COR2) ENDIF pol=strlowcase(polar) polval=1001 ; total brightness IF pol EQ 'tbr' or pol EQ 'dbl' THEN imtyp='img' ELSE imtyp='seq' IF keyword_set(VERBOSE) THEN strn = ' Using:' ELSE strn = '' usep0=0 ishi=0 prog='Sequ' CASE cam OF 'cor1': BEGIN if pol EQ 'tbr' THEN BEGIN s=GETENV('SECCHI_P0')+'/'+strcompress(strlowcase(sc),/remove_all)+'/cor1/' pol_str='pTBr' ENDIF ELSE BEGIN s=GETENV('secchi')+'/'+src+'/L0/'+strcompress(strlowcase(sc),/remove_all)+'/seq/cor1/' pp=strcompress(pol,/remove_all) IF (float(pp) EQ 0) THEN pol_str='p000' ELSE pol_str='p'+pp polval=float(pol) ENDELSE pref='c1' ; filename prefix xaxis=512 yaxis=512 END 'cor2': BEGIN stdexptime=4. ; am I using this? dblexptime=6. ; am I using this? prog='Norm' if pol EQ 'tbr' THEN BEGIN s=GETENV('SECCHI_P0')+'/'+strcompress(strlowcase(sc),/remove_all)+'/cor2/' pol_str='pTBr' imtyp='pol' usep0=1 ENDIF ELSE if pol EQ 'dbl' THEN BEGIN s=GETENV('secchi')+'/'+src+'/L0/'+strcompress(strlowcase(sc),/remove_all)+'/img/cor2/' pol_str='dbTB' prog='Doub' ENDIF ELSE BEGIN s=GETENV('secchi')+'/'+src+'/L0/'+strcompress(strlowcase(sc),/remove_all)+'/seq/cor2/' pp=strcompress(pol,/remove_all) IF (float(pp) EQ 0) THEN pol_str='p000' ELSE pol_str='p'+pp polval=float(pol) ENDELSE pref='c2' ; filename prefix xaxis=1024 yaxis=1024 perc=0.08 ; for darkpix detection END 'hi1': BEGIN IF KEYWORD_SET(polar) THEN pol=strupcase(polar) ELSE pol='' s=GETENV('secchi')+'/'+src+'/L0/'+strcompress(strlowcase(sc),/remove_all)+'/img/hi_1/' pref='h1' pol='tbr' pol_str='pTBr' polval=0 xaxis=1024 yaxis=1024 ishi=1 perc=0.08 END 'hi2': BEGIN pol='' s=GETENV('secchi')+'/'+src+'/L0/'+strcompress(strlowcase(sc),/remove_all)+'/img/hi_2/' pref='h2' pol='tbr' pol_str='pTBr' polval=0 xaxis=1024 yaxis=1024 ishi=1 perc=0.04 END ELSE: BEGIN PRINT,'Unrecognized telescope code: '+cam RETURN END ENDCASE IF keyword_set(OUTSIZE) THEN outs=outsize ELSE outs=1024 ;if xaxis lt outs then outs=xaxis IF keyword_set(DOREBIN) THEN outs=512 img_type='' matchpos=1 maxshift=0. ; maximum pixel shift smoofac=17 IF keyword_set(SMOOFILT) THEN IF smoofilt GT 1 THEN smoofac=smoofilt smoofilt=keyword_set(SMOOFILT) IF keyword_set(NOSHIFT) or keyword_set(NOSMOOTHPEAKS) or ~ishi THEN matchpos=0 IF ~ishi and ~keyword_set(DARKPIX) THEN nosmoothpeaks=1 IF matchpos THEN smoofilt=1 IF keyword_set(NOSMOOTHPEAKS) THEN smoofilt=0 ut0=nrlanytim2utc(date0) date=utc2yymmdd(ut0,/yyyy) dt2=strmid(date,0,4)+'-'+strmid(date,4,2)+'-'+strmid(date,6,2) ; stop IF keyword_set(SAVEDIR) THEN BEGIN ; in case savedir= '.' cd,savedir spawn,'pwd',outp,/sh savedir=outp[0] ENDIF ; Find appropriate files... IF keyword_set(FILES) THEN BEGIN fnames = files ng=n_elements(files) good = indgen(ng) print,'Using ',trim(ng),' files for median.' ENDIF ELSE BEGIN ; Check if data exists yet... ; ; s=GETENV('secchi')+'/lz/L0/'+strcompress(strlowcase(sc),/remove_all)+'/'+imgdir ; got here ************** ; IF ~file_exist(s) THEN BEGIN ; message,s+' does not exist; returning',/info ; message,'Not doing '+tel+sc+' '+Date0+' '+polar,/info ; wait,10 ; return ; ENDIF ; path=s+date+'/' ; f=file_search(path+'./*fts') ; sz=size(f) ; IF (sz(0) NE 0) THEN BEGIN ; CD,path,curr=orig ; f=file_search((path+'./*fts') ; sz=size(f) ; IF (sz(0) EQ 0) THEN BEGIN ; PRINT,'No files in directory for '+cam+' telescope on '+date ; CD,orig ; RETURN ; ENDIF ; ENDIF ELSE BEGIN ; PRINT,'No directory for '+cam+' telescope on '+date ; return ; ENDELSE ; help,dt2 ;; use beacon during sidelobe operations ; spw=strmid(f,17,1) ; z=where(spw ne 7,cnt) ; if cnt eq 0 then begin ; nobeacon=0 ; n_median_min = 2 ; endif else $ nobeacon=1 summary=scc_read_summary(DATE=dt2,SPACECRAFT=sc,TELESCOPE=tel,TYPE=imtyp, TOTALB=usep0, nobeacon=nobeacon, source=src) if nobeacon eq 0 then begin if cam eq 'hi2' or cam eq 'hi1' then begin if cam eq 'hi2' then exta=3 else exta=2 ut0=nrlanytim2utc(date0) utm1=ut0 & utp1=ut0 utm1.mjd=utm1.mjd-exta utp1.mjd=utp1.mjd;+exta date1=utc2str(utm1,/date_only) date3=utc2str(utp1,/date_only) summary=scc_read_summary(DATE=[date1,date3],SPACECRAFT=sc,TELESCOPE=tel,TYPE=imtyp, TOTALB=usep0, nobeacon=nobeacon, source=src) if cam eq 'hi1' and (n_elements(summary) ge 2 ) then summary=summary(where(summary.compr eq 'RICE')) endif IF (n_elements(summary) ge 2 )then begin xaxis=min(summary.xsize) yaxis=min(summary.xsize) endif endif PRINT,'Found a total of ',strcompress(n_elements(summary)),' files for the day in the summary file.' IF (n_elements(summary) LE 5 and nobeacon GT 1) THEN BEGIN PRINT,'Based on the summary file, there do not appear to be enough files for making a median for today on this telescope.' CD,orig return ENDIF IF datatype(summary) EQ 'INT' THEN BEGIN ; CD,orig return ENDIF ; Either an individual COR polarizer angle is picked, or a HI image (no polarizer) or total brightness is picked. if nobeacon ge 1 then good = where(summary.VALUE EQ polval and summary.XSIZE mod outs EQ 0 and summary.YSIZE mod outs EQ 0 and summary.PROG EQ prog, ng) if nobeacon eq 0 then good = where(summary.VALUE EQ polval and summary.XSIZE mod xaxis EQ 0 and summary.YSIZE mod yaxis EQ 0, ng) if cam eq 'cor2' and pol_str eq 'dbTB' then begin tmp=where(summary[good].exptime ge dblexptime-.5 and SUMMARY[good].exptime le dblexptime+.5) good=good[tmp] endif print,'Found ',trim(ng),' appropriate images for ',dt2,' from summary file.' IF ng EQ 0 THEN BEGIN print,'ERROR: There are not enough acceptable images to make a good daily median image. Returning.' wait,5 ; CD,orig return ENDIF ELSE begin if cam eq 'cor2' and pol_str eq 'pTBr' then $ fnames=getenv('SECCHI_P0')+'/'+strlowcase(sc)+'/cor2/'+strmid(summary[good].FILENAME,0,4)+'/'+strmid(summary[good].FILENAME,0,8)+'/'+summary[good].FILENAME else $ fnames=sccfindfits(summary[good].FILENAME) endelse ENDELSE print,'Screening roll for acceptability...' ; - Make sure roll is not outside min and max value ut=anytim2utc(dt2) ;uta=replicate(ut,5) ;uta.mjd=uta.mjd + indgen(5)-2 ; spawn,'pwd',/SH ; First read in all headers for screening and position info. chk=where(file_exist(fnames),ng) fnames=fnames[chk] jnk=sccreadfits(fnames,h,/nodata) IF keyword_set(ROLL) THEN medianscroll=roll $ ELSE BEGIN scrolls1=-1*get_stereo_roll(h.date_obs,sc,/verbose,/post_conjunction) ww=where(scrolls1 GT 300,nww) IF nww GT 0 THEN scrolls1[ww]=scrolls1[ww]-360. ww=where(scrolls1 LT -170,nww) IF nww GT 0 THEN scrolls1[ww]=scrolls1[ww]+180. if n_elements(scrolls1) gt 1 then medianscroll=median(scrolls1) else medianscroll=scrolls1 ENDELSE minr=medianscroll-1 maxr=medianscroll+1 help, minr,maxr sc_rolls=h.crota nmiss=h.nmissing IF keyword_set(FILES) THEN nrg=ng ELSE BEGIN ;if cam eq 'hi1' or cam eq 'hi2' then ;sc_rolls=((-1*get_stereo_roll(anytim2utc(h.date_cmd),sc,/verbose))+(-1*get_stereo_roll(anytim2utc(h.date_end),sc,/verbose)))/2 wok=where( sc_rolls GT minr and sc_rolls LT maxr, nrg, COMPLEMENT=badroll,NCOMPLEMENT=nbad) IF nbad GT 0 THEN BEGIN print,'Omitting ',fnames[badroll],'because SC_ROLL=',sc_rolls[badroll] wait,5 ENDIF if wok[0] gt -1 then begin fnames=fnames[wok] h=h[wok] endif wok=where( nmiss eq 0, nrg, COMPLEMENT=badmiss,NCOMPLEMENT=nbad) IF nbad GT 0 THEN BEGIN print,'Omitting ',fnames[badmiss],'because nmissing =',nmiss[badmiss] wait,5 ENDIF if wok[0] gt -1 then begin fnames=fnames[wok] h=h[wok] endif ENDELSE ; make sure images distributed throughout the day IF nrg GT 40 THEN BEGIN indx=findgen(40)*nrg/40. fnames=fnames[indx] h=h[indx] ENDIF n=n_elements(fnames) hdr0=h[n/2] if n LT n_median_min then begin print,'ERROR: There are not enough acceptable images to make a good daily median image. Returning.' wait,5 ; CD,orig return endif ELSE PRINT,'Using ',trim(n),' files for median.' IF (matchpos) and (ishi) THEN BEGIN dlm=get_delim() hisrch=getenv('SCC_DATA')+dlm+'hi'+dlm+'pnt_'+strupcase(tel+sc)+'_'+dt2+'*fts' hifn=file_search(hisrch) IF hifn[0] EQ '' THEN BEGIN message,hisrch+' not found; setting matchpos=0',/info wait,2 matchpos=0 ENDIF ELSE BEGIN hihdrs=hi_mread_pointing(hifn) strtimes=intersection(strmid(hihdrs.filename,0,15),strmid(h.filename,0,15),ind1=wf) hihdrs=hihdrs[wf] ; we only want to shift toward the sun, so we will use the min crval1 value, ; and otherwise the median value for hdr0. w0=where(abs(hihdrs.crval1) EQ min(abs(hihdrs.crval1))) ;hdrn=h[w0[0]] ;hdr0=str_copy_tags(hdrn,hihdrs[w0[0]]) sumdif=float(hihdrs[w0[0]].nx/hdr0.naxis1) help,sumdif hdr0.crval1a = hihdrs[w0[0]].crval1a hdr0.crval1 = hihdrs[w0[0]].crval1 hdr0.filename = h[w0[0]].filename SMOOFILT hdr0.crval2a = median(hihdrs.crval2a) hdr0.pc1_1a = median(hihdrs.pc1_1a) hdr0.pc1_2a = median(hihdrs.pc1_2a) hdr0.pc2_1a = median(hihdrs.pc2_1a) hdr0.pc2_2a = median(hihdrs.pc2_2a) hdr0.cdelt1a = median(hihdrs.cdelt1a*sumdif) hdr0.cdelt2a = median(hihdrs.cdelt2a*sumdif) hdr0.pv2_1a = median(hihdrs.pv2_1a) hdr0.crval2 = median(hihdrs.crval2) hdr0.pc1_1 = median(hihdrs.pc1_1) hdr0.pc1_2 = median(hihdrs.pc1_2) hdr0.pc2_1 = median(hihdrs.pc2_1) hdr0.pc2_2 = median(hihdrs.pc2_2) hdr0.cdelt1 = median(hihdrs.cdelt1*sumdif) hdr0.cdelt2 = median(hihdrs.cdelt2*sumdif) hdr0.pv2_1 = median(hihdrs.pv2_1) hdr0.xcen = median(hihdrs.xcen) hdr0.ycen = median(hihdrs.ycen) hdr0.crota = median(hihdrs.crota) hdr0.ins_x0 = median(hihdrs.ins_x0) hdr0.ins_y0 = median(hihdrs.ins_y0) hdr0.ins_r0 = median(hihdrs.ins_r0) hdr0.ravg = median(hihdrs.ravg) hdr0.naxis=2 hdr0.ipsum=1 device,window_state=ws IF ws[2] THEN wset,2 ELSE BEGIN window,2,xsi=600,ysi=800 plotprep ENDELSE !p.multi=[0,1,2,0,1] plot,hihdrs.crval1,title='CRVAL1 for '+strupcase(tel+sc)+'_'+dt2 plot,hihdrs.crval2,title='CRVAL2 for '+strupcase(tel+sc)+'_'+dt2 !p.multi=0 ENDELSE ENDIF ELSE $ hdr0=h[n/2] get_utc,utn datacube= fltarr(outs,outs,n) crotas= fltarr(n) biases= fltarr(n) exptimes= dblarr(n) offsetcrs= fltarr(n) if n gt 1 then cadences= fltarr(n-1) else cadences= fltarr(n) readtimes= dblarr(n) cleartimes= dblarr(n) pitches= fltarr(n) yaws= fltarr(n) nmissings= lonarr(n) nsats= lonarr(n) cosmicss= lonarr(n) midpoint=n/2 n_used=n ; used later bimsum=bytarr(outs,outs) nbim=0 ; This is a temporary bug fix for the COR total-B summary files last_char=strmid(fnames[0],strlen(fnames[0])-1,1) if (last_char) NE 's' then fnames=fnames+'s' ; end bug fix noexptime=~keyword_set(DOEXPTIMECORRECTION) IF keyword_set(VERBOSE) THEN quiet=0 ELSE quiet=1 pntroll=medianscroll IF pntroll LT 0 THEN titleangle = 360+pntroll ELSE titleangle = pntroll postd='' IF abs(medianscroll) GT 16 THEN postd = 'r'+STRING(round(titleangle),format='(I3.3)') IF abs(medianscroll) GT 16 THEN stop ; roll of B was up to 15 until June 2007 fname0 = 'd'+pref+strupcase(sc)+'_'+pol_str+'_'+strmid(date,2,6)+postd+'.fts' fac=0 ;wait,5 FOR m=0,n-1 DO BEGIN if file_exist(fnames[m]) THEN $ fn=fnames[m] ELSE $ if file_exist(fnames[m-1]) THEN BEGIN fn=fnames[m-1] PRINT,'Could not find file: ',fnames[m] PRINT,'Will use previous file instead...' ENDIF ELSE BEGIN PRINT,'OK, I can not find the previous file either. I give up.' CD,orig return ENDELSE ;print,fn IF (pol EQ 'tbr' and cam EQ 'cor2') THEN BEGIN ; tbr is already prepped im = sccreadfits(fn,hdr) ; im is type float IF hdr.naxis1 EQ 2048 THEN BEGIN hdr.dstart1=1 hdr.dstop1 =2048 ; some early TBr images are incorrect ENDIF im = SCC_PUTIN_ARRAY(im,hdr,outs,/new) ss='' k=0 WHILE ss NE 'EXPTIMEs' DO BEGIN k=k+1 ss=strmid(hdr.comment[k],0,8) ENDWHILE print,hdr.comment[k] parts=strsplit(hdr.comment[k],' ,',/extract) recexp=avg(float(parts[1:3])) effexp=1.0 ; TBr from SECCHI_P0 are effectively 1 second dsatval=30000. ; approximate for 6sec exposure, from 1 image example (20071005_025230_0B4c2B.fts) ; collect dark pixel stats IF ~keyword_set(NOSMOOTHPEAKS) THEN BEGIN ;perc=0.1 ; 0.5 im1=smoothpeaks(im,smoofac,perc*median(im),/display,WHEREs=w,/dark, DEBUG=~quiet) xyouts,10,10,hdr.filename,/device IF (SMOOFILT) THEN BEGIN ; correct bad pixels before align_image smoomsg='pixels darker than median filter(,'+trim(smoofac)+')+'+trim(perc)+'*median replaced' message,smoomsg,/info im=im1 ENDIF ; making mask of dark pixels bim=uintarr(outs,outs) bim[w]=1 bimsum=bimsum+bim nbim=nbim+1 ENDIF ; ENDIF ELSE BEGIN secchi_prep,fn, hdr,im0, /NOCALFAC, /calimg_off, EXPTIME_OFF=noexptime, OUTSIZE=outs, SILENT=quiet, /nowarp ; image is corrected for IP including IPSUM ; Do not correct for EXPTIME in secchi_prep because of weighting done in hi_correction.pro recexp=hdr.exptime ;help,recexp ; wait,2 IF noexptime THEN im0=im0/recexp effexp=1.0 dsatval=hdr.dsatval IF dsatval LE 0 THEN BEGIN IF strmid(pref,0,1) EQ 'h' THEN ccdsat=14000. ELSE ccdsat=15000. dsatval = long(ccdsat * hdr.n_images * (2^(hdr.IPSUM - 1))^2) IF strmid(pref,0,1) NE 'h' THEN dsatval = dsatval<60000 ENDIF ; collect dark pixel stats IF ~keyword_set(NOSMOOTHPEAKS) THEN BEGIN ;perc=0.1 ; 0.5 im1=smoothpeaks(im0,smoofac,perc*median(im0),/display,WHEREs=w,/dark, DEBUG=~quiet) xyouts,10,10,hdr.filename,/device IF (SMOOFILT) THEN BEGIN ; correct bad pixels before align_image smoomsg='pixels darker than median filter(,'+trim(smoofac)+')+'+trim(perc)+'*median replaced' message,smoomsg,/info im0=im1 ENDIF ; making mask of dark pixels bim=uintarr(outs,outs) bim[w]=1 bimsum=bimsum+bim nbim=nbim+1 ENDIF IF (matchpos) THEN BEGIN ut0=anytim2utc(hdr.date_avg) IF hdr.ravg EQ -881 and utn.mjd-ut0.mjd GT 21 THEN BEGIN openw,evlun,'errmsg.txt',/get_lun printf,evlun,'******* mk_daily_med Error *******' printf,evlun,'' printf,evlun,'HI pointing file not found for ',hdr.date_avg printf,evlun,'Not making ',fname0 printf,evlun,'' close,evlun free_lun,evlun print,'' spawn,'cat errmsg.txt',/sh wait,5 ;maillist='nathan.rich@nrl.navy.mil lynn.simpson.ctr@nrl.navy.mil' spawn,'mailx -s '+hdr.obsrvtry+'_'+date+'_mk_daily_med nathan.rich@nrl.navy.mil < errmsg.txt',/sh return ENDIF im=hi_align_image(hdr,im0,hdr0,/zerom,/verbose) ; hdr is modified in hi_align_image.pro maxshift=maxshift>hdr.jitrsdev ENDIF ELSE im=im0 ENDELSE thistai=anytim2tai(hdr.date_obs) IF (m GT 0) THEN cadences[m-1]= thistai-lasttai lasttai=thistai readtimes[m]= hdr.readtime cleartimes[m]= hdr.cleartim pitches[m]= hdr.crval1 yaws[m]= hdr.crval2 nmissings[m]= hdr.nmissing cosmicss[m]= hdr.cosmics crotas[m]=hdr.crota biases[m]=hdr.biasmean exptimes[m]=recexp offsetcrs[m]=hdr.offsetcr ; bias that has been subtracted nsats[m]=hdr.datasat datacube[*,*,m]=im ;--Datacube is normalized to 1 second ; ; stop if (m EQ midpoint) THEN outhdr0 = hdr if (m EQ 0) THEN dateobstai=thistai if (m EQ 0) THEN date_obs=hdr.date_obs ENDFOR ; This is experimental but what I'm trying to do here is remove from the datacube any images ; that would not make a good median image. So I start by removing images that are rolled more ; than 5-sigma of the median roll value for that day. Later I will remove pixels that are ; freakishly bright compared to the others. Note that 5-sigma might not be tight enough... we'll see... ;crota_stdv=stddev(hdrs.crota) ;crota_med=median(hdrs.crota) ;good=where( abs(hdrs.crota - crota_med) LT (5*crota_stdv) ) crota_stdv=stddev(crotas) crota_med=median(crotas) if n_elements(good) gt 2 then good=where( abs(crotas - crota_med) LT (5*crota_stdv) ) help,good if (n_elements(good) NE n) THEN BEGIN if (n_elements(good) GT 0) THEN datacube=datacube[*,*,good] ELSE BEGIN message,'ERROR! Something went wrong with range of CROTA in selected images...' ENDELSE ENDIF ; stop if n_Elements(good) gt 1 then begin exptime_stdv =stdev(exptimes[good],exptime_mean) bias_stdv =stdev(biases[good],bias_mean) offsetcr_stdv =stdev(offsetcrs[good],offsetcr_mean) crota_stdv =stdev(crotas[good],crota_mean) cadence_stdv =stdev(cadences[good],cadence_mean) readtime_stdv =stdev(readtimes[good],readtime_mean) cleartime_stdv=stdev(cleartimes[good], cleartime_mean) crval1_stdv =stdev(pitches[good], crval1_mean) crval2_stdv =stdev(yaws[good], crval2_mean) cosmics_stdv =stdev(cosmicss[good], cosmics_mean) endif nsat =min(nsats) help,exptime_mean if n_Elements(good) gt 1 then dsatval=float(dsatval/exptime_mean) ; normalize this also to 1 second datendtai =anytim2tai(hdr.date_end) date_avg =utc2str(tai2utc(dateobstai + (datendtai-dateobstai)/2)) IF (matchpos) or n_Elements(good) eq 1 THEN BEGIN outhdr0=hdr0 crota_out=hdr0.crota crval1out=hdr0.crval1 crval2out=hdr0.crval2 rotmsg=' shifted to match '+hdr0.filename v1msg=rotmsg v2msg=rotmsg ENDIF ELSE BEGIN crota_out=crota_mean crval1out=crval1_mean crval2out=crval2_mean v1msg='=avg over N_IMAGES; stdev='+trim(crval1_stdv) v2msg='=avg over N_IMAGES; stdev='+trim(crval2_stdv) rotmsg='=avg over N_IMAGES; stdev='+trim(crota_stdv) ENDELSE outhdr0.nmissing=total(nmissings[good]) median_array=fltarr(outs,outs) if n_Elements(good) gt 1 then begin for j=0,outs-1 DO BEGIN for k=0,outs-1 DO BEGIN pixel_row = datacube[j,k,*] ; Get rid of zero-value pixels ; stop w = where(pixel_row[0,0,*] GT 0) if (n_elements(w) GT 1 ) THEN median_array[j,k]=median(pixel_row[0,0,w]) ; Now try to eliminate extreme-value pixels. ; Note that 10-sigma might be too high... ; pix_med = median(pixel_row) ; pix_sd = stddev(pixel_row) ; ok=where( abs(pixel_row - pix_med) LT (10*pix_sd) ) endfor endfor endif else median_array=datacube ; ; Do final smoothing to get rid of planets and comets ; IF keyword_set(SMOOFILT) THEN BEGIN ;IF 0 THEN BEGIN message,'Doing final smoothing with median_filter...',/info boxs=30 te2=smoothpeaks(median_array,boxs,0.05,/bright,/dark,/display, DEBUG=~quiet) median_array=te2 ENDIF ;IF keyword_set(DOREBIN) THEN median_array=rebin(temporary(median_array),512,512) ; This is handled with putin_array previously. -nr ;if (xaxis GT 1024) THEN median_array=rebin(median_array,1024,1024) maxmin,median_array maxmedian=max(median_array) dsatval=maxmedian bscale=1.0 bzero=0.0 bmsg='Data is type FLOAT' IF ~noexptime or n_elements(good) eq 1 THEN median_array_out=median_array>0 ELSE $ IF dsatval LT 64000/exptime_mean THEN BEGIN ; COR ; This is copied from reduce_level_1.pro: ;bscale = (scalemax-scalemin)/65536 ;bzero=bscale*32769 ;help,bscale,bzero ;nz=where(b NE 0) ;bout=fltarr(hdr.naxis1,hdr.naxis2) ;bout[nz] = ROUND(((b[nz]scalemin)-bzero)/bscale) ;bout = FIX(bout) bscale= 1./exptime_mean bzero=bscale*32768 median_array_out= fix(round( ((median_array>000,outhdr0,verbose=verbose,satmax=dsatval) outhdr1 = rem_tag2(outhdr0,'TIME_OBS') outhdr = struct2fitshead(outhdr1,/allow_crota,/dateunderscore) ;IF nsat LE 0 THEN it=where(median_array_out GE dsatval, nsat) ;fxaddpar,outhdr,'DSATVAL',dsatval ;fxaddpar,outhdr,'DATASAT',nsatload help,dsatval,nsat,bscale,bzero get_utc,dte,/ecs fxaddpar,outhdr,'DATE-OBS',date_obs,' Start of first exposure in N_IMAGES' fxaddpar,outhdr,'DATE-END',hdr.date_end,' End of last exposure in N_IMAGES' fxaddpar,outhdr,'DATE-AVG',date_avg,' Midpoint between OBS and END' if n_Elements(good) gt 1 then fxaddpar,outhdr,'CADENCE', median(cadences),' sec (median); stdev='+trim(cadence_stdv) if n_Elements(good) gt 1 then fxaddpar,outhdr,'READTIME',readtime_mean,' sec; stdev='+trim(readtime_stdv) if n_Elements(good) gt 1 then fxaddpar,outhdr,'CLEARTIM',cleartime_mean,' sec; stdev='+trim(cleartime_stdv) fxaddpar,outhdr,'CRVAL1',crval1out,v1msg fxaddpar,outhdr,'CRVAL2',crval2out,v2msg if n_Elements(good) gt 1 then fxaddpar,outhdr,'COSMICS',cosmics_mean,' stdev='+trim(cosmics_stdv) fxaddpar,outhdr,'DATE',dte IF ishi and ~noexptime THEN fxaddpar,outhdr,'DISTCORR','T','hi_desmear has been applied' ; need a flag so scc_mkframe knows not to run hi_prep on this image fxaddpar,outhdr,'BUNIT','DN/sec' fxaddpar,outhdr,'BSCALE',bscale fxaddpar,outhdr,'BZERO',bzero,bmsg fxaddpar,outhdr,'HISTORY','Data values normalized to 1 sec' fxaddpar,outhdr,'EXPTIME', 1.0,'effective exptime; actual avg='+trim(exptime_mean)+'; stdev='+trim(exptime_stdv) if n_Elements(good) gt 1 then fxaddpar,outhdr,'BIASMEAN',bias_mean,'=avg over N_IMAGES' if n_Elements(good) gt 1 then fxaddpar,outhdr,'BIASSDEV',bias_stdv,' Over N_IMAGES' if n_Elements(good) gt 1 then fxaddpar,outhdr,'OFFSETCR',offsetcr_mean,'=avg over N_IMAGES; stdev='+trim(offsetcr_stdv) fxaddpar,outhdr,'CROTA',crota_out,'deg; '+rotmsg fxaddpar,outhdr,'SC_ROLL',medianscroll fxaddpar,outhdr,'N_IMAGES',n_used,' scc_mk_daily_med.pro' ;stop sdir1 = GETENV('SECCHI_BKG')+'/'+strlowcase(sc)+'/daily_med/' sdir2 = sdir1+strmid(date,0,6) if not file_exist(sdir2) THEN BEGIN PRINT,'Directory ',sdir2, ' does not exist.' cmd='mkdir -p -m 775 '+sdir2 PRINT,cmd spawn,cmd,/sh ENDIF IF keyword_set(SAVEDIR) THEN sdir = savedir ELSE sdir = sdir2 CD,sdir PRINT,'Writing daily median fits file: '+sdir+'/'+fname0 print,' using '+string(n_used)+' images' wait,2 fxaddpar,outhdr,'FILENAME',fname0 fxaddpar,outhdr,'HISTORY',version fxaddpar,outhdr,'COMMENT','%SCC_MK_DAILY_MED: header values from middle image unless noted' IF matchpos THEN fxaddpar,outhdr,'COMMENT',smoomsg IF matchpos THEN fxaddpar,outhdr,'JITRSDEV',maxshift,'Max hi_align_image pixel shift in x or y' wkwd=[keyword_set(DOEXPTIMECORRECTION), ~matchpos, keyword_set(FILES), keyword_set(NOSMOOTHPEAKS), keyword_set(SMOOFILT)] kwds=['/DOEXPTIMECORRECTION','/NOSHIFT','/FILES','/NOSMOOTHPEAKS','/SMOOFILT'] wkw=where(wkwd) IF wkw[0] GE 0 THEN fxaddpar,outhdr,'COMMENT','Keywords used: '+arr2str(kwds[wkw]) WRITEFITS,fname0,median_array_out,outhdr ;stop IF keyword_set(DARKPIX) and ~keyword_set(NOSMOOTHPEAKS) THEN BEGIN ;if 0 then begin ; record mask of dark pixels, one per month dpnm= pref+strupcase(sc)+'_darkpixels_'+strmid(date,0,6)+'.fts' dpfi='/net/earth/secchi/data/darkpix/'+strlowcase(sc)+'/'+dpnm IF file_exist(dpfi) THEN BEGIN dpim=readfits(dpfi,dph) lastn=fxpar(dph,'N_IMAGES') ENDIF ELSE BEGIN lastn=0 dpim=0L dph=outhdr fxaddpar,dph,'comment','each count indicates pixel in one original image was darker ' fxaddpar,dph,'comment','than median filter(,'+trim(fac)+')-0.1*median' ENDELSE fxaddpar,dph,'n_images',lastn+nbim writefits,dpfi,long(bimsum)+dpim,dph ENDIF ;IF not keyword_set(FILES) THEN CD,orig ;IF floor(r) NE last_r THEN GOTO, beginning RETURN END