; file: fitsimages2fitscube.pro = write cube.fits from loose fits images ; init: Feb 7 2013 Rob Rutten LA Deil ; last: Jan 9 2021 Rob Rutten Deil ; note: version gongimages2fitscube.pro in gong_sdo.pro; nowhere used ;+ pro fitsimages2fitscube,inpath,filestring,cubefile,$ timefile=timefile,swap=swap,$ xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,$ itstart=itstart,itend=itend ; convert a seqence of fits file images to a cube in a fitsfile ; ; INPUTS: ; inpath = string with path to input fits files (end with /) ; filestring = string with selective part of file names, eg _1700 ; cubefile = string with output path/filename ending with .fits ; ; OPTIONAL KEYWORD INPUTS: ; timefile: string filename to date_obs as TAI, include 'tai' ; swap=1/0: byteswap for wrong endian (default 0) ; xmin,xmax,ymin,ymax,itstart,itend: select partial cube ; ; OUTPUTS: ; fits file in cubefile, timefile ; ; METHOD: ; assoc so can treat large file lists ; ; HISTORY: ; Feb 7 2013 RR: start ; Nov 10 2013 RR: refound via Google, retrieved, corrected ; Nov 10 2013 RR: changed to assoc ; Nov 10 2013 RR: corrected to maintain data type ; Oct 11 2015 RR: swap ; Sep 16 2020 RR: uncompress .fz, timefile ;- ; answer a no-parameter query if (n_params() lt 3) then begin sp,fitsimages2fitscube return endif if (n_elements(swap) eq 0) then swap=0 if (n_elements(timefile) eq 0) then timefile='/tmp/fitsimagestimestai.dat' if (n_elements(xmin) eq 0) then xmin=0 if (n_elements(xmax) eq 0) then xmax=0 if (n_elements(ymin) eq 0) then ymin=0 if (n_elements(ymax) eq 0) then ymax=0 if (n_elements(itstart) eq 0) then itstart=0 if (n_elements(itend) eq 0) then itend=0 ; check timefile on including tai if (strmatch(timefile,'*tai*') ne 1) then begin print,' ##### timefile name should include "tai"' return endif ; remove earlier timefile spawn,'rm -f '+cubefile+' '+timefile ; uncompress files when rice-compressed (.fz) spawn,'funpack -D '+inpath+'/*'+filestring+'*.fits.fz' ; get filelist files=findfile(inpath+'/*'+filestring+'*') ntfull=n_elements(files) if (itend eq 0) then itend=ntfull-1 nt=itend-itstart+1 ; get file size inheader=headfits_rr(files[0]) nxfull=fxpar(inheader,'naxis1') nyfull=fxpar(inheader,'naxis2') if (xmax eq 0) then xmax=nxfull-1 if (ymax eq 0) then ymax=nyfull-1 nx=xmax-xmin+1 ny=ymax-ymin+1 ; set endian, get data type, define outheader bigendian=1 bitpix=fxpar(inheader,'bitpix') mkhdr,outheader,abs(bitpix)/8,[nx,ny,itend-itstart+1] sizeoutheader=size(outheader) ; fits header = Nx36 "card images" = Nx2880 bytes outheadersize=(1+fix(sizeoutheader[1]/36.))*2880 ; open output file for assoc, write header get_lun, unit_out if (bigendian) then openw,unit_out,cubefile,/swap_if_little_endian $ else openw,unit_out,cubefile if (bitpix eq -32) then outassoc=assoc(unit_out,fltarr(nx,ny),outheadersize) if (bitpix eq 16) then outassoc=assoc(unit_out,intarr(nx,ny),outheadersize) if (bitpix eq 8) then outassoc=assoc(unit_out,bytarr(nx,ny),outheadersize) if (outheadersize ne 0) then begin rec=assoc(unit_out, bytarr(outheadersize)) rec[0]=byte(outheader) endif ; set up timing array if (timefile ne '') then timarr=dblarr(nt) ; loop over timestep it for it=itstart,itend do begin ; read next fitsfile image=readfits(files[it],head) ; assoc image into output cubefile if (swap ne 0) then image=swap_endian(image) outassoc[it-itstart]=image[xmin:xmax,ymin:ymax] ;; print,'==== time step ',ntostr(it) ; add time into timing array if (timefile ne '') then begin datetime=fxpar(head,'date_obs') ; GONG uses date-obs if (datatype(datetime) ne 'STR') then datetime=fxpar(head,'date-obs') if (datatype(datetime) ne 'STR') then begin print,' ##### no DATE_OBS or DATE-OBS keyword in header' return endif timarr[it-itstart]=anytim2tai(datetime) endif endfor ; end loop over files ; write timefile writecol,timefile,timarr,fmt='(E21.12)' ; free the input and output files free_lun,unit_out print,' ===== fits2cubefile wrote ',cubefile end ; =============== test hyper-C =============== cd,'/home/rutten/data/ALMA/2018-04-12-ca/gong/' indir='.' selstring='20180412' outfile='gongall.fits' timefile='gongalltimestai.dat' fitsimages2fitscube,indir,selstring,outfile,timefile=timefile ; inspect showex,outfile readcol,timefile,times,format='D',/silent print,anytim2utc(times[0:10],/ccsds) ; 20s GONG cadence, each station 1 min end