; file: sdo_makereffitscube.pro ; init: Dec 29 2013 Rob Rutten Deil ; last: Jul 6 2020 Rob Rutten Deil ;+ pro sdo_makereffitscube,datadir,$ clipmin=clipmin,clipmax=clipmax,no_intscale=no_intscale,$ refwav=refwav,border=border,itend=itend,verbose=verbose ; PURPOSE: ; convert level2 ("1.5") SDO images at reference wavelength into fitscube ; with timing info in header. It only crops as set by border parameter. ; If verbose outputs shifts needed to undo the 1-px rotation sawtooth ; but its application (if there notrack=1) and crossalign (nocrossalign=0) is ; done in sdo_makeslavefitscube.pro also for the reference sequence. ; ; CALL EXAMPLE: ; see end of this program file ; ; INPUTS: ; datadir = string path/dir containing dir level2 ; ; OPTIONAL KEYWORD INPUTS: ; border = border with to cut x and y borders for unequal cutouts (def 15) ; refwav: string specifying AIA or HMI reference wavelength, default '171' ; clipmin, clipmax: intensity clipping in sdo_intscale.pro ; no_intscale = 1: do not intscale target AIA and HMIcont intensities ; itend=value: earlier stop for smaller-volume tests ([0] = all = default) ; verbose=1/[0]: print and plot image shifts ; ; OUTPUT: ; fitscube file still containing the JSOC 1-px sawtooth shifts and drifts ; header keywords XCEN and YCEN are from HMI when available ; ; METHOD: ; uses assoc to permit larger cubes than memory size ; ; HISTORY: ; Dec 29 2013 RR: start ; Jan 31 2014 RR: crop grey edges ; Feb 25 2016 RR: aiahmi out; also EUV > HMI ; Apr 9 2016 RR: no shifting (left for sdo_makeslavefitscube.pro) ; Mar 30 2017 RR: clipmin, clipmax ; Jan 10 2018 RR: option no_intscale ; Sep 3 2018 RR: set XCEN, YCEN keywords to the HMI values ;- ; answer no-parameter query if (n_params() lt 1) then begin sp,sdo_makereffitscube return endif ; defaults for keywords if (n_elements(refwav) eq 0) then refwav='171' if (n_elements(clipmin) eq 0) then clipmin=0 if (n_elements(clipmax) eq 0) then clipmax=0 if (n_elements(no_intscale) eq 0) then no_intscale=0 if (n_elements(border) eq 0) then border=15 if (n_elements(itend) eq 0) then itend=0 if (n_elements(verbose) eq 0) then verbose=0 ; define datadir subdirs level2dir=datadir+'/level2' cubesdir=datadir+'/cubes' ; make dirs if not yet existing spawn,'mkdir -p '+cubesdir if (verbose ne 0) then spawn,'mkdir -p '+cubesdir+'/plots_align/' ; define permitted level2 SDO (AIA + HMI) filename strings wavs=['94','131','171','193','211','304','335','1600','1700',$ 'mag','cont','dop'] nwavs=n_elements(wavs) fitscubenames=['aia94.fits','aia131.fits','aia171.fits',$ 'aia193.fits','aia211.fits','aia304.fits',$ 'aia335.fits','aia1600.fits','aia1700.fits',$ 'hmimag.fits','hmicont.fits','hmidop.fits'] ; set channel names for plots channels=wavs for i=0,8 do channels[i]='AIA '+wavs[i] channels[9]='HMI magnetogram' channels[10]='HMI continuum' channels[11]='HMI Dopplergram' ; check wav specification; get iwav iwav=-1 wav=refwav for iw=0,nwavs-1 do if (wavs[iw] eq wav) then iwav=iw if (iwav eq -1 ) then begin print,' ##### FATAL ERROR: refwav invalid name = ',wav return endif ; get refwav file list (undo compression if fits files fpacked into .fz) ;RR wav=refwav here spawn,'funpack -D '+datadir+'/level2/*'+wav+'.fits.fz' files=findfile(datadir+'/level2/*'+wav+'.fits*') nfiles=n_elements(files) if (nfiles eq 1) then message,' ##### FATAL ERROR: no files for refwav' if (itend eq 0) then itend=nfiles-1 ; get header info first refwav image read_sdo, files[0],index0,image0,/nodata nx0=index0.naxis1 ny0=index0.naxis2 cenxpx0=nx0/2. cenypx0=ny0/2. t0=index0.t_obs ; becomes starttim in output header ; determine refwav cadence, goes into header ; and becomes the timing for all cross-aligned SDO files ; average of max 100 files (not nfiles since segmentation may have jumps) itcad=min([100,nfiles-1]) ; segmentation occurs only for more than 100 read_sdo,files[itcad],indexcad,/nodata tcad=indexcad.t_obs cadence=(anytim2tai(tcad)-anytim2tai(t0))/float(itcad) ; get XCEN, YCEN = image center in (X,Y) arcsec from solar center ; Sep 3 2018 Deil: take the HMI values if present since the whole ; pipeline shifts to HMI while HMI measures them properly for ; each image from the limb fit done for every HMI image spawn,'funpack -D '+datadir+'/level2/*cont.fits.fz' filescont=findfile(datadir+'/level2/*cont.fits') ncont=n_elements(filescont) ; 1 if no file present spawn,'funpack -D '+datadir+'/level2/*mag.fits.fz' filesmag=findfile(datadir+'/level2/*mag.fits') nmag=n_elements(filesmag) ; HMI cont present or HMI mag present if (ncont gt 1) then begin read_sdo, filescont[0],indexcont0,imagecont0,/nodata xcenarc0=indexcont0.xcen ycenarc0=indexcont0.ycen endif else if (nmag gt 1) then begin read_sdo, filesmag[0],indexmag0,imagemag0,/nodata xcenarc0=indexmag0.xcen ycenarc0=indexmag0.ycen endif ; no HMI cont or mag: instead take the first refwav XCEN, YCEN if (ncont lt 2 and nmag lt 2) then begin xcenarc0=index0.xcen ycenarc0=index0.ycen endif ;; ; check ;; print,' ####### xcenarc0 ='+trimd(xcenarc0) ;; STOP ; define cropped image center in (X,Y) = arcsec from solar center xcenarc=xcenarc0 ; remains the same since I use equal border cropping ycenarc=ycenarc0 ; define outheader with timing and location parameters for future use bigendian=1 ;RR my Toshiba laptop bitpix=16l ;RR output = integer cube mkhdr,outheader,abs(bitpix)/8,$ [nx0-2*border,ny0-2*border,itend+1] sxaddpar,outheader,'starttim',t0,'first image mid-exposure' sxaddpar,outheader,'cadence',cadence,'image cadence [s]' sxaddpar,outheader,'cdelt1',index0.cdelt1,'x-scale [arcsec/px]' sxaddpar,outheader,'cdelt2',index0.cdelt2,'y-scale [arcsec/px]' ;; sxaddpar,outheader,'crpix1',index0.crpix1-border,$ ;??? -border? ;; 'disk-center offset in px from leftmost of FOV = nr 1' ;; sxaddpar,outheader,'crpix2',index0.crpix2-border,$ ;??? -border? ;; 'disk-center offset in px from bottom of FOV = nr 1' ;RR Jul 2 2020 do not copy since image size changes so wrong sxaddpar,outheader,'rsun_obs',index0.rsun_obs,'solar radius [arcsec]' sxaddpar,outheader,'xcen',xcenarc,$ 'FOV center X [arcsec from sun center]' sxaddpar,outheader,'ycen',ycenarc,$ 'FOV center Y [arcsec from sun center]' sxaddpar,outheader,'channel',channels[iwav],$ 'SDO diagnostic name [string]' ; get size outheader sizeoutheader=size(outheader) outheadersize=(1+fix(sizeoutheader[1]/36.))*2880 ;RR fits header = Nx36 "card images" = Nx2880 bytes ; open output cube file for assoc, write header get_lun, unit_out if (bigendian) then openw,unit_out,cubesdir+'/'+fitscubenames[iwav],$ /swap_if_little_endian else openw,unit_out,$ cubesdir+'/'+fitscubenames[iwav] outassoc=assoc(unit_out,$ intarr(nx0-2*border,ny0-2*border),$ outheadersize) rec=assoc(unit_out, bytarr(outheadersize)) rec[0]=byte(outheader) ; start loop over timesteps ;RR Apr 9 2016 no image changes, get done in sdo_makeslavefitscube shiftx=fltarr(itend+1) shifty=fltarr(itend+1) for it=0,itend do begin ; read next image read_sdo,files[it],index,image,/uncomp_delete ; rescale (is this needed? gets redone in sdo_makeslavefitscube) im=sdo_intscale(index,image,wav,clipmin=clipmin,clipmax=clipmax,$ no_intscale=no_intscale) ; no shifting done here, only get the shifts for \verbose printoutput ;RR actual shifting in sdo_makeslavefitscube run also for reference if (verbose ne 0) then begin ;; ; undo the whole-pixel tracking roundoff ;; if (it eq 0) then missing=avg(im) else begin if (it ne 0) then begin newpos=rot_xy(xcenarc0,ycenarc0,tstart=t0,tend=index.t_obs) cenxpxit=newpos(0)/index.cdelt1+index.crpix1-0.5 cenypxit=newpos(1)/index.cdelt2+index.crpix2-0.5 shiftx[it]=cenxpxit-cenxpx0 shifty[it]=cenypxit-cenypx0 ;; ; check ;; imshift=shift_img(im,-[shiftx[it],shifty[it]],missing=missing) ;; showex,im,imshift ;; STOP endif endif ; optionally print shifts per it if (verbose ne 0) then print,' --- it = ',ntostr(it),$ ' index.crpix1,2 ='+trimd(index.crpix1)+trimd(index.crpix2)+$ ' shifts x,y =', ntostr([shiftx[it],shifty[it]],format='(2f10.3)') ; crop the image im=im[border:nx0-1-border,border:ny0-1-border] ; write nonshifted but scaled and cropped image (crop size needed) outassoc[it]=im endfor ; free the output file free_lun,unit_out ; optionally plot shifts if (verbose ne 0) then begin psfilename=cubesdir+'/plots_align/shifts-'+fitscubenames[iwav] psfilename=str_replace(psfilename,'.fits','.eps') openpsplot,psfilename,thick=2,fontsize=9,xsize=8.8,ysize=8.8 yrange=[min([min(shiftx),min(shifty)]),$ max([max(shiftx),max(shifty)])] yrange=[yrange[0]-0.1*(yrange[1]-yrange[0]),$ yrange[1]+0.1*(yrange[1]-yrange[0])] plot,indgen(itend+1),shiftx,$ xrange=[0,itend],yrange=yrange,xstyle=1,ystyle=1,$ xtitle='time step',ytitle='shift in px',$ position=[0.2,0.1,0.95,0.90] oplot,indgen(itend+1),shifty,linestyle=2 itlabel=(itend+1)/2 labeloffset=(yrange[1]-yrange[0])*0.02 xyouts,itlabel,shiftx[itlabel]+labeloffset,'x' xyouts,itlabel,shifty[itlabel]+labeloffset,'y' xyouts,0.2,0.92,/norm,'shifts '+fitscubenames[iwav],charsize=1.2 closepsplot,psfilename,opengv=verbose ; print success if (verbose ne 0) then print,$ ' ===== wrote ',cubesdir+'/'+fitscubenames[iwav] endif end ; =============== main for testing per IDLWAVE H-c ========================= datadir='/home/rutten/data/SDO/2014-06-14-small/center/' ;; datadir='/home/rutten/data/SDO/2019-11-11-transit/midpointfull/target' sdo_makereffitscube,datadir,/verbose,/no_intscale ;; ,itend=20 ; inspect showex,datadir+'/cubes/aia171.fits' end