; file: findsplinealign.pro ; init: Jan 6 2016 Rob Rutten Deil ; last: Jan 14 2016 Rob Rutten Deil ; note: uses isml routines not in newer IDL nor in MacOS IDL ;+ pro findsplinealign,alignlabel=alignlabel,aligndir=aligndir,$ splinesmooth=splinesmooth,cutbad=cutbad,$ inforfile=inforfile,inrevfile=inrevfile,$ outforfile=outforfile,outrevfile=outrevfile,$ plotdir=plotdir,$ angfor=angfor,shiftfor=shiftfor,scalefor=scalefor,$ nxfor=nxfor,nyfor=nyfor,$ angrev=angrev,shiftrev=shiftrev,scalerev=scalerev,$ nxrev=nxrev,nyrev=nyrev,$ verbose=verbose ; NAME: ; findsplinealign ; PURPOSE: ; put spline approximations through output of findalignfitscubes.pro ; CALL: ; see above ; OPTIONAL KEYWORD INPUTS: ; alignlabel = string identifying comparison (no spaces, default '') ; path and file name defaults: ; aligndir: /tmp/align/' ; inforfile: aligndir+'align_arrays/formet-'+alignlabel+'.dat' ; inrevfile: aligndir+'align_arrays/revmet-'+alignlabel+'.dat' ; outforfile: aligndir+'align_arrays/forspline-'+alignlabel+'.dat' ; outrevfile: aligndir+'align_arrays/revspline-'+alignlabel+'.dat' ; plotdir: aligndir+'align_plots/ + similar .ps file names ; splinesmooth = 3-elem vector: shift, rotate, scale; default [0,100,100] ; value 0: use IDL's estimation algorithm to derive smoothness ; value ne 0: typically 1 = through every point, 100 = straight line ; play with this parameter according to your data! ; cutbad = multiplier to rms to cut off outliers (bad seeing?) (default 5) ; verbose = 1/0: more/less terminal messaging (default 0) ; OUTPUTS: ; arrays in optional variables and as columns in outforfile and outrevfile: ; angfor,shiftfor,scalefor,nxfor,nyfor = forward parameter values ; angrev,shiftrev,scalerev,nxrev,nyrev = reverse parameter values ; plots of Metcalf input and spline output, outliers marked ; OPTIONAL OUTPUTS: ; plots on screen if verbose=1 ; RESTRICTIONS: ; none that I am aware off ; REMARKS: ; see also explanation file: rjlib/00-README/align-pipeline.txt ; HISTORY: ; Jan 9 2016 RR: start ;- ; keyword defaults if (n_elements(splinesmooth) eq 0) then splinesmooth=[0,100,100] if (n_elements(cutbad) eq 0) then cutbad=5. if (n_elements(alignlabel) eq 0) then alignlabel='' if (n_elements(aligndir) eq 0) then aligndir='/tmp/align/' if (n_elements(inforfile) eq 0) then $ inforfile=aligndir+'align_arrays/formet-'+alignlabel+'.dat' if (n_elements(inrevfile) eq 0) then $ inrevfile=aligndir+'align_arrays/revmet-'+alignlabel+'.dat' if (n_elements(outforfile) eq 0) then $ outforfile=aligndir+'align_arrays/forspline-'+alignlabel+'.dat' if (n_elements(outrevfile) eq 0) then $ outrevfile=aligndir+'align_arrays/revspline-'+alignlabel+'.dat' if (n_elements(plotdir) eq 0) then plotdir=aligndir+'align_plots/' if (n_elements(verbose) eq 0) then verbose=0 ; check splinesmooth if (n_elements(splinesmooth) ne 3) then begin print,' ##### abort: splinesmooth not 3-elem: shift, rotate, scale' return endif ; flag spaces in aligndir (as Catherine Fischer's Mac disk) file1=strjoin(strsplit(aligndir,' ',/extract),'\ ') ; make directories if not yet existing spawn,'mkdir -p '+aligndir spawn,'mkdir -p '+aligndir+'align_arrays/' spawn,'mkdir -p '+plotdir ; read forward Metcalf-fit parameters from file readcol,inforfile,$ metcalftimesfor,ntfor,$ angmetfor,xshiftmetfor,yshiftmetfor,xscalemetfor,yscalemetfor,$ nxmetfor,nymetfor ; read reverse Metcalf-fit parameters from file readcol,inrevfile,$ metcalftimesrev,ntrev,$ angmetrev,xshiftmetrev,yshiftmetrev,xscalemetrev,yscalemetrev,$ nxmetrev,nymetrev ; check metcalftimes if (total(metcalftimesrev) ne total(metcalftimesfor)) then begin print,' ##### abort: metcalftimesfor differs from metcaftimesrev' return endif else metcalftimes=metcalftimesfor ; check nt if (ntrev[0] ne ntfor[0]) then begin print,' ##### abort: ntrev ne ntrev' return endif else nt=ntfor[0] ; get and check ntmet sizemet=size(metcalftimes) ntmet=sizemet[1] if (ntmet lt 2 or ntmet gt nt) then begin print,' ##### abort: ntmet = '+nostr(ntmet)+'; nt = '+ntostr(nt) return endif ; check metcalf sizes if (min(nxmetfor) ne max(nxmetfor) or min(nymetfor) ne max(nymetfor) or $ min(nxmetrev) ne max(nxmetrev) or min(nymetrev) ne max(nymetrev))$ then begin print,' ##### abort: metcalf nx and ny sizes are not time-constant' return endif ; make indexed input array metarr=[[[xscalemetfor],[yscalemetfor],$ [angmetfor],$ [xshiftmetfor],[yshiftmetfor]],$ [[xscalemetrev],[yscalemetrev],$ [angmetrev],$ [xshiftmetrev],[yshiftmetrev]]] metmuckarr=metarr ; declare output arrays angfor=fltarr(nt) shiftfor=fltarr(nt,2) scalefor=fltarr(nt,2) angrev=fltarr(nt) shiftrev=fltarr(nt,2) scalerev=fltarr(nt,2) ; initiate metit=indgen(ntmet) allit=indgen(nt) xtitle='time step' thick=1 senselabel=['forward','reverse'] diagname=['xscale','yscale','angle','xshift','yshift'] ; this order re excess allitbad=[-1] ; splinesmooth should be set as [shift,rotate,scale] splinesm=[splinesmooth[2],splinesmooth[2],$ splinesmooth[1],$ splinesmooth[0],splinesmooth[0]] ; iterate spline determination to find and remove outliers (bad seeing?) niter=2 ;RR keep it simple; for more cutbad must be increased per iiter for iiter=0,niter-1 do begin ; start loop over forward,reverse for isense=0,1 do begin ; start loop over metcalf parameters splinall=fltarr(5,nt) splinmet=fltarr(ntmet) for idiag=0,4 do begin ; assign metdata=reform(metarr[*,idiag,isense]) metmuckdata=reform(metmuckarr[*,idiag,isense]) ; get spline approximation splinsol=imsl_cssmooth(metcalftimes,metmuckdata,$ smpar=float(ntmet)*splinesm[idiag]) splinmet[*]=imsl_spvalue(metcalftimes[metit],splinsol) splinall[idiag,*]=imsl_spvalue(allit,splinsol) ; non-final iteration: replace outliers by spline-fit value if (iiter ne niter-1) then begin ; find outliers ("bad" points) in this parameter diff=abs(metmuckdata-splinmet) momdiff=moment(diff,sdev=sdev) itbad=where(diff gt cutbad*(iiter/2.+1)*sdev,nbad) ; add to earlier ones for application in subsequent parameters if (nbad ne 0) then begin allitbad=[allitbad,itbad] if (allitbad[0] eq -1) then allitbad=allitbad[1:nbad] endif ; replace all bad it values so far by spline fit value in Metcalf array if (allitbad[0] ne -1) then $ metmuckarr[allitbad,idiag,isense]=splinmet[allitbad] endif ; final iteration: make plot (against array index since no timing info) if (iiter eq niter-1) then begin psfilename=plotdir+senselabel[isense]+'-'+$ diagname[idiag]+'-spline-'+alignlabel+'.ps' xrange=[-nt/10.,(nt-1)+nt/10.] datarange=[min(metdata),max(metdata)] dataspan=datarange[1]-datarange[0] yrange=[datarange[0]-0.1*dataspan,datarange[1]+0.1*dataspan] openpsplot,psfilename,fontsize=8,thick=thick plot,allit,reform(splinall[idiag,*]),$ position=[0.20,0.15,0.95,0.90],/normal,$ xticklen=0.03,yticklen=0.03/1.6,$ xtitle=xtitle,ytitle=diagname[idiag],$ xrange=xrange,xstyle=1,yrange=yrange,ystyle=1 ; overplot Metcalf input (before outlier removal) oplot,metcalftimes,metdata,psym=2,symsize=0.7 ; mark the bad points by triangles if (allitbad[0] ne -1) then $ oplot,metcalftimes[allitbad],metdata[allitbad],psym=5,symsize=1.5 ; add alignlabel xyouts,0.4,0.92,/norm,senselabel[isense]+' '+alignlabel ; plot done closepsplot,psfilename,opengv=verbose endif ; end plot output endfor ; end loop over idiag ; final iteration: assign and write the spline outputs if (iiter eq niter-1) then begin if (isense eq 0) then begin scalefor=[[reform(splinall[0,*])],[reform(splinall[1,*])]] angfor=reform(splinall[2,*]) shiftfor=[[reform(splinall[3,*])],[reform(splinall[4,*])]] nxfor=intarr(nt)+nxmetfor[0] nyfor=intarr(nt)+nymetfor[0] writecol,outforfile,$ shiftfor[*,0],shiftfor[*,1],angfor,scalefor[*,0],scalefor[*,1],$ nxfor,nyfor,fmt='(5F9.4,2I6)' if (verbose ne 0) then print,' ===== wrote '+outforfile endif if (isense eq 1) then begin scalerev=[[reform(splinall[0,*])],[reform(splinall[1,*])]] angrev=reform(splinall[2,*]) shiftrev=[[reform(splinall[3,*])],[reform(splinall[4,*])]] nxrev=indgen(nt)*0+nxmetrev[0] nyrev=indgen(nt)*0+nymetrev[0] writecol,outrevfile,$ shiftrev[*,0],shiftrev[*,1],angrev,scalerev[*,0],scalerev[*,1],$ nxrev,nyrev,fmt='(5F9.4,2I6)' if (verbose ne 0) then print,' ===== wrote '+outrevfile endif endif endfor ; end loop isense = forward and reverse endfor ; end iteration loop to remove outliers ; print how many outliers also if verbose=0) if (allitbad[0] ne -1) then $ print, ' ===== it values bad outliers (repeats = how often detected): ',$ allitbad[sort(allitbad)] end ; ================= test per hyper-C ========================== aligndir='/media/rutten/SSTDATA/alldata/SST/2014-06-21-quiet/align_sdo_sst/' alignlabel='AIA1700-SST8542_0' findsplinealign,aligndir=aligndir,alignlabel=alignlabel ; inspect (copy to other terminal; cd to aligndir/align_plots) ; killgv' ; gvallps xshift end