; file: stx_findlocation.pro = use SDO to find where STX was pointing ; init: May 5 2017 Rob Rutten Deil from sst_findlocation.pro ; last: Aug 12 2022 Rob Rutten Deil ; site: rridl/sdolib ;+ pro stx_findlocation,im_stx,datetime,$ x_stx,y_stx,angle_stx,px_stx,$ skipmetcalf=skipmetcalf,$ xyrange=xyrange,xystep=xystep,$ anglerange=anglerange,ccps=ccps,ccshiftps=ccshiftps,$ minangle=minangle,maxangle=maxangle, anglestep=anglestep,$ sdowav=sdowav,trimboxstx=trimboxstx,$ smearstx=smearstx,histopstx=histopstx,flatten=flatten,$ muckdarksdo=muckdarksdo,muckdarkstx=muckdarkstx,$ muckbrightsdo=muckbrightsdo,muckbrightstx=muckbrightstx,$ norinse=norinse,shiftfor=shiftfor,shiftrev=shiftrev,$ nxfor=nxfor,nyfor=nyfor,nxrev=nxrev,nyrev=nyrev,$ writesdoforward=writesdoforward,writesplinip=writesplinp,$ blink=blink,checkmetcalf=checkmetcalf,show=show,verbose=verbose ; find precise solar X,Y location and orientation of an STX image ; (STX = solar telescope X taking e.g., granulation images) ; - apply image mucking for difficult comparisons (examples at bottom) ; - rotate STX image to input (X,Y), resample to AIA pixels (0.6 arcsec) ; - find best-matching SDO cutout by rastered cross-correlation matching ; - anglerange mode (slow): find best angle_stx in trial range ; (last resort if you know only datetime and scale, bad (X,Y) and angle) ; - optional fine-tuning STX angle and pixel scale per Metcalf ; advised comparison pair: STX continuum versus HMI continuum ; alternatives: STX profile longest wavelength versus HMI continuum ; STX magnetogram versus HMI magnetogram ; ; EXAMPLES: ; at end of this file ; ; INPUTS: ; im_stx = STX image (2D array of any type) ; datetime = string in mandatory format '2014.06.19_16:03' ; x_stx, y_stx: initial estimate STX (X,Y) pointing (arcsec) ; angle_stx: initial FOV angle (deg NESW = counterclock from North) ; px_stx: initial STX pixel size (arcsec) ; ; OPTIONAL KEYWORD INPUTS: ; skipmetcalf: skip (slow) Metcalf improvement of STX image angle and scale ; xyrange = half X,Y search range in 0.6 arcsec AIA px (default 50) ; xystep = stepsize in AIA px (integer, default 1 = 0.6", maybe 3) ; anglerange = 1/0: run outer loop over range of trial angles (very slow) ; then: specify result plotfiles ccps, ccshiftps ; specify angle range with minangle, maxangle, anglestep (degrees) ; (for HMI granulation anglestep not above 10, default 1) ; sdowav = SDO diagnostic to correlate with (default 'continuum') ; trimboxstx: 4-elem int vector [xmin,ymin,xmax,ymax] to select a subfield ; this should be centered to maintain field-of-vie location ; use also to ignore non-straight borders (such as black or grey rims) ; smearstx: smear STX image (default -1.5 = 1.5*pxratio) ; histopstx: hist_opt STX image with given value (clip limit < 1, top only) ; flatten: subtract boxcar average from each, value = width in coarsest px ; use for spots or limb in the field ; muckdark sdo/ssx: set intensities below threshold (units mean=1) to mean ; muckbright sdo/ssx: set intensities above threshold to mean ; norinse = 1/0: no final rinse in findalignimages.pro ; writesdoforward = 1/0: write forward SDO result ; writesplinip = 1/0: splinip for forward SDO result write ; blink: blink duration for [SDO cutout - STX trimbox] (default 10 seconds) ; checkmetcalf: showex final Metcalf blink cutouts (may replace blink) ; show = 0/1/2/3: as findalignimages.pro (1 = start result showexes) ; verbose: 0 = no printout, 1 = slight, 2 = more ; ; OUTPUTS angle-range mode: ; psplots and printout specifying best angle_stx in trial range ; ; OUTPUTS single-angle mode: ; x_stx, y_stx: improved (X,Y) coordinates center of STX field of view ; angle_stx: improved STX orientation angle NESW (counterclockwise) ; px_stx: improved STX pixel size in arcsec ; ; OPTIONAL KEYWORD OUTPUTS ; nxfor,nyfor = central-cut common size in finest px ; shiftrev = im2 [xshift,yshift] in finest-px units for STX > SDO ; nxrev,nyrev = central-cut size in finest px, accommodating STX rotation ; RESTRICTIONS: ; needs image at good STX seeing and no post-eclipse SDO defocus ; WARNING: the stx_ parameters are sticky when rerunning within session ; ; HISTORY: ; May 5 2017 RR: start ; May 19 2017 RR: improved Metcalf per findalignimages.pro ; May 3 2020 RR: anglerange, interpolate straddle pair ; May 18 2020 RR: use sdo_getimage.pro ; Aug 1 2022 RR: checked for sst_getlocation_sstred.pro, still OK ;- ; answer no-parameter query if (n_params(0) lt 6) then begin sp,stx_findlocation return endif ; defaults for keywords if (n_elements(skipmetcalf) eq 0) then skipmetcalf=0 if (n_elements(xyrange) eq 0) then xyrange=50 ; NB default if (n_elements(xystep) eq 0) then xystep=1 if (n_elements(anglerange) eq 0) then anglerange=0 if (n_elements(ccps) eq 0) then ccps='/tmp/cc.ps' if (n_elements(ccshiftps) eq 0) then ccshiftps='/tmp/ccshift.ps' if (n_elements(minangle) eq 0) then minangle=-120 if (n_elements(maxangle) eq 0) then maxangle=+120 if (n_elements(anglestep) eq 0) then anglestep=1 if (n_elements(sdowav) eq 0) then sdowav='continuum' if (n_elements(trimboxstx) eq 0) then trimboxstx=-1 if (n_elements(smearstx) eq 0) then smearstx=-1.5 ; NB default if (n_elements(histopstx) eq 0) then histopstx=0 if (n_elements(muckdarksdo) eq 0) then muckdarksdo=0 if (n_elements(muckdarkstx) eq 0) then muckdarkstx=0 if (n_elements(muckbrightsdo) eq 0) then muckbrightsdo=0 if (n_elements(muckbrightstx) eq 0) then muckbrightstx=0 if (n_elements(norinse) eq 0) then norinse=0 ; changed default if (n_elements(flatten) eq 0) then flatten=0 ; changed default if (n_elements(shiftfor) eq 0) then shiftfor=0 if (n_elements(shiftrev) eq 0) then shiftrev=0 if (n_elements(nxfor) eq 0) then nxfor=0 if (n_elements(nxrev) eq 0) then nxrev=0 if (n_elements(nyfor) eq 0) then nyfor=0 if (n_elements(nyrev) eq 0) then nyrev=0 if (n_elements(writesdoforward) eq 0) then writesdoforward=0 if (n_elements(writesplinip) eq 0) then writesplinip=0 if (n_elements(blink) eq 0) then blink=10 if (n_elements(checkmetcalf) eq 0) then checkmetcalf=0 if (n_elements(show) eq 0) then show=0 if (n_elements(verbose) eq 0) then verbose=0 ; check format of datetime if (strmid(datetime,4,1) ne '.' or $ strmid(datetime,7,1) ne '.' or $ strmid(datetime,10,1) ne '_' or $ strmid(datetime,13,1) ne ':') then begin print, ' ##### stx_findlocation abort: datetime not YYYY.MM.DD_HH:MM:SS' return endif ; more checks if (fix(xystep) ne xystep) then begin print,' ##### stx_findlocation abort: xystep not integer' return end ; pixel sizes sdo_px=0.6 ; AIA px size holds also for HMI after aia_prep below ; ====== treat STX image ; NB: it may have flat borders along its x and y axes but should have ; no derotation triangles ; histop, muck, smear as requested and resample to SDO px if (smearstx lt 0) then smearstx=abs(smearstx)*sdo_px/px_stx reformimage,im_stx,stxmuck,$ congridfactor=[px_stx/sdo_px,px_stx/sdo_px],$ histopttop=histopstx,muckdark=muckdarkstx,muckbright=muckbrightstx,$ smear=smearstx,flatten=flatten,missingvalue=-1,splinip=1 ; make dimensions even (since nx/2, ny/2 below) sizeim=size(stxmuck) nx=fix(sizeim[1]/2.)*2 ny=fix(sizeim[2]/2.)*2 stxmuck=stxmuck[0:nx-1,0:ny-1] ;; ; check ;; sv,im_stx & sv,stxmuck ;; STOP ;; wdelall ; ====== get SDO image (copied from sdo_featurelocator.pro) ; define permitted SDO (AIA + HMI) level1 fitsfilename strings wavs=['94','131','171','193','211','304','335','1600','1700',$ ; AIA 'magnetogram','continuum','Dopplergram'] ; HMI nwavs=n_elements(wavs) hmiseries=['hmi.M_45s','hmi.Ic_45s','hmi.V_45s'] ; check sdowav specification iwav=-1 for iw=0,nwavs-1 do if (wavs[iw] eq sdowav) then iwav=iw if (iwav eq -1 ) then begin print,' ===== ERROR: sdowav invalid = ',sdowav print,' ===== valid: ',wavs return endif ; set aia_prep file label lev2wav=sdowav if (iwav eq 9) then lev2wav='mag' if (iwav eq 10) then lev2wav='cont' if (iwav eq 11) then lev2wav='dop' ; get old image written by sdo_getimage or new straddle-interpolated one oldsdoimage='/tmp/sdo_'+datetime+'_'+wavs[iwav]+'.fits' if (file_test(oldsdoimage) ne 0) then begin if (verbose) then print,' ----- found '+oldsdoimage+' and use that' read_sdo,oldsdoimage,indexim,im if (n_elements(im) eq 0) then begin print,' ##### stx_findlocation abort: wrong old image in /tmp' return endif endif else sdo_getimage,datetime,wavs[iwav],im,indexim,/intpolpair ; set standard clips clipmin=-1 clipmax=-1 ; magnetogram: clip to see weaker fields if (iwav eq 9) then $ im=sdo_intscale(indexim,im,lev2wav,clipmin=-1000,clipmax=+1000) ; pretty up if (iwav lt 9) then im=histo_opt_rr(im) if (iwav eq 7 or iwav eq 8) then im=sqrt(im>1) sdoim=im xcen=2048 ycen=2048 ; start off for unknown angle_stx: extra loop over range of trial angles if (anglerange eq 1) then begin nangles=fix((maxangle-minangle)/anglestep+0.5)+1 angles=minangle+indgen(nangles)*anglestep maxccangle=fltarr(nangles) maxccshiftangle=fltarr(nangles) xstxangle=fltarr(nangles) ystxangle=fltarr(nangles) endif else begin nangles=1 endelse ; ===== start loop over range angle trials for iangle=0,nangles-1 do begin if (anglerange eq 1) then angle_stx=angles[iangle] ; rotate SDO to STX orientation (positive = CCW) ; (in this way avoid corner triangles that you get the other way) reformimage,sdoim,sdoim2,rotate=angle_stx,splinip=1 ;; ; check ;; sv,sdoim & sv,sdoim2 ;; STOP ; ===== double loop over X and Y grid position ; setup X and Y grid theta=angle_stx*!pi/180. xrot0=x_stx*cos(theta)-y_stx*sin(theta) yrot0=x_stx*sin(theta)+y_stx*cos(theta) xpx0=fix(xrot0/sdo_px) ; AIA px also for HMI after aia_prep ypx0=fix(yrot0/sdo_px) xshiftim=intarr(4096,4096) yshiftim=intarr(4096,4096) ccmaxim=fltarr(4096,4096) ccshiftim=fltarr(4096,4096) ; start double loop over X and Y grid for ix0=xpx0-xyrange,xpx0+xyrange,xystep do begin for iy0=ypx0-xyrange,ypx0+xyrange,xystep do begin ; check that this location is still on disk if (sqrt(ix0^2+iy0^2) gt 2048.) then goto, OFFDISK ; get rotated SDO image cutout with size of coarse STX image selsdo=sdoim2[ix0+xcen-nx/2:ix0+xcen+nx/2-1,$ iy0+ycen-ny/2:iy0+ycen+ny/2-1] ; muck this selected SDO cutout as requested reformimage,selsdo,sdomuck,flatten=flatten,$ muckdark=muckdarksdo,muckbright=muckbrightsdo,splinip=1 ; find cross-correlation peak value and location thisshift=findimshift_rr(sdomuck,stxmuck,niter=0,$ maxcc=maxcc,$ croptriangles=1,cropborders=1) if (verbose ge 2) then $ print,' -----ix0,iy0 ='+trimd([ix0,iy0])+$ ' maxcc ='+trimd(maxcc)+' thisshift ='+trimd(thisshift) ; fill result images xshiftim[ix0+xcen,iy0+ycen]=thisshift[0] yshiftim[ix0+xcen,iy0+ycen]=thisshift[1] ccmaxim[ix0+xcen,iy0+ycen]=maxcc ccshiftim[ix0+xcen,iy0+ycen]=$ maxcc/(sqrt(thisshift[0]^2+thisshift[1]^2)+1) OFFDISK: endfor ; end loop over ix0 endfor ; loop over iy0 if (max(ccshiftim) eq 0) then begin print,' ##### stx_findlocation abort: no peak in ccshiftim' STOP endif ;; ; inspect (hit "quit this" to continue) ;; showex,ccmaxim,ccshiftim ; zoom in to center: ccshiftim peak = 1 px if (anglerange eq 1) then begin maxccangle[iangle]=max(ccmaxim) maxccshiftangle[iangle]=max(ccshiftim) ; normalized by shift length iangmax=where(max(ccmaxim) eq ccmaxim) ; get x_stx and y_stx estimates for this angle (as below) cctop=max(ccshiftim,maxlocation) xpxmax=maxlocation mod 4096 ; convert 1D to 2D index value ypxmax=maxlocation/4096 xshift=xshiftim[maxlocation] yshift=yshiftim[maxlocation] stx_rotpx=xpxmax+xshift stx_rotpy=ypxmax+yshift stx_rotx=(stx_rotpx-xcen)*sdo_px stx_roty=(stx_rotpy-ycen)*sdo_px xstxangle[iangle]=stx_rotx*cos(theta)+stx_roty*sin(theta) ystxangle[iangle]=-stx_rotx*sin(theta)+stx_roty*cos(theta) ; define marker flag for growth/reduction flag='' if (iangle gt 0) then begin if (maxccangle[iangle] gt 1.2*previous) then flag=' ++' if (maxccangle[iangle] gt 1.5*previous) then flag=' ++++' if (maxccangle[iangle] lt 0.67*previous) then flag=' ----' if (maxccangle[iangle] lt 0.83*previous) then flag=' --' endif previous=maxccangle[iangle] ; print results for this angle if (verbose ne 0) then print,' ----- iangle ='+trimd(iangle)+$ ' angle_stx ='+trimd(angle_stx)+$ ' max cc ='+trimd(maxccangle[iangle])+$ ' max cc/shift ='+trimd(maxccshiftangle[iangle])+$ ' max x_stx ='+trimd(xstxangle[iangle])+$ ' max y_stx ='+trimd(ystxangle[iangle])+flag endif endfor ; end loop over angles when trial range ; if anglerange print and plot final result and quit if (anglerange eq 1) then begin ; plot and print maxcc ianglemaxcc=where(maxccangle eq max(maxccangle)) openpsplot,ccps,thick=2,fontsize=9,xsize=9,ysize=6 plot,angles,maxccangle,$ xtitle='angle_stx',ytitle='max cc',psym=-6,symsize=1.3,/ynozero xyouts,0.2,0.95,/norm,'max cc at angle = '+trimd(angles[ianglemaxcc]) closepsplot,ccps,opengv=1 print,' ----- maxcc reached at iangle ='+trimd(ianglemaxcc)+$ ' angle_stx ='+trimd(angles[ianglemaxcc])+$ ' x_stx ='+trimd(xstxangle[ianglemaxcc])+$ ' y_stx ='+trimd(ystxangle[ianglemaxcc]) ; plot and print ccshift (narrower peak, needs finer sampling) ianglemaxccshift=where(maxccshiftangle eq max(maxccshiftangle)) openpsplot,ccshiftps,thick=2,fontsize=9,xsize=9,ysize=6 plot,angles,maxccshiftangle,$ xtitle='angle_stx',ytitle='max cc/shift',psym=-6,symsize=1.3,/ynozero xyouts,0.2,0.95,/norm,'max ss/shift at angle = '+$ trimd(angles[ianglemaxccshift]) closepsplot,ccshiftps,opengv=1 print,' ----- maxccshift reached at iangle ='+trimd(ianglemaxccshift)+$ ' angle_stx ='+trimd(angles[ianglemaxccshift])+$ ' x_stx ='+trimd(xstxangle[ianglemaxccshift])+$ ' y_stx ='+trimd(ystxangle[ianglemaxccshift]) ; quit return endif ; =========== rest program is for call with single angle that should be close ; get max = max cc where shift minimal (cc divided by shift length) cctop=max(ccshiftim,maxlocation) xpxmax=maxlocation mod 4096 ; convert 1D to 2D index value ypxmax=maxlocation/4096 xshift=xshiftim[maxlocation] yshift=yshiftim[maxlocation] stx_rotpx=xpxmax+xshift stx_rotpy=ypxmax+yshift stx_rotx=(stx_rotpx-xcen)*sdo_px stx_roty=(stx_rotpy-ycen)*sdo_px x_stx=stx_rotx*cos(theta)+stx_roty*sin(theta) y_stx=-stx_rotx*sin(theta)+stx_roty*cos(theta) if (verbose ne 0) then $ print,' ----- raster cc max: (x,y) px ='+trimd([stx_rotpx,stx_rotpy])+$ ' shift ='+trimd([xshift,yshift])+$ ' (X,Y) ='+trimd([x_stx,y_stx],1) ;;; check cc images ;; sv,ccmaxim[xpx0+xcen-xyrange:xpx0+xcen+xyrange,$ ;; ypx0+ycen-xyrange:ypx0+ycen+xyrange]>.5*cctop ;; sv,ccmaxim[stx_rotpx-2*xyrange:stx_rotpx+2*xyrange,$ ;; stx_rotpy-2*xyrange:stx_rotpy+2*xyrange]>.5*cctop ;; sv,ccshiftim[stx_rotpx-2*xyrange:stx_rotpx+2*xyrange,$ ;; stx_rotpy-2*xyrange:stx_rotpy+2*xyrange]>.5*cctop ; get finally selected SDO cutout finsdo=sdoim2[stx_rotpx-nx/2:stx_rotpx+nx/2-1,$ stx_rotpy-ny/2:stx_rotpy+ny/2-1] ; blink muck result unless Metcalf blinks follow if (verbose ne 0 and blink ne 0 and skipmetcalf ne 0) then begin reformimage,finsdo,finsdomuck,flatten=flatten,$ muckdark=muckdarksdo,muckbright=muckbrightsdo,splinip=1 print,' ----- blink: selected rotated-SDO patch versus resampled STX' flickrr,finsdomuck,stxmuck,duration=blink endif ; ===== angle and pixel improvement per Metcalf if (skipmetcalf eq 0) then begin ; now work on intensity-scaled nonrotated SDO cutout sdo_cut=sdoim[x_stx/sdo_px+xcen-nx:x_stx/sdo_px+xcen+nx,$ y_stx/sdo_px+ycen-ny:y_stx/sdo_px+ycen+ny] ;; ; check or write as test image for e.g. findalignimages.pro ;; sv,sdo_cut ;; STOP ;; ;; wdelall ; get align parameters per Metcalf findalignimages,sdo_cut,im_stx,sdo_px,px_stx,angle_stx,$ px2asym,shiftfor,shiftrev,nxfor,nyfor,nxrev,nyrev,$ trimboxim2=trimboxstx,$ smearim2=smearstx,histopim2=histopstx,flatten=flatten,$ muckdarkim1=muckdarksdo,muckdarkim2=muckdarkstx,$ muckbrightim1=muckbrightsdo,muckbrightim2=muckbrightstx,$ norinse=norinse,$ blink=blink,checkmetcalf=checkmetcalf,show=show,verbose=verbose ; apply shifts x_stx=x_stx+shiftfor[0]*px_stx y_stx=y_stx+shiftfor[1]*px_stx ; write forward sdo-transformed image if (writesdoforward) then begin reformimage,sdo_cut,sdo_for,congridfactor=1./[px_stx/0.6,px_stx/0.6],$ shift=shiftfor,rotate=angle_stx,cutcentralx=nxfor,cutcentraly=nyfor,$ splinip=writesplinip writefits,'sdo2stx_'+sdowav+'_forward.fits',sdo_for endif endif ; end of optional angle and pixel improvememt ; print final result if (verbose ne 0) then $ print,' ===== result: (X,Y) ='+trimd([x_stx,y_stx],1)+$ ' angle_stx ='+trimd(angle_stx,2)+$ ' px_stx ='+trimd(px_stx,4) end ; =============== main for testing per IDLWAVE H-c ====================== ; ------------------------------------------------old standard test ;RR rather difficult sunspot with dark umbra and many bright points in moat ;RR also in sst_findlocation.pro tests calling this program cd,'/media/rutten/RRHOME/alldata/SST/2016-09-05-demo' im_stx=readfitscube('sst/ha_wb.fits',trange=[18,18]) datetime='2016.09.05_09:54:34.626' histopstx=0.1 ; muck moat BPs, tops get cut and then smeared xyrange=30 x_stx=-117.2 y_stx=18.6 angle_stx=98.22 px_stx=0.0568 flatten=5 blink=5 verbose=1 show=1 skipmetcalf=0 writesdoforward=1 writesplinip=1 stx_findlocation,im_stx,datetime,$ x_stx,y_stx,angle_stx,px_stx,$ skipmetcalf=skipmetcalf,$ xyrange=xyrange,xystep=xystep,$ anglerange=anglerange,ccps=ccps,ccshiftps=ccshiftps,$ minangle=minangle,maxangle=maxangle,anglestep=anglestep,$ sdowav=sdowav,trimboxstx=trimboxstx,$ smearstx=smearstx,histopstx=histopstx,$ muckdarksdo=muckdarksdo,muckdarkstx=muckdarkstx,$ muckbrightsdo=muckbrightsdo,muckbrightstx=muckbrightstx,$ norinse=norinse, shiftfor=shiftfor,shiftrev=shiftrev,$ nxfor=nxfor,nyfor=nyfor,nxrev=nxrev,nyrev=nyrev,$ flatten=flatten,writesdoforward=writesdoforward,$ writesplinip=writesplinp,$ blink=blink,verbose=verbose,show=show end