; file: sdo_writepairspline.pro ; init: Feb 24 2016 Rob Rutten Deil ; last: Dec 12 2020 Rob Rutten Deil ; note: findimshift_tiled.pro finds shift im2 to im1 but all pairs here are ; reversed because for some historical reason sdo_makeslavefitscube.pro ; applies negative shift values. So pair[0] is put on pair[1]. ; note: muck parameter choices marked by !!!! ;+ pro sdo_writepairspline,datadir,wavpair,trimbox=trimbox,$ aligncenter=aligncenter,xalpair=xalpair,delaymin=delaymin,$ boxcarlabel=boxcarlabel,blink=blink,verbose=verbose ; NAME: ; sdo_writepairspline.pro ; ; PURPOSE: ; measure, spline shifts and plot shifts between two SDO wavelengths ; ; INPUTS: ; datadir = string path/dir with /cubes, usually .../center ; wavpair: [wavnr1, wavnr2] = RR_channelnrs, 1 to be put on 2 ; ; OPTIONAL KEYWORD INOUTS: ; trimbox = align only over selected subfield [xmin,ymin,xmax,ymax] ; aligncenter = 1/0: make checkplot xalpair from cubesxal ; xalpair = RR paircode for the final check plot when /aligncenter ; delaymin = delaymin: sampling delay first after second in minutes ; boxcarlabel = identifier for boxcar-smeared input ('' if none) ; blink = 1/0: yes/no showex image and shifted image for it=0 ; verbose = 1/0: yes/no align plots on screen, printout ; ; OUTPUTS: ; files in subdir ../centerdrifts with plots and spline structures ; ; EXAMPLES: ; below at end of file ; ; HISTORY: ; Feb 24 2016 RR: start ; Apr 15 2016 RR: better mucks, spline outlier removal ; Apr 19 2016 RR: spline IMSL > Vitas for Mac-challenged ; Aug 24 2016 RR: refined AIA mucking ; Aug 28 2016 RR: trimbox option ; Feb 21 2017 RR: revamp EUV cross-alignment ; Sep 21 2017 RR: added mgaussnorm ; Apr 8 2020 RR: spline fits to splin2diter.pro ; Jun 12 2020 RR: boxcarlabel, time-smear and delay euvanchor ; Jun 17 2020 RR: heightdiff, 95% conf values ; Jun 28 2020 RR: image pair mucking > sdo_muckimagepair.pro ;- ; RR channel numbers and AIA telescopes ; 94,131,171,193,211,304,335,1600,1700,4500,mag,cont,dop ; 1 2 3 4 5 6 7 8 9 10 11 12 13 ; 4 1 3 2 2 4 1 3 3 3 0 0 0 ; answer no-parameter query if (n_params(0) lt 2) then begin sp,sdo_writepairspline return endif ; get current dir path; \ flag spaces in path if present (bloody Mac?) ; https://groups.google.com/forum/#!topic/comp.lang.idl-pvwave/Lo10H5XtN80 cd,current=thisdir thisdir=strjoin(strsplit(thisdir,' ',/extract),'\ ') ; default keyword values if (n_elements(datadir) eq 0) then datadir=thisdir+'/center/' if (n_elements(trimbox) eq 0) then trimbox=[-1,-1,-1,-1] if (n_elements(aligncenter) eq 0) then aligncenter=0 if (n_elements(xalpair) eq 0) then xalpair=[2,9] if (n_elements(delaymin) eq 0) then delaymin=0 if (n_elements(boxcarlabel) eq 0) then boxcarlabel='' if (n_elements(blink) eq 0) then blink=0 if (n_elements(verbose) eq 0) then verbose=0 ; make driftsdir if not there already driftsdir=datadir+'/../driftscenter'+boxcarlabel+'/' spawn,'mkdir -p '+driftsdir ; define permitted level2 fitsfilename strings (RR_channelnr = index+1) sdowavs=['94','131','171','193','211','304','335',$ '1600','1700','4500','mag','cont','dop'] nsdowavs=n_elements(sdowavs) sdotel=[4,1,3,2,2,4,1,3,3,3,0,0,0] ; SDO telescopes (AIA 1-4, HMI=0) ; check wavpair for iwavnr=0,1 do begin if (wavpair[iwavnr] lt 1 or wavpair[iwavnr] gt 13) then begin print,' ##### sdo_writepairspline abort: wrong wavpair = ', wavpair return endif endfor if (wavpair[1] eq wavpair[0]) then begin print,' ##### sdo_writepairspline abort: wav 2 = wav 1' return endif ; get wav names wavnames=strarr(2) for iwav=0,nsdowavs-1 do begin if (iwav eq wavpair[0]-1) then wavnames[0]=sdowavs[iwav] if (iwav eq wavpair[1]-1) then wavnames[1]=sdowavs[iwav] endfor ; get the two cube files ;RR add "/quote" for \ flagged spaces in paths (eg Mac "My Passport" disk) ; select cubessmear if boxcarlabel non-empty cd,datadir if (boxcarlabel ne '') then begin if (file_test('cubessmear'+boxcarlabel)) then $ cubesdir='/cubessmear'+boxcarlabel+'/' else begin print,' ##### sdo_writepairspline abort: boxcarlabel present but no files' STOP endelse endif else cubesdir='/cubes/' if (aligncenter ne 0 and wavpair[0] eq xalpair[0] and $ wavpair[1] eq xalpair[1]) then cubesdir='/cubesxal/' files1=file_search(datadir+cubesdir+'*'+wavnames[0]+'.fits',$ count=nfiles1,/quote) files2=file_search(datadir+cubesdir+'*'+wavnames[1]+'.fits',$ count=nfiles2,/quote) ; check cube file existence if (nfiles1 ne 1 or nfiles2 ne 1) then begin print,' ##### sdo_writepairspline abort: pair '+$ wavnames[0]+'-'+wavnames[1]+' has no or too many cube files; skip' goto,SKIPTHISPAIR endif ; now we have the files for this wavpair infile1=files1[0] infile2=files2[0] ; get cube1 dimensions and more from fits header inhead1=headfits_rr(infile1) inhead1size=(1+fix(n_elements(inheader)/36.))*2880 nx1=fxpar(inhead1,'naxis1') ny1=fxpar(inhead1,'naxis2') nt1=fxpar(inhead1,'naxis3') bitpix=fxpar(inhead1,'bitpix') cadence=fxpar(inhead1,'cadence') xcen=fxpar(inhead1,'xcen') ycen=fxpar(inhead1,'ycen') ; smear wavpairs 304-1600, 304-1700, 304-mag temporally ; irrespective of boxcarmin smearing in sdo_writeallpairsplines ; idea: 304 sits above network indirectly, possibly non-E retardation smearlabel='' spearpx=0 if (wavpair[0] eq 6 and $ (wavpair[1] eq 8 or wavpair[1] eq 9 or wavpair[1] eq 11)) then begin smearmin=5 ; temporal boxcar width in minutes ; !!!!! RR-muck choice smearit=fix((smearmin*60.)/cadence) ; for even nr boxcarfitscube adds 1 if (smearmin gt 0.001) then begin ; check that smearit doesn't exceed nt/2 if (smearit gt nt1/2) then begin smearit=fix(nt1/2.) smearmin=smearit*cadence/60. print,' ##### sdo_writepairspline: smearit = '+trimd(smearit)+$ ' exceeds half nt = '+trimd(nt1)+' smearmin set to'+trimd(smearmin,1) endif smearlabel='tsm'+trimd(smearmin,1) rannr=fix(abs(randomn(seed))*10000) smearfile1=$ '/tmp/sdosmeareuv'+strmid(string(rannr+1E5,format='(i6)'),2)+'.fits' smearfile2=$ '/tmp/sdosmearuv'+strmid(string(rannr+1E5,format='(i6)'),2)+'.fits' boxcarfitscube,infile1,smearfile1,boxcar=[0,0,smearit] boxcarfitscube,infile2,smearfile2,boxcar=[0,0,smearit] infile1=smearfile1 infile2=smearfile2 endif endif ; set endian bigendian=1 ; open input file 1 for assoc get_lun, unit_in1 if (bigendian) then openr,unit_in1,infile1,/swap_if_little_endian $ else openr,unit_in1,infile1 if (bitpix eq -32) then inassoc1=assoc(unit_in1,fltarr(nx1,ny1),inhead1size) if (bitpix eq 16) then inassoc1=assoc(unit_in1,intarr(nx1,ny1),inhead1size) if (bitpix eq 8) then inassoc1=assoc(unit_in1,bytarr(nx1,ny1),inhead1size) ; get cube2 dimensions and file datatype from the fits header inhead2=headfits_rr(infile2) inhead2size=(1+fix(n_elements(inheader)/36.))*2880 nx2=fxpar(inhead2,'naxis1') ny2=fxpar(inhead2,'naxis2') nt2=fxpar(inhead2,'naxis3') bitpix2=fxpar(inhead2,'bitpix') ; check they have the same bitpix if (bitpix2 ne bitpix) then begin print,' ##### sdo_writepairspline abort: '+$ wavnames[0]+'-'+wavnames[1]+' cubes do not have same bitpix' return endif ; open input file 2 for assoc get_lun, unit_in2 if (bigendian) then openr,unit_in2,infile2,/swap_if_little_endian $ else openr,unit_in2,infile2 if (bitpix eq -32) then inassoc2=assoc(unit_in2,fltarr(nx2,ny2),inhead2size) if (bitpix eq 16) then inassoc2=assoc(unit_in2,intarr(nx2,ny2),inhead2size) if (bitpix eq 8) then inassoc2=assoc(unit_in2,bytarr(nx2,ny2),inhead2size) ; check common dimensions (they should be same size already) if (nx1 ne nx2) then print,' ===== OOPS?? nx1 # nx2 wavpair '+trim(wavpair) if (ny1 ne ny2) then print,' ===== OOPS?? ny1 # ny2 wavpair '+trim(wavpair) if (nt1 ne nt2) then print,' ===== OOPS?? nt1 # nt2 wavpair '+trim(wavpair) nx=min([nx1,nx2]) ny=min([ny1,ny2]) nt=min([nt1,nt2]) ; set trimbox if requested if (trimbox[0] ne -1) then begin selxmin=trimbox[0] selymin=trimbox[1] selxmax=trimbox[2] selymax=trimbox[3] if (selxmax gt nx-1 or selymax gt ny-1) then begin print,' ###### sdo_writepairspline abort: selxmax or selymax too large' return endif end else begin selxmin=0 selymin=0 selxmax=nx-1 selymax=ny-1 endelse ; convert delaymin to delayit, define delaylabel delayit=fix((delaymin*60.)/cadence+0.5) if (delaymin ne 0) then delaylabel='tdt'+trimd(delaymin,1) else delaylabel='' ; set up shifts determination shifts=fltarr(nt-delayit,2) confx=fltarr(nt-delayit) confy=fltarr(nt-delayit) timearr=dblarr(nt-delayit) arcsecpx=0.6 pxtile=50 precision=0.05 heightdiff=0 ; overwritten per muck pair below with best estimate timediff=delayit*cadence ;RR Dec 11 2020 I don't timediff anymore, would also clash in derot ; start loop to get shift per time step for it=0,nt-1-delayit do begin ; read image pair im1in=inassoc1[it+delayit] im2in=inassoc2[it] ; cut im1=im1in[selxmin:selxmax,selymin:selymax] im2=im2in[selxmin:selxmax,selymin:selymax] ; muck this image pair astile=pxtile*arcsecpx sdo_muckimagepair,im1,im2,wavpair,$ muckim1,muckim2,astile,tilecut,heightdiff ; find shift by center field tiling thisshift=$ findimshift_tiled(muckim1,muckim2,$ timediff=timediff,heightdiff=heightdiff,$ arcsecpx=arcsecpx,pxtile=pxtile,$ tilecut=tilecut,precision=precision) ; optional showex blinker for first it (zoom in) if (blink ne 0 and it eq 0) then begin print,' ##### hit showex "quit this" to continue' wait,2 shiftim1=shift_img(im1in,-thisshift) showex,sqrt(histo_opt_rr(shiftim1>0,1.E-3)),$ sqrt(histo_opt_rr(im2in>0,1.E-3)) endif ; done shifts[it,0]=thisshift[0] shifts[it,1]=thisshift[1] confx[it]=thisshift[2] confy[it]=thisshift[3] endfor ; end it loop over image pairs with time ; done with reading: free the input and output files free_lun,unit_in1,unit_in2 ;; ; INSERT: write shifts file for tests ;; shiftsfile=driftsdir+'/shifts-'+ntostr(wavpair[0])+ntostr(wavpair[1])+'.dat' ;; writecol,shiftsfile,timearr,shifts[*,0],shifts[*,1],$ ;; fmt='(D20.2,2f15.3)' ; ======== spline fit to temporal behavior ; set time variables starttime=fxpar(inhead1,'starttim') cadence=fxpar(inhead1,'cadence') starttai=anytim2tai(starttime) timearr=starttai+indgen(nt-delayit)*cadence stoptai=starttai+(nt-delayit-1)*cadence stoptime=anytim2utc(stoptai,/ccsds) ; rework for plot labeling starttim2=str_replace(starttime,'T',' ') startdate=strmid(starttime,0,10) starthour=strmid(starttime,11,5) stophour=strmid(stoptime,11,5) timestring=startdate+' '+starthour+' - '+stophour ; setup spline computation shiftx=reform(shifts[*,0]) shifty=reform(shifts[*,1]) ; spline fit with 1x outlier removal (Apr 8 2020 shunted to splin2diter) niter=2 cutsdev=3. nt=n_elements(shiftx) smooth1=5.E5*nt/25. ; Apr 12 2020 5.E6 smooth2=5.E6*nt/25. ; Apr 12 2020 5.E8 splin2diter,timearr,shiftx,shifty,timearrleft,$ ax,bx,cx,dx,ay,by,cy,dy,$ smooth1=smooth1,smooth2=smooth2,niter=niter,cutsdev=cutsdev,$ verbose=verbose ; get spline curves at all sample points splinx=splinapp_get(timearr,timearrleft,ax,bx,cx,dx) spliny=splinapp_get(timearr,timearrleft,ay,by,cy,dy) ; find and print remaining rms deviations from splines match,timearr,timearrleft,indstay diffx=abs(shiftx[indstay]-splinx[indstay]) diffy=abs(shifty[indstay]-spliny[indstay]) momx=moment(diffx,sdev=sdevx) momy=moment(diffy,sdev=sdevy) if (verbose) then begin print,' -----'+trimd(nt-n_elements(indstay))+' outliers removed' print,' ----- rms remaining points from spline in x'+trimd(sdevx) print,' ----- rms remaining points from spline in y'+trimd(sdevy) endif ; ========= make spline plots ; plot label pairname=wavnames[0]+'-'+wavnames[1] plotlabel=pairname+' tels '+$ trim(sdotel[wavpair[0]-1])+'-'+trim(sdotel[wavpair[1]-1])+' ' if (aligncenter ne 0 and $ wavpair[0] eq xalpair[0] and wavpair[1] eq xalpair[1]) then $ plotlabel='check xal: '+plotlabel ; time arrays for original and left-over samples sdt0=anytim2utc(timearr,/external) jdarr0=julday(sdt0.month,sdt0.day,sdt0.year,sdt0.hour,$ sdt0.minute,sdt0.second) sdt=anytim2utc(timearrleft,/external) jdarr=julday(sdt.month,sdt.day,sdt.year,sdt.hour,sdt.minute,sdt.second) ; define plot axes dummy=label_date(date_format='%H:%I') xtitle='time [UT]' ytitle='offset [px = 0.6 arcsec]' thick=1 timerange=[min(jdarr0),max(jdarr0)] timespan=timerange[1]-timerange[0] xrange=[timerange[0]-0.15*timespan,timerange[1]+0.15*timespan] ; shifts range may have terrible outliers so set scale for remaining datarange=[min([shiftx[indstay]-confx[indstay],$ shifty[indstay]-confy[indstay]]),$ max([shiftx[indstay]+confx[indstay],$ shifty[indstay]+confy[indstay]])] dataspan=datarange[1]-datarange[0] yrange=[datarange[0]-0.15*dataspan,datarange[1]+0.15*dataspan] ; start shifts plot (formerly "full_shifts") psfilename=driftsdir+'shifts_'+pairname+'.ps' if (aligncenter ne 0 and wavpair[0] eq xalpair[0] $ and wavpair[1] eq xalpair[1]) then $ psfilename=driftsdir+'check_xal_'+pairname+'.ps' openpsplot,psfilename,fontsize=8,thick=thick plot,jdarr0,splinx,$ thick=2,$ position=[0.15,0.15,0.95,0.9],/normal,$ xticklen=0.03,yticklen=0.03/1.6,$ xtickformat='label_date',xtickunits='Time',$ xtickinterval=fix((xrange[1]-xrange[0])/5.),xminor=2,$ xtitle=xtitle,ytitle=ytitle,$ xrange=xrange,xstyle=1,yrange=yrange,ystyle=1 oplot,jdarr0,spliny,thick=2 ; label curves xyouts,jdarr0[0]-0.08*timespan,splinx[0],charsize=1.2,'x' xyouts,jdarr0[0]-0.08*timespan,spliny[0],charsize=1.2,'y' xyouts,jdarr0[nt-1]+0.06*timespan,splinx[nt-1],charsize=1.2,'x' xyouts,jdarr0[nt-1]+0.06*timespan,spliny[nt-1],charsize=1.2,'y' ; overplot original data points oplot,jdarr0,shiftx,psym=2,symsize=4./sqrt(nt) oplot,jdarr0,shifty,psym=6,symsize=4./sqrt(nt) ; overplot 95% confidence bars from tile statistics per time step np=n_elements(jdarr0) for ip=0,np-1 do begin plots,[jdarr0[ip],jdarr0[ip]],[shiftx[ip]-confx[ip],shiftx[ip]+confx[ip]],$ noclip=0,color=128 plots,[jdarr0[ip],jdarr0[ip]],[shifty[ip]-confy[ip],shifty[ip]+confy[ip]],$ noclip=0,color=64,linestyle=1,thick=2 ; make distinct for overlap endfor ; mark removed outliers by adding box around them match,timearr,timearrleft,indstay xshiftgone=shiftx yshiftgone=shifty xshiftgone[indstay]=-9999 yshiftgone[indstay]=-9999 oplot,jdarr0,xshiftgone,psym=6,symsize=8./sqrt(nt) oplot,jdarr0,yshiftgone,psym=6,symsize=8./sqrt(nt) ; specify mean of tile confidence limits ; mean of grey bars = +- 95% confidence limit per it step obtained from tiles xyouts,0.18,0.18,/normal,$ 'cfx'+trimd(mean(confx),2)+' cfy'+trimd(mean(confy),2) ; add rms deviation from spline ; = standard deviation in offset from curve per sample xyouts,0.92,0.18,/normal,alignment=1,$ 'sdx'+trimd(sdevx,2)+' sdy'+trimd(sdevy,2) ; add boxcarlabel when nonempty if (boxcarlabel ne '') then xyouts,0.55,0.82,/normal,'box'+boxcarlabel ; add delaylabel when nonempty at center if (delaylabel ne '') then xyouts,0.70,0.82,/normal,delaylabel ; add smearlabel when nonempty (default for euvanchor pairs) if (smearlabel ne '') then xyouts,0.85,0.82,/normal,smearlabel ; add figtitle figtitle=plotlabel+timestring xyouts,0.5,0.95,figtitle,/normal,alignment=0.5 ; plot done closepsplot,psfilename,opengv=verbose ; ========= write spline structures as fits files ;RR I didn't understand how to read a structure with unknown dimensions ;RR and so I write out all in longhand = clumsy - but it works ; write spline coefficients on files (format correction Maria Loukitcheva) outfilex=driftsdir+$ 'splinx-'+ntostr(wavpair[0])+ntostr(wavpair[1])+'.dat' writecol,outfilex,timearrleft,ax,bx,cx,dx,fmt='(5E21.12)' outfiley=driftsdir+$ 'spliny-'+ntostr(wavpair[0])+ntostr(wavpair[1])+'.dat' writecol,outfiley,timearrleft,ay,by,cy,dy,fmt='(5E21.12)' SKIPTHISPAIR: ; clean up spawn,'rm -f /tmp/sdosmear*.fits' end ; =============== main for testing per IDLWAVE H-c ====================== ;; datadir='/home/rutten/data/SDO/2014-06-14-small/center/' ;; datadir='/home/rutten/data/SDO/2014-06-14-longer/center/' ;; datadir='~/data/SDO/2019-11-11-transit/complete_notrack_1700/center' ;; datadir='/home/rutten/data/SDO/2019-11-11-transit/midpointfull/target' ;; datadir='/media/rutten/RRHOME/alldata/SDO/2019-11-11-transit/full_long/center' datadir='/media/rutten/RRHOME/alldata/SolO/2020-05-30-first/sdo/center' ;RR select pair below ;; wavpair=[9,11] ; 1700 > mag ; exact ;; wavpair=[8,9] ; 1600 > 1700 ; exact although many differences ;; wavpair=[6,8] ; 304 > 1600 ; traditional anchor EV > UV ;; wavpair=[6,9] ; 304 > 1700 ; new anchor EUV > UV ;; wavpair=[6,11] ; 304 > mag ; newer euvanchor EUV > HMI ;; wavpair=[5,6] ; 211 > 304 ; easy for quiet areas ;; wavpair=[3,6] ; 171 > 304 ; difficult ;; wavpair=[1,5] ; 94 > 221 ; OK ;; wavpair=[7,5] ; 335 > 221 ; OK ;; wavpair=[4,5] ; 193 > 211 ; strange ;; wavpair=[2,1] ; 131 > 94 ; OK ;; wavpair=[3,5] ; 171 > 211 ;; wavpair=[5,7] ; 211 > 335 ;; wavpair=[7,1] ; 335 > 94 ;; wavpair=[1,6] ; 94 > 304 ; tel pair as check ;; wavpair=[1,4] ; 94 > 131 ;; wavpair=[1,3] ; 94 > 171 ; same tel pair as 304 > 1600 ;; wavpair=[4,3] ; 193 > 171 ;; wavpair=[2,5] ; 131 > 221 ; Pradeep test ;; wavpair=[2,3] ; 131 > 171 ; classic finalcheck ;; wavpair=[2,9] ; 131 > 1700 ; old finalcheck wavpair=[2,6] xalpair=wavpair ; OOPS: Y O.5 py, X -0.5 px aligncenter=1 delaymin=0 sdo_writepairspline,datadir,wavpair,delaymin=delaymin,$ /verbose,blink=0,aligncenter=aligncenter,xalpair=xalpair end