; file: sdo_sst_writeblinkfile.pro ; init: Jun 7 2017 Rob Rutten Deil ; last: Jun 17 2017 Rob Rutten Deil ; todo: SST fcubes > make icubes maintaining zero, treat ns = Stokes ; note: Aug 12 2022 this seems old and no longer used ;+ pro sdo_sst_writeblinkfile,sdofiles,sstfiles,$ outfilestart=outfilestart,spfile=spfile,$ sdodir=sdodir,sstdir=sstdir,$ sstspecsavs=sstspecsavs,sstlcs=sstlcs,bytscale=bytscale,$ xrange=xrange,yrange=yrange,trange=trange,$ wavlabels=wavlabels,verbose=verbose ; NAME: ; sdo_sst_writeblinkfile.pro ; PURPOSE: ; assemble aligned SDO and SST files into a single "La Palma" cubefile ; CALL: ; see above ; INPUTS: ; sdofiles = string array desired SDO files or 'all' ; sstfiles = string array desired crispex files ; OPTIONAL KEYWORD INPUTS: ; outfilestart: output path + filename start (default sdodir+'sdo_sst') ; spfile = 1/0: also write CRISPEX .sp. file (can be very slow...) ; sdodir: path of dir with SDO input files (default 'sdo2sst/') ; sstdir: path of dir with SST input files (default 'sst/') ; sstspecsavs: string array list spectfile.WWWW.sav files with spect_pos ; sstlcs: wavelength indices of line center (default each = nwavs/2) ; bytscale=0: maintain intensities, outfilestart=icube (default) ; ## use for analysis with crispex.pro and figure production ; ## ? TODO: SST fcubes converted to icubes maintaining zero ; bytscale=1: bytscale each (sub-)image separately ; ## optimum greyscales for exploration with showex.pro ; bytscale=2: bytscale to min and max of first full image per wav ; ## fast fixed greyscales for showex.pro but maybe clipping ; bytscale=3: bytscale to min max of full cube per wav ; ## the correct way, but slow and outlier spikes may spoil ; bytscale=4: bytscale to min max of range-selected subcube ; ## dangerous for small subcubes, they won't intercompare ; xrange, yrange, trange: partial cube cutout (2-element index arrays) ; wavlabels=1/0: add wavelength labels in images ; verbose=1/0: print progress ; OUTPUTS: ; file outfile, name has xyt, nxy, nwavs specification ; RESTRICTIONS: ; output in SST format = "La Palma" little_endian, not fits file ; HISTORY: ; Jun 11 2017 RR: start (but never used myself) ;- ; answer no-parameter query if (n_params() lt 2) then begin print,' sdo_sst_writeblinkfile,sdofiles,sstfiles,' print,' [outfilestart=outfilestart,spfile=spfile,' print,' sdodir=sdodir,sstdir=sstdir,' print,' sstspecsavs=sstspecsavs,sstlcs=sstlcs,bytscale=bytscale,' print,' xrange=xrange,yrange=yrange,trange=trange,' print,' wavlabels=wavlabels,verbose=verbose]' sp,sdo_sst_writeblinkfile return endif ; defaults for keywords if (n_elements(sdodir) eq 0) then sdodir='sdo2sst/' if (n_elements(sstdir) eq 0) then sstdir='sst/' if (n_elements(outfilestart) eq 0) then outfilestart=sdodir+'/sdo_sst' if (n_elements(spfile) eq 0) then spfile=0 if (n_elements(sstspecsavs) eq 0) then sstspecsavs=0 if (n_elements(sstlcs) eq 0) then sstlcs=0 if (n_elements(bytscale) eq 0) then bytscale=0 if (n_elements(wavlabels) eq 0) then wavlabels=0 if (n_elements(xrange) eq 0) then xrange=[0,-1] if (n_elements(yrange) eq 0) then yrange=[0,-1] if (n_elements(trange) eq 0) then trange=[0,-1] if (n_elements(verbose) eq 0) then verbose=0 ; set wall-clock timer (seconds) timestart=systime(1) ; print opener if (verbose) then print,' ===== sdo_sst_writeblinkfile.pro starts' ; -------- define SDO files and wav lists ; define permitted SDO channel strings in mice display order allsdowavs=['cont','mag','dop','1700','1600','304','171','193',$ '335','211','131','94'] nallsdowavs=n_elements(allsdowavs) ; when sdofiles='all' select all SDO files that are present if (sdofiles[0] eq 'all') then begin sdowavs=strarr(nallsdowavs) ifile=-1 sdofiles=strarr(1) for iwav=0,nallsdowavs-1 do begin infile=findfile(sdodir+'/*'+allsdowavs[iwav]+'*.fits') if (infile ne '') then begin ifile=ifile+1 sdowavs[ifile]=allsdowavs[iwav] if (ifile eq 0) then sdofiles[0]=infile else sdofiles=[sdofiles,infile] endif endfor nsdofiles=n_elements(sdofiles) sdowavs=sdowavs[0:nsdofiles-1] ; verbose: print wavs of found files if (verbose) then begin print,' ===== '+ntostr(nsdofiles)+' SDO channels found: ' print,' '+ntostr(sdowavs) endif endif else begin ; otherwise use the specified list of SDO files sdofiles=sdodir+sdofiles nsdofiles=n_elements(sdofiles) sdowavs=strarr(nsdofiles) selwav=0 for ifile=0,nsdofiles-1 do begin for iwav=0,nallsdowavs-1 do begin if (strmatch(sdofiles[ifile],'*'+allsdowavs[iwav]+'*.fits')) then $ sdowavs[ifile]=allsdowavs[iwav] endfor endfor endelse ; get SDO geometry and datatype from the first SDO fits header header_sdo=headfits_rr(sdofiles[0]) header_sdosize=(1+fix(n_elements(header_sdo)/36.))*2880 bitpix_sdo=fxpar(header_sdo,'bitpix') nx_sdo=fix(fxpar(header_sdo,'naxis1')) ny_sdo=fix(fxpar(header_sdo,'naxis2')) nt_sdo=fix(fxpar(header_sdo,'naxis3')) ; check if (nx_sdo eq 0 or ny_sdo eq 0 or nt_sdo eq 0) then begin print,' ###### sdo_sst_writeblinkfile abort: '+$ ' no nx, ny, nt in header '+infile return endif ; set SDO wavs nwavsfile=intarr(nsdofiles)+1 ; SDO files have nwavs=1 labels=sdowavs ; ------- get dimensions etc for SST files nsstfiles=n_elements(sstfiles) nw_sst=intarr(nsstfiles) bitpix_sst=intarr(nsstfiles) for ifilesst=0,nsstfiles-1 do begin sstfile=sstdir+'/'+sstfiles[ifilesst] ; get header lp_readheader,sstfile,dims=ndim_sst,header=header_sst,datatype=data_sst,$ nx=nx_sst,ny=ny_sst,nt=nz_sst,endian=endian_sst ; split nz_sst into nw_sst and nt_sst nt_sst=nt_sdo nw_sst[ifilesst]=nz_sst/nt_sst nwavsfile=[nwavsfile,nw_sst[ifilesst]] ; add this SST set ; set bitpix_sst if (data_sst eq 1) then bitpix_sst[ifilesst]=8 if (data_sst eq 2) then bitpix_sst[ifilesst]=16 if (data_sst eq 4) then bitpix_sst[ifilesst]=-32 ; check file validity if (ndim_sst ne 3 or nx_sst eq 0 or ny_sst eq 0 or nz_sst eq 0) then begin print,' ##### sdo_sst_writeblinkfile abort: wrong SST file '+sstfile return sp,sdo_sst_writeblinkfile endif ; check quality to SDO image size if (nx_sst ne nx_sdo or ny_sst ne ny_sdo) then begin print,' ##### sdo_sst_writeblinkfile abort: '+$ 'SST file '+sstfile+' differs in nx,ny from SDO files' return sp,sdo_sst_writeblinkfile endif ; get and strip SST wavelengths of this file for labels (@ add Stokes) if (wavlabels ne 0) then begin restore,sstdir+'/'+sstspecsavs[ifilesst] nsstwavs=n_elements(spect_pos) sstlabels=strarr(nsstwavs) if (nsstwavs eq 0) then begin print,' ##### sdo_sst_writeblinkfile: no spect_pos in sstspecsav' print,' so no SST labels' endif if (sstlcs[0] eq 0) then sstlc=fix(nsstwavs/2) $ else sstlc=sstlcs[ifilesst] if (spect_pos[sstlc] ne fix(spect_pos[sstlc])) then $ print,' ##### sdo_sst_writeblinkfile: sstlc = '+trimd(sstlc)+$ ' wrong? Specify keyword sstlcs. No SST labels for now' sstwavs=fix((spect_pos-spect_pos[sstlc])*1000.)/1000. species='' if strmatch(sstfile,'*6563*') then species='Ha' if strmatch(sstfile,'*8542*') then species='Ca' if strmatch(sstfile,'*6302*') then species='Fe' for ilabel=0,nsstwavs-1 do begin if (sstwavs[ilabel] ge 0) then $ sstlabels[ilabel]=species+'+'+trim(sstwavs[ilabel]) else $ sstlabels[ilabel]=species+trim(sstwavs[ilabel]) endfor labels=[labels,sstlabels] endif endfor ; end loop over SST files for parameter definition ; set xrange, yrange, trange for desired subcube or the full monty xmin=xrange[0] xmax=xrange[1] ymin=yrange[0] ymax=yrange[1] tmin=trange[0] tmax=trange[1] if (xmax eq -1) then xmax=nx_sdo-1 if (ymax eq -1) then ymax=ny_sdo-1 if (tmax eq -1) then tmax=nt_sdo-1 nx_out=xmax-xmin+1 ny_out=ymax-ymin+1 nt_out=tmax-tmin+1 ; ---------- start big loop over all "wavelength" cubes to combine nallfiles=nsdofiles+nsstfiles nallwavcubes=nsdofiles+fix(total(nw_sst)) ; SST split per wavelength ifile=-1 for iwavcube=0,nallwavcubes-1 do begin ; first set up assoc for the SDO files if (iwavcube lt nsdofiles) then begin ifile=ifile+1 ; each SDO file is a separate wavcube iwavfilestart=iwavcube get_lun, unit_sdo ;RR check machine endian ;RR fits files are big endian alas, but my intel machine is little endian if (BYTE(1,0,1) eq 1) then my_endian='l' else my_endian='b' if (my_endian eq 'b') then $ openr,unit_sdo,sdofiles[iwavcube] else $ openr,unit_sdo,sdofiles[iwavcube],/swap_if_little_endian if (bitpix_sdo eq -32) then $ inassoc=assoc(unit_sdo,fltarr(nx_sdo,ny_sdo),header_sdosize) if (bitpix_sdo eq 16) then $ inassoc=assoc(unit_sdo,intarr(nx_sdo,ny_sdo),header_sdosize) if (bitpix_sdo eq 8) then $ inassoc=assoc(unit_sdo,bytarr(nx_sdo,ny_sdo),header_sdosize) if (verbose) then print,' ===== file '+ntostr(iwavcube+1)+$ ' = '+sdofiles[iwavcube] ;; ; check ;; sv,inassoc[0] ;; STOP endif ; ----- switch to SST with split of files into cubes per wav if (iwavcube ge nsdofiles) then begin ; find whether to open new SST file if (iwavcube eq nsdofiles) then begin free_lun,unit_sdo ifilesst=0 ifile=ifile+1 iwavfilestart=iwavcube newsstfile=1 iwavsst0=nsdofiles endif else newsstfile=0 if (iwavcube eq iwavsst0+nw_sst[ifilesst]) then begin free_lun,unit_sst iwavsst0=iwavsst0+nw_sst[ifilesst] ifilesst=ifilesst+1 ifile=ifile+1 iwavfilestart=iwavcube if (ifilesst eq nsstfiles) then goto, ALLDONE newsstfile=1 endif ; open SST file if (newsstfile eq 1) then begin if (verbose) then $ print,' ===== file '+ntostr(nsdofiles+ifilesst+1)+' = '+$ sstfiles[ifilesst] sstfile=sstdir+'/'+sstfiles[ifilesst] header_sstsize=512 get_lun, unit_sst if (endian_sst eq my_endian) then $ openr,unit_sst,sstfile else $ openr,unit_sst,sstfile,/swap_if_little_endian if (bitpix_sst[ifilesst] eq -32) then $ inassoc=assoc(unit_sst,fltarr(nx_sst,ny_sst),header_sstsize) if (bitpix_sst[ifilesst] eq 16) then $ inassoc=assoc(unit_sst,intarr(nx_sst,ny_sst),header_sstsize) if (bitpix_sst[ifilesst] eq 8) then $ inassoc=assoc(unit_sst,bytarr(nx_sst,ny_sst),header_sstsize) ;; ; check ;; sv,inassoc[0] ;; STOP ; check first image size of this new SST file against SDO im0=inassoc[0] sizeim0=size(im0) nx_im0=sizeim0[1] ny_im0=sizeim0[2] if (nx_im0 ne nx_sdo or ny_im0 ne ny_sdo) then begin print, ' ##### sdo_sst_writeblinkfile abort: '+$ ' first image cube '+trimd(iwavcube)+' [nx,ny] inequal to SDO' return endif endif ; end new SST file endif ; -------- end setup next SST wavcube ; ========== from here on for both SDO and SST wavcubes ; if bytscale > 1 set min and max (1st image or whole cube or subcube) if (bytscale gt 1) then begin cubemin=1.e10 cubemax=-1.e10 itstart=0 itend=nt_sdo-1 if (bytscale eq 2) then itend=itstart if (bytscale eq 4) then begin itstart=tmin itend=tmax endif for it=itstart,itend do begin image=inassoc[(it+itstart)*nwavsfile[ifile]+iwavcube-iwavfilestart] if (bytscale eq 4) then image=image[xmin:xmax,ymin:ymax] if (min(image) lt cubemin) then cubemin=min(image) if (max(image) gt cubemax) then cubemax=max(image) endfor ; maintain zero value for SDO magnetograms and Doppler maps if (iwavcube lt nsdofiles) then begin if (sdowavs[iwavcube] eq 'mag' or sdowavs[iwavcube] eq 'dop') then begin minzero=min([cubemin,-cubemax]) maxzero=max([cubemax,-cubemin]) limit=max([abs(minzero),maxzero]) cubemax=limit cubemin=-limit endif endif endif ; initial SDO file: set up and start image output file if (iwavcube eq 0) then begin ; set output parameters nz_out=long(nt_out)*nallwavcubes ; all 3rd-axis images ;RR ?? accommodate SST float fcubes too? Needed for Stokes if (bytscale eq 0) then bitpix_out=16 else bitpix_out=8 blank=bytarr(nx_out,ny_out) ; define outfile filename including specification of dimensions outfiletail='_xyt_'+ntostr(xmin)+'_'+ntostr(ymin)+'_'+ntostr(tmin)+$ '_nxyt_'+ntostr(nx_out)+'_'+ntostr(ny_out)+'_'+ntostr(nt_out)+$ '_wavs_'+ntostr(nsdofiles) for jsst=0,nsstfiles-1 do outfiletail=outfiletail+'_'+ntostr(nw_sst[jsst]) if (bitpix_out eq 8) then outfileext='.bcube' if (bitpix_out eq 16) then outfileext='.icube' outfile=outfilestart+outfiletail+outfileext ; define header_out (La Palma format crispex image file: [nx,ny,nt*nw*ns] ;RR ns=1, @@ still to treat ns=4 (full Stokes) if (bitpix_out eq 8) then im0=bytarr(nx_out,ny_out) if (bitpix_out eq 16) then im0=intarr(nx_out,ny_out) header_out=make_lp_header(im0,nt=nz_out) header_out=str_replace(header_out,' datatype',$ 'stokes=[I], ns=1, datatype') header_outsize=512 ; open output file for assoc (write little_endian La Palma file) get_lun, unit_out openw,unit_out,outfile,/swap_if_big_endian if (bitpix_out eq 8) then $ outassoc=assoc(unit_out,bytarr(nx_out,ny_out),header_outsize) if (bitpix_out eq 16) then $ outassoc=assoc(unit_out,intarr(nx_out,ny_out),header_outsize) ; write header_out rec=assoc(unit_out, bytarr(header_outsize)) rec[0]=byte(header_out) ; end startup output at first SDO file endif ; set wavlabels overlay image for this wavcube if (wavlabels ne 0) then begin while !d.window ne -1 do wdelete,!d.window window,xsize=nx_out,ysize=ny_out,/pixmap,retain=2 !p.font=-1 tv,blank xyouts,ny_out/30,ny_out/30,labels[iwavcube],color=255,$ /device,charsize=sqrt(ny_out)/10.,charthick=sqrt(ny_out)/10. imwavlabels=tvrd() endif ; get and rewrite all images of this iwavcube iw-distributed into outfile for it=0,nt_out-1 do begin image=inassoc[(it+tmin)*nwavsfile[ifile]+iwavcube-iwavfilestart] image=image[xmin:xmax,ymin:ymax] if (bytscale eq 1) then image=bytscl(image) if (bytscale gt 1) then image=bytscl(image,min=cubemin,max=cubemax) if (wavlabels ne 0) then image[where(imwavlabels gt 1)]=max(image) izim=it*nallwavcubes+iwavcube outassoc[izim]=image ;; ; check ;; sv,outassoc[0] ;; STOP endfor endfor ; end big loop over all iwav subcubes per wavelength ALLDONE: ; free the output file free_lun,unit_out if (verbose ne 0) then $ print,' ===== sdo_sst_writeblinkfile wrote image file: '+outfile ; write optional _sp spectrum file (may be very slow) if (spfile ne 0) then begin print,' ===== writing also spectral file: '+$ outfilestart+outfiletail+'_sp'+outfileext lp_im2spfile,outfile,nlp=nallwavcubes,ns=1,verbose=0 endif ; print elapsed time timedone=systime(1) timelaps=ntostr((timedone-timestart)/60.,format='(F11.1)') if (verbose ne 0) then print,$ ' ===== sdo_sst_writeblinkfile took '+timelaps+' minutes' end ; ==================== test per IDLWAVE super C ====================== ;; cd,'/media/rutten/SSTDATA/alldata/SST/2014-06-21-quiet' ;; sdofiles=['aia1600_algbo.fits','aia304_algbo.fits','aia171_algbo.fits'] ;; sdofiles='all' ;; sdofiles=['aia304_algbo.fits','aia171_algbo.fits'] ;; sstdir='crispex' ;; sstfiles=['crispex.6563.08:02:31.time_corrected.aligned.icube',$ ;; 'crispex.8542.08:02:31.time_corrected.icube'] ;; sstspecsavs=['spectfile.6563.idlsave',$ ;; 'spectfile.8542.idlsave'] ;; ; small subcube for fast tests ;; xrange=[400,599] ;; yrange=[400,599] ;; trange=[90,99] ; 0.009 min ;; ; subcube = contrail cutout ;; xrange=[300,700] ;; yrange=[400,700] ;; trange=[50,150] cd,'/media/rutten/SSTDATA/alldata/SST/2016-09-05-demo' sdofiles=['aia1700_alfor_px.fits','aia1600_alfor_px.fits',$ 'aia304_alfor_px.fits','aia171_alfor_px.fits'] sstfiles=['crispex.6563.09:48:31.time_corrected.aligned.icube',$ 'crispex.8542.09:48:31.stokesI.icube'] sstspecsavs=['spectfile.6563.idlsave',$ 'spectfile.8542.idlsave'] bytscale=4 wavlabels=1 spfile=0 verbose=1 ; delete previous ones except full ones if (n_elements(xrange) eq 0) then xmin=-1 else xmin=xrange[0] spawn,'rm -f sdo2sst/sdo_sst_xyt_'+ntostr(xmin)+'_*cube' ; do the work sdo_sst_writeblinkfile,sdofiles,sstfiles,$ outfilestart=outfilestart,spfile=spfile,$ sdodir=sdodir,sstdir=sstdir,$ sstspecsavs=sstspecsavs,sstlcs=sstlcs,bytscale=bytscale,$ xrange=xrange,yrange=yrange,trange=trange,$ wavlabels=wavlabels,verbose=verbose ; get result file(s) outfile=findfile('sdo2sst/sdo_sst_xyt_'+ntostr(xmin)+'*cube') if (n_elements(outfile) eq 2) then begin imfile=outfile[0] spfile=outfile[1] endif else imfile=outfile ; show with showex showex,imfile ; try crispex (only useful when bytscale=0) ;; crispex_2016,imfile,spfile end