; raw2sgy.pro ; IDL procedure to convert raw radargrams produced by COSHARPS (JPL and ; Smithsonian processors) or from the PDS (SHOC and Smithsonian processors) ; to SEG-Y format for use in seismic data processing & interpretation software. ; SEG-Y specifications call for a 3600-byte header at the beginning of the file ; and 240-byte headers at the beginning of every trace (i.e., frame). ; Smithsonian (UPB & FPB) and SHOC output is stored as an [nframes,nsamples] ; (i.e., all t=time[0] samples, followed by all t=time[1] samples, etc.) ; whereas JPL output is stored as an array [nsamples,nframes] (i.e., all ; samples for frame f=1, followed by all samples for frame f=2, etc.). ; Smithsonian "raw" decoded data (RPB) is stored as [nsamples,nframes]. ; In contrast, the SEGY structure is [header,[header+nsamples, nframes]] ; (i.e., file header followed by a trace header, the first trace samples, ; the second trace header, the second trace samples, etc.) ; 2007 Jan 25 created by Than Putzig ; 2007 Jan 31 Than Putzig, implemented multichannel (RGB) stuff ; 2007 Feb 08 Than Putzig, added call to orb2stn.pro ; 2007 Apr 17 v01 Than Putzig, added prodir, outdir ; 2007 May 22 v02 Than Putzig, updated for new directory and file naming; ; automated dataset issue, process version ; 2007 Oct 08 v03 Than Putzig, for UPA, FPA, and FPS (SHOC/PDS) formats ; Optional override for LINEID ; 2007 Oct 31 Than Putzig, Allow OBS or UPA prefix for UPA products ; 2007 Nov 19 v04 Than Putzig, Allow 'product' to be file name ; 2008 Jan 16 Than Putzig, fixed problem with OBS allowance ; 2008 Jan 23 v05 Than Putzig, Use multilook FPAs/QDAs (.mlk in place of .foc); ; added wti to orb2stn call ; 2008 Apr 17 v06 Than Putzig, fixed problem with NTR, added TRO, ODAT ; 2008 May 21 v07 Than Putzig, Added FOC switch, modified mlk handling since mlk ; files are not complex after all. ; 2008 Jun 16 v08 Than Putzig, Use assoc files for temporary storage and ; handle frames in FPI=10000 chunks for speed; ; added RDR keyword for FPS preprocessing ; 2008 Jun 27 v09 Than Putzig, Modified calls to RDR2SAN to obtain RANGE_SHIFTs ; for multi-mode observations ; 2009 Jun 27 Than Putzig, Fixed bug with filename inputs ; 2009 Jul 30 v10 Than Putzig, Added /xy keyword for orb2stn ; 2009 Nov 17 v11 Than Putzig, Added FPB functions ; 2010 Apr 08 Than Putzig, Fixed handling when product includes e.g. _1_01 ; 2010 May 10 v12 Than Putzig, Allow lower-case name for raw FPS data; ; pass lineid as pid to rdr2san ; 2010 Jun 02 v13 Than Putzig, Applied NTR, TRO to temporary files ; 2011 Jul 19 v14 Than Putzig, Correct obs>2147401000 problem (2^31=2147483648): ; LINEID now written to 3200-byte text header and ; 32-bit binary header word written as lineid/1000 ; 2011 Nov 04 v15 Than Putzig, Changed default to FPB with POW=1 ; 2011 Nov 09 v15 Than Putzig, Added SUPDIR ; 2012 May 24 v15 Than Putzig, fixed problem with PRODUCT variable mangling ; 2012 May 24 v16 Than Putzig, added prompting for missing files ; 2012 May 30 v17 Than Putzig, allow proc='QDA' ; 2013 Feb 06 v18 Than Putzig, change default PROC to FPS for external use ; 2013 Feb 06 v19 Than Putzig, set default PROC back to FPB for internal use ; 2014 Jan 06 v20 Than Putzig, to read PDS variant of Smithsonian FPB products, ; fix partial-file read/write and trace seq no. ; 2014 Jan 26 v21 Than Putzig, use pixlatlon rather than Orbit files (UPB,FPB) ; 2014 Apr 22 v22 Than Putzig, Set default PROC by test of FPB archive existence ; 2014 May 29 v23 Than Putzig, Modify default PROC if PRODUCT is a file (PIF); ; fix basename of pixlatlon for PIF cases. ; 2014 Nov 13 v23 Than Putzig, Fix default PDS PROC if PRODUCT is a file (PIF); ; added /radii keyword for orb2stn. ; 2015 Feb 27 v24 Than Putzig, Fixed data sample format header (IEEE, not IBM). ; 2015 Jun 17 v25 Than Putzig, Updates to handle COSHARPS FPB v9.1 (PDS-style). ; Added indir parameter. ; 2015 Jul 03 v26 Than Putzig, Updates to handle UTS (simulations). Will call ; sim2raw.pro if needed. ; 2015 Jul 31 v27 Than Putzig, Added proc=RPB handling (Grima InSight support). ; 2016 Nov 29 v28 Than Putzig, to pass clat and clon to orb2stn ; 2016 Dec 01 v29 Than Putzig, workaround for strsplit count keyword in GDL pro raw2sgy,product,odat,prodir=prodir,supdir=supdir,indir=indir,outdir=outdir,$ type=type,modes=modes,ns=ns,proc=proc,ntr=ntr,tro=tro,ibs=ibs,obs=obs,$ fpi=fpi,dns=dns,dso=dso,dsf=dsf,dst=dst,ct=ct,r=r,g=g,b=b,lineid=lineid,$ rdr=rdr,nostn=nostn,display=display,nowrite=nowrite,png=png,ali=ali,$ pow=pow,foc=foc,pad=pad,wti=wti,null=null,diag=diag,odiag=odiag,xy=xy,$ proj=proj,clat=clat,clon=clon,pds=pds,radii=radii,help=help pron='RAW2SGY' prom='['+pron+']: ' version='29' if n_elements(product) ne 0 then begin ; process PRODUCT parameter and use to set default PROC sproduct=strcompress(product,/rem) ; Test whether PRODUCT is a file name. PIF=file_test(sproduct,/regular) endif else PIF=0 if PIF then begin aproc=['FPB','UPB','RPB','UPA','FPA','QDA','FPS','UTS'] sproc=strmid(file_basename(sproduct),0,3) ; get first three chars tproc=strlowcase(strmid(file_basename(sproduct),0,2)); get first two chars wproc=where(sproc eq aproc,cproc) dproc=(cproc eq 1)?aproc[wproc]:((tproc eq 'r_')?'FPS':'FPB') ;if n_elements(PDS) eq 0 and cproc gt 2 then PDS=1 if n_elements(PDS) eq 0 and (cproc eq 0 or wproc eq 6) then PDS=1 endif else begin ;dproc='FPS' ; SHOC Focused Processor is default v18 ;dproc='FPB' ; Focused Processor "B" is default v19 fpbgpd='/data/g/WUSHARPS/Archive/FPB_PROD/' dproc=(file_test(fpbgpd,/directory))?'FPB':'FPS' endelse if n_elements(proc) eq 0 then proc=dproc debug=0 arcv=(proc eq 'UPB')?'d':((proc eq 'FPB')?'g':$ ((proc eq 'QDA' or proc eq 'UTS')?'j1':'e')) ppfx=(proc eq 'UTS')?'QDA':proc dsupdir='/data/'+arcv+'/WUSHARPS/Archive/'+ppfx+'_PROD/' if n_elements(outdir) eq 0 then outdir='.' IF N_PARAMS() EQ 0 OR KEYWORD_SET(HELP) THEN BEGIN PRINT,"" PRINT,"Version ",version PRINT,pron +",PRODUCT{,ODAT,PRODIR=,SUPDIR=,INDIR=,OUTDIR=,TYPE=,PROC=," PRINT," LINEID=,MODES=,NS=,IBS=,OBS=,NTR=,TRO=,FPI=,DNS=,DSO=,DSF=," PRINT," DST=,CT=,/PDS,/XY,PROJ=,CLAT=,CLON=,/R,/G,/B,/RDR,/NOSTN," PRINT," /NOWRITE,/DISPLAY,/FOC,/PAD,/POW,/ALI,/PNG,/DIAG,/ODIAG,/HELP}" IF KEYWORD_SET(HELP) THEN BEGIN PRINT," Parameter Description {default value}" PRINT," --------------------------------------------------------" PRINT," PRODUCT raw-file name or CO-SHARPS product ID (REQ.) {}" PRINT," If CO-SHARPS product, e.g., '219001000_1_01'," PRINT," then if '_1_01' not included, will use highest" PRINT," available dataset issue and/or process version" PRINT," ODAT returned array of data {}" PRINT," SUPDIR CO-SHARPS processed data super-directory" PRINT," {'"+dsupdir+"'}" PRINT," PRODIR CO-SHARPS processed data parent directory" PRINT," "+$ "{'SUPDIR/"+ppfx+"_PROD_???00'}" PRINT," or FPS product directory {.}" PRINT," INDIR input directory {}" PRINT," OUTDIR output directory {"+$ outdir+"}" PRINT," PROC processor (FPB,UPB,RPB,UPA,FPA,QDA,FPS,UTS) {'"+$ proc+"'}" PRINT," TYPE UPB/FPB product type {'unfocpow/focpow'}" PRINT," default is 'unfocgram' if 'unfocpow' not found " PRINT," and /R, /G, /B assumed if TYPE is 'unfocgram' " PRINT," LINEID Override output LINEID {extracted from PRODUCT}" PRINT," MODES e.g., [5,6,7] {ALL}" PRINT," NS is number of samples per frame" PRINT," {U/F/RPB: 3600; U/FPA,QDA: 4096; FPS:1000}" PRINT," For UTS, default is taken from file name" PRINT," IBS is input bytes per sample {UTS:1;FOC:8;other:4}" PRINT," OBS is output bytes per sample {4}" PRINT," NTR is total number of traces to output {ALL}" PRINT," TRO output trace offset {0}" PRINT," PID is override for product ID {from PRODUCT}" PRINT," FPI frame processing interval (for memory) {10000}" PRINT," DNS is number of samples to display {ALL}" PRINT," DSO is sample offset to begin display {0}" PRINT," DSF is frame scale factor for display {1}" PRINT," DST is time scale factor for display {1}" PRINT," NULL null value, also passed to RDR2SAN {-99.}" PRINT," CT color table to use for display {3, red temp}" PRINT," /PDS PRODUCT is from PDS {FPB: 0; FPS: 1}" PRINT," /XY produce X and Y projected values from orb2stn " PRINT," PROJ IDL projection number (106=polar stereo) {106}" PRINT," CLAT Center latitude for projection {90}" PRINT," CLON Center longitude for projection {0}" PRINT," /R extract red (forward) channel (UPB only) {0}" PRINT," /G extract green (nadir) channel (UPB only) {1}" PRINT," /B extract blue (aft) channel (UPB only) {0}" PRINT," /RDR process FPS input through RDR2SAN first {0}" PRINT," /NOSTN no seq trace # file output {PDS: 1; others: 0}" PRINT," /NOWRITE to disable writing of output file {0}" PRINT," /DISPLAY to display input radargram {0}" PRINT," /PNG to save a PNG file of displayed radargram {0}" PRINT," /FOC use complex .foc files instead of real .mlk {0}" PRINT," /PAD retain padding in .foc files (diagnostic) {0}" PRINT," /POW to use power scaling {R/FPB: 1; others: 0}" PRINT," /ALI to use Ali-style scaling {0}" PRINT," /DIAG to print diagnostic messages {0}" PRINT," /ODIAG to print orb2stn diagnostic messages {0}" ENDIF PRINT, " /HELP prints parameter details" PRINT,"" RETURN STOP ENDIF ; set processor flags FPS=(proc eq 'FPS') UPB=(proc eq 'UPB') FPB=(proc eq 'FPB' or proc eq 'RPB') RPB=(proc eq 'RPB') UTS=(proc eq 'UTS') PDS=(n_elements(PDS) ne 0) ? PDS:((FPS) ? 1:0) xy=keyword_set(xy) ; orb2stn xy flag radii=keyword_set(radii) ; orb2stn radii flag ; set the no-sequential-trace-number-file flag nostn=(n_elements(nostn) ne 0) ? nostn:((FPS and not xy and not radii) ? 1:0) rdr=keyword_set(rdr) ; flag to run rdr2san display=keyword_set(display) or keyword_set(png) ; radargram display flag nowrite=keyword_set(nowrite) ; radargram write flag png=keyword_set(png) ; png write flag ali=keyword_set(ali) ; use Ali-style scaling foc=keyword_set(foc) ; use .foc, not .mlk pad=keyword_set(pad) ; include padded frames ;pow=keyword_set(pow) ; use power scaling if n_elements(pow) eq 0 then pow = (FPB) ? 1 : 0 ; default power scaling diag=keyword_set(diag) ; diagnostics flag odiag=keyword_set(odiag) ; orb diagnostics flag if n_elements(null) eq 0 then null=-99. ; set default frame-handling interval if n_elements(fpi) eq 0 then fpi=10000 ; Flag JPL and Smithsonian (SI) processors JPL=(proc eq 'UPA' or proc eq 'FPA' or proc eq 'QDA') SI=(proc eq 'UPB' or proc eq 'FPB' or proc eq 'RPB') ; set write trace interval for orb2stn routine ; - use 2 for FPA because .etm files are out of sync with .mlk files ;; Update 2008 May 20: ;; Not actually true -- .mlk files are NOT complex, they're all real. ;;wti=(proc eq 'FPA') ? 2 : 1 ;;wti=1 if n_elements(wti) eq 0 then wti=(UTS)?8:1 mode=(JPL) ? 'Mode_' : '' ; location file suffix lsfx=(JPL or UTS) ? '.dat.etm' : ((FPS) ? '.latlon' : ((PDS) ? '.tab' : '.txt')) ;otxt=(SI) ? ((PDS) ? 'GEOM' : 'Orbit_' ) : '' ;otxt=(SI) ? ((PDS) ? '_GEOM' : '_pixlatlon_' ) : '' otxt=(SI) ? ((PDS) ? '_geom' : '_pixlatlon_' ) : '' nm=n_elements(modes) ; count # of modes ; default input product type if n_elements(type) eq 0 then $ type = (UPB) ? 'unfocpow' : ((FPB) ? ((PDS) ? 'rgram' : 'focpow') : '') ; set input raw data suffix case proc of 'UPB': sfx='.raw' 'UTS': sfx='.raw' 'FPB': sfx=(PDS)?'.IMG':'.raw' 'RPB': sfx=(PDS)?'.IMG':'.raw' 'UPA': sfx='.dat.cmp' 'FPA': sfx='.dat'+((foc)?'.foc':'.mlk') 'QDA': sfx='.dat'+((foc)?'.foc':'.mlk') 'FPS': sfx=(rdr)?'_A.DAT':'.dat' endcase simpng=0 ; flag for PNG simulation input ; ; process PRODUCTs given at command line if PIF then begin rfn=sproduct ; set raw file search string sfxpos=strpos(sproduct,sfx) if sfxpos eq -1 and PDS then begin ; try lowercase suffix for PDS products sfx=strlowcase(sfx) sfxpos=strpos(sproduct,sfx) endif if sfxpos eq -1 and UTS then begin ; look for PNG as simulation input sfx='.png' sfxpos=strpos(sproduct,sfx) if sfxpos ne -1 then simpng=1 endif if sfxpos eq -1 then begin print,'Expected suffix '+sfx+' not found in '+sproduct+'.' title='Please select appropriate raw data file(s):' sproduct=dialog_pickfile(title=title,filter='*'+sfx,/fix_filter) if sproduct eq '' then begin print,'No raw data file selected; exiting.' return stop endif rfn=sproduct sfxpos=strpos(sproduct,sfx) endif ; strip directory and suffix from product name sproduct=file_basename(strmid(sproduct,0,sfxpos)) ; run RDR2SAN on FPS product file if rdr keyword given if rdr then $ rdr2san,rfn,/binary,/location,null=null,diag=diag,dir=outdir,$ pid=lineid ; run sim2raw if PNG simulation input if simpng eq 1 then sim2raw,rfn,diag=diag pprod=sproduct ; strip to first- or second-to-last underscore for FPS and SI products if SI or FPS then begin pprod=strmid(sproduct,0,strpos(sproduct,'_',/reverse_search)) ;if not FPS then pprod=strmid(pprod,0,strpos(pprod,'_',/reverse_search)) if not PDS then pprod=strmid(pprod,0,strpos(pprod,'_',/reverse_search)) endif ;obn=pprod+otxt+((SI and not PDS)?'*':'')+lsfx ; set loc file search string obn=((UTS)?'QDA_':pprod+otxt)+'*'+lsfx ; set loc file search string if rdr then obn=outdir+'/'+obn ; get prefix of input product pfx=strmid(sproduct,0,strpos(sproduct,'_')) ; had problems parsing a non-standard prefix (FFF); ; check string length instead of specific strings ;;if not FPS and (pfx eq 'FPA' or pfx eq 'UPA' or $ ;; pfx eq 'UPB' or pfx eq 'FPB' or pfx eq 'OBS') then $ ;;if not FPS and strlen(pfx) eq 3 then begin ; Assuming here that <=3 chars is a bona fide prefix, including PDS data if strlen(pfx) le 3 or RPB then begin opfx=pfx sproduct=strmid(sproduct,strpos(sproduct,'_')+1) ; strip prefix endif else opfx=proc endif else begin if n_elements(indir) eq 0 then begin ; process PRODUCTs stored in CO-SHARPS database ; get dataset issue and process version if not given in PRODUCT ;spp=strsplit(strcompress(sproduct,/rem),'_',count=count,/extract) spp=strsplit(strcompress(sproduct,/rem),'_',/extract) count=n_elements(spp) ddir=ppfx+'_PROD_'+string(long64(spp[0])/10000000*100,format='(I05)') if n_elements(prodir) eq 0 then begin if n_elements(supdir) eq 0 then supdir=dsupdir if strmid(supdir,0,1,/rev) ne '/' then supdir=supdir+'/' ; Use current directory for PDS data ;prodir = (PDS) ? '.' : supdir+ddir prodir = (FPS) ? '.' : supdir+ddir endif if diag then print, "product directory is ", prodir ;if count lt 3 and PDS eq 0 then begin if count lt 3 and FPS eq 0 then begin if diag then print,"searching for ", prodir+'/'+ppfx+'_'+sproduct+'*' edir=file_basename(file_search(prodir+'/'+ppfx+'_'+sproduct+'*',$ /test_dir,count=ndir)) if edir[0] eq '' then begin print, prodir+'/'+ppfx+'_'+sproduct+'* not found. Exiting.' return stop endif ;spd=strsplit(edir[ndir-1],'_',count=nspd,/extract) spd=strsplit(edir[ndir-1],'_',/extract) nspd=n_elements(spd) if count eq 1 then sproduct=sproduct+'_'+spd[2] pdir=ppfx+'_'+sproduct+((nspd eq 4)?('_'+spd[3]):'') if SI then begin sproduct=sproduct+'_'+spd[3] if diag then print, "sproduct set to ", sproduct endif ; Use current directory for PDS data ;endif else pdir=(PDS eq 0)?ppfx+'_'+sproduct:'.' endif else pdir=(FPS eq 0)?ppfx+'_'+sproduct:'.' indir=prodir+'/'+pdir endif ; set input file prefix ;;prefix=indir+'/'+proc+'_'+sproduct+'_' ; use less specific prefix (first four digits of orbit) to accommodate ; Ali-processed FPAs (e.g. 6660_01_0_Mode_001.dat.foc) orbit=strmid(sproduct,0,4) rpref=indir+'/*'+orbit+'*' opref=((rdr)? outdir:(indir))+'/*'+orbit+'*' ;;ttype=type+((type ne '')?'_':'') ; set raw file search string ;;rfn=rpref+ttype+((nm ne 0) ? $ sim=(UTS)?'*-1024-Sim':'' rfn=rpref+type+((nm ne 0) ? $ (mode+string(modes,format='(I03)')):((UPB or FPB)?'*':''))+sim+sfx ; set location file search string obn=opref+otxt+'*'+((nm ne 0) ? $ (mode+string(modes,format='(I03)')):((UPB or FPB)?'???':''))+lsfx opfx=proc endelse ; Find raw file(s) rf=file_search(rfn,count=ni) if ni eq 0 then begin ; obsolete with looser prefix used above ; if JPL then strput, rfn, 'OBS', strpos(rfn[0],proc,/reverse_search) else $ if sfx eq '.dat.mlk' then begin sfx='.dat.foc' foc=1 strput, rfn, 'foc', strpos(rfn[0],'mlk') endif if type eq 'unfocpow' then begin strput, rfn, 'gr*', strpos(rfn[0],'pow') type='unfocgram' endif if PDS then begin rfn=strlowcase(rfn) sfx=strlowcase(sfx) endif if UTS then begin ; look for PNG as simulation input sfxpos=strpos(rfn,sfx) sfx='.png' rfn=strmid(rfn,0,sfxpos)+sfx simpng=1 endif rf=file_search(rfn,count=ni) if ni eq 0 then begin print,'Product '+rfn+' not found.' title='Please select appropriate raw data file(s):' rfn=dialog_pickfile(title=title,filter='*'+sfx,/fix_filter) rf=file_search(rfn,count=ni) if ni eq 0 then begin print, 'No valid product '+rfn+' found; exiting.' return stop endif endif endif ; determine maximum range shift in FPS input files if rdr and PIF eq 0 then begin mxsft=intarr(ni) for n=0,ni-1 do begin rdr2san,rf[n],range_shift mxsft[n]=max(range_shift) endfor shift=-1*max(mxsft) endif for n=0,ni-1 do begin if rf[n] eq '' then begin PRINT,"Input product "+rfn[n]+" not found." RETURN STOP endif ; run RDR2SAN on FPS product file if rdr keyword given if rdr and PIF eq 0 then begin ; set the trace number offset for latlon file if n gt 0 then begin sfxpos=strpos(rf[n-1],sfx) prevloc=outdir+'/'+file_basename(strmid(rf[n-1],0,sfxpos)+'.latlon') spawn, "tail -1 "+prevloc+" | awk '{ print $2 }'",fno fno=long(fno) endif else fno=0 rdr2san,rf[n],/binary,/location,null=null,dir=outdir,diag=diag,fno=fno,$ pid=lineid,shift=shift endif ; run sim2raw if PNG simulation input if simpng eq 1 and PIF eq 0 then sim2raw,rf[n],diag=diag endfor ; modify FPS *_A.DAT names with RDR output *.dat names if rdr then begin strput, rf, '.dat ', strpos(rf[0],sfx) rf=file_basename(strtrim(rf)) for n=0,ni-1 do begin RIF=file_test(outdir+'/'+rf[n],/regular) if RIF eq 0 then begin print,'RDR2SAN failed to produce '+rf[n]+'. Exiting.' return stop endif endfor endif ; modify UTS *.png names with sim2raw output *.raw names if simpng eq 1 then begin strput, rf, '.raw', strpos(rf[0],sfx) rf=file_basename(strtrim(rf)) for n=0,ni-1 do begin SIF=file_test(rf[n],/regular) if SIF eq 0 then begin print,'SIM2RAW failed to produce '+rf[n]+'. Exiting.' return stop endif endfor endif ; Find location file(s) if nostn eq 0 then begin if diag then print, "searching for ", obn ob=file_search(obn,count=ni,/fold_case) if ni eq 0 and (JPL or UTS) then begin strput,obn,'OBS',strpos(obn[0],proc,/reverse_search) ob=file_search(obn,count=ni) endif if ni eq 0 then begin print,'No location file '+obn+' found.' title='Please select appropriate location file(s):' obn=dialog_pickfile(title=title,filter='*'+lsfx,/fix_filter) ob=file_search(obn,count=ni) if ni eq 0 then begin print, 'No valid location file '+obn+' found; exiting.' return stop endif endif for n=0,n_elements(ob)-1 do begin if ob[n] eq '' then begin PRINT,"WARNING: Input location file "+obn[n]+" not found." endif endfor endif nc=(type eq 'unfocgram') ? 3 : 1 ; set # of channels omode=(nm eq 0) ? '' : strjoin('_'+string(modes,format='(I03)')) isiz=lon64arr(ni) for n=0,ni-1 do begin ; get file size in bytes rfi=file_info(rf[n]) isiz[n]=rfi.size endfor ; set number of samples per frame if n_elements(ns) eq 0 then begin if UTS then spin=strsplit(file_basename(rf[0]),'-',/extract) case proc of 'UPB': ns=3600 'FPB': ns=3600 'RPB': ns=3600 'UPA': ns=4096 'FPA': ns=4096 'QDA': ns=4096 'FPS': ns=1000 'UTS': ns=fix(spin[2]) endcase endif if n_elements(ibs) eq 0 then ibs=(foc)?8:((UTS)?1:4) ; input bytes/sample if n_elements(obs) eq 0 then obs =4 ; output bytes/sample if n_elements(r) eq 0 then r=0 ; red channel flag if n_elements(b) eq 0 then b=0 ; blue channel flag if n_elements(g) eq 0 then g=(r eq 0 and b eq 0) ; green channel flag if n_elements(ct) eq 0 then ct=3 ; color table if n_elements(dns) eq 0 then dns=ns else display=1 ; display # of samples if n_elements(dsf) eq 0 then dsf=1 else display=1 ; display scale factor if n_elements(dst) eq 0 then dst=1 else display=1 ; display scale factor if n_elements(dso) eq 0 then dso=0 else display=1 ; display sample offset ; SEGY filename otype=((type ne '')?'_':'')+type ;oname=outdir+"/"+((FPS and PIF) ? "":opfx+"_")+sproduct+((PIF)?"":otype+omode) oname=outdir+"/"+opfx+"_"+sproduct+((PIF) ? "":otype+omode) ofile=oname+".sgy" tfile=oname+".temp" ; SEQ TR NO filename ;ostfn=outdir+"/"+((FPS and PIF) ? "":opfx+"_")+sproduct+'_latlon'+omode+".txt" ostfn=outdir+"/"+((PIF)?pprod:(opfx+"_"+sproduct))+$ ((xy or radii) ? '_location':'_latlon')+omode+".txt" if rdr then begin spawn,"cat "+obn+" > "+ostfn,dummy if file_test(ostfn) then file_delete, obn endif if nostn eq 0 and file_test(ostfn,/write) then file_delete, ostfn ; ; line id e.g., 176902000 or 176902004 (mode 4) ; Jul 2011: Product id exceeded 32-bit at obs > 2147401000 ; lineid is now a 64-bit integer ; superceded code: ;if n_elements(lineid) eq 0 then lineid=((FPS) ? $ ; long(strmid(file_basename(rf[0]),2,$ ; strpos(file_basename(rf[0]),'_',3)-2))*1000 : $ ; long(strmid(sproduct,0,strpos(sproduct,'_'))))+((nm eq 1) ? modes[0] : 0) ; new code: if n_elements(lineid) eq 0 then begin if FPS then begin ;if PDS then begin rfb=file_basename(rf[0]) linestr=strmid(rfb,2,strpos(rfb,'_',3)-2) endif else begin ;if PDS and FPB then $ ; linestr=strmid(sproduct,0,strpos(sproduct,'_')) +$ ; strmid(sproduct,7,strpos(sproduct,'_',/reverse_search))$ ;else $ linestr=strmid(sproduct,0,strpos(sproduct,'_')) endelse lineid=long64(linestr)*((FPS)?1000:1) ; if a particular mode given, then add it to the lineid if nm eq 1 then lineid=lineid+modes[0] endif ; convert lineid to string for use in SEGY text header slineid=strcompress(lineid) if diag then print, "LineID set to ", slineid ; get series number from lineid series=lineid/100000000 tnt=total(isiz,/int)/ns/ibs/nc ; total # of traces ; For foc files (which are padded to integer multiple of 1024 frames), set ; ntr to number of lines in etm files, decimated by 4 for 5000 series obs, ; provided this value is smaller than that obtained from file sizes. if foc and pad eq 0 then tnt=total(file_lines(ob),/int)/((series eq 5)? 4:1) ; set number of traces to read if n_elements(ntr) eq 0 then ntr=tnt else ntr=long(ntr) ; # traces to read if n_elements(tro) eq 0 then tro=0l else tro=long(tro) ; trace offset ch=(nc eq 1 or r eq 1) ? 0 : ( (g eq 1) ? 1 : 2 ) ; channel to extract if diag then begin print, '' print, 'Number of input channels is '+strcompress(nc,/rem) if nc eq 3 then print, 'Extract from channel '+strcompress(ch+1,/rem)+$ ((ch eq 0) ? ' (RED)' : ((ch eq 1) ? ' (GREEN)' : ' (BLUE)')) print,"Will read "+strcompress(ntr,/rem)+" traces (frames) from "+$ strcompress(ni,/rem)+" file(s)" endif ; SEGY format calls for an integer sample interval in microseconds -- but radar ; comes in sub-microsecond sampling. So we give 10 X nanoseconds instead si=(FPS) ? 0.0750 : 0.0375 ; sample interval in us (37.5 ns) sih=fix(10000*si) ; sample interval for header (nsX10) lfh=3600l ; length of SEGY record header (bytes) lth=240 ; length of SEGY trace header (bytes) if nowrite eq 0 then begin openw,olun,ofile,/get_lun ; Open SEGY output file osiz=ntr*(lth+ns*obs)+lfh ; output file size oss=" of "+string(osiz,format='(I12)')+" bytes " ; +systime() ; CREATE SEGY file header fhead=bytarr(lfh) a=assoc(olun,fhead) a[0]=fhead ; write blank file header ; write lineid as string to textual header (first 3200 bytes) a=assoc(olun,byte(slineid)) ; byte offset description a[0]=byte(slineid) ; 0-12 line number (SHARAD product ID) ; write long-word file headers a=assoc(olun,[0l],3200) ; byte offset description a[1]=[lineid/1000] ; 3205-3208 "line number" (observation number) a[2]=[lineid/1000] ; 3209-3212 reel number (observation number) ; write short-word file headers a=assoc(olun,[0],3200) a[6]=[1] ; 3213-3214 data traces per record a[7]=[0] ; 3215-3216 auxillary traces per record a[8]=[sih] ; 3217-3218 sample interval (bogus) a[10]=[ns] ; 3221-3222 samples per trace ;a[12]=[1] ; 3225-3226 data sample format (4-byte IBM float) a[12]=[5] ; 3225-3226 data sample format (4-byte IEEE float) a[13]=[4] ; 3227-3228 Fold (hardwired to nominal pre/postsum) a[27]=[1] ; 3255-3256 Units (1=meters, 2=feet) if diag then begin print,"Wrote header to "+ofile fio=file_info(ofile) ;;csiz=strcompress(fio.size,/rem) csiz=string(fio.size,format='(I12)') print," "+csiz+oss endif ; CREATE SEGY trace header ltr=ns*obs+lth ; trace length with header thead=intarr(lth/2) ; initialize trace header ; write short-word trace headers ; byte offset description thead[14]=1 ; 29-30 trace id code (1=seismic) thead[57]=ns ; 115-116 number of samples per trace thead[58]=sih ; 117-118 sample interval (bogus) c=assoc(olun,intarr(ltr/2),lfh) ; associate output file as integer d=assoc(olun,lonarr(ltr/4),lfh) ; associate output file as long e=assoc(olun,fltarr(ltr/4),lfh) ; associate output file as float endif t=0l ; initialize sequential trace number openw,tlun,tfile,/get_lun ; open raw data file tdat=assoc(tlun,fltarr(tnt,ns)) ; associate raw data file tdat[0]=fltarr(tnt,ns) ; initialize raw data array for n=0,ni-1 do begin ; loop over input files nf=isiz[n]/ns/ibs/nc ; # of frames in current input file inty=(ibs eq 1)?'bytarr':'fltarr' ; input file array type if foc and pad eq 0 then nf=file_lines(ob[n])/((series eq 5)? 4:1) mfi=(nf-1)/fpi ; maximum frame interval index if diag then print,"Reading frames from "+file_basename(rf[n]) openr,ilun,rf[n],/get_lun ; open input file if n eq 0 and ali then begin if diag then print, "Determining Ali-style scaling factor" ; determining scale factor has to be done over intervals of FPI frames ; to avoid failures due to insufficient memory for large observations rds=fltarr(mfi+1) for fin=0l,mfi do begin ; loop over frame intervals fbi=fin*fpi ; begin interval frame fei=(fin lt mfi) ? fpi-1 : nf-fin*fpi-1 ; end interval frame bof=long64(fin)*fpi*ns*ibs*nc ; byte offset to data block ; Associate data block within JPL input files and entire file for ; other processors. Note that JPL processor output is complex ; (but not strictly IDL complexarr) and the position of samples ; and frames are switched relative to other processors. ex=execute("rdi=(JPL)?((foc)?assoc(ilun,"+inty+"(nc,ns,2,fei+1),bof):"+$ "assoc(ilun,"+inty+"(nc,ns,fei+1),bof)):" +$ "((RPB)?assoc(ilun,"+inty+"(nc,ns,nf)):"+$ "assoc(ilun,"+inty+"(nc,nf,ns)))") rds[fin]=max(abs( ( (JPL) ? ( (foc) ? $ complex(rdi[ch,*,0,0:fei,0],rdi[ch,*,1,0:fei,0]):$ rdi[ch,*,0:fei,0]): $ ((RPB)?rdi[ch,*,fbi:fbi+fei,0]:$ rdi[ch,fbi:fbi+fei,*,0]))))*1.e-2 if diag then $ print, " read ", strcompress(fin*fpi+fei+1,/rem), " frames" endfor rds=max(rds) if diag then print,"Ali-style scaling factor is ",rds endif if diag then print, "Transferring data to temporary file" for fin=0l,mfi do begin ; loop over frame intervals fbi=fin*fpi ; begin interval frame fei=(fin lt mfi) ? fpi-1 : nf-fin*fpi-1 ; end interval frame bof=long64(fin)*fpi*ns*ibs*nc ; byte offset to data block ; Associate data block within JPL input files and entire file for ; other processors. Note that JPL processor output is complex ; (but not strictly IDL complexarr) and the positions of samples ; and frames are switched relative to other processors. ex=execute("rdi=(JPL)?((foc) ? assoc(ilun,"+inty+"(nc,ns,2,fei+1),bof):"+$ "assoc(ilun,"+inty+"(nc,ns,fei+1),bof)):" +$ "((RPB)?assoc(ilun,"+inty+"(nc,ns,nf)):"+$ "assoc(ilun,"+inty+"(nc,nf,ns)))") ; only extract from intervals with frames of interest if t+fbi le tro+ntr and t+fbi+fei gt tro then begin if diag then $ print, " writing frames ", strcompress(t+fbi+1,/rem), " to ", $ strcompress(t+fbi+fei+1,/rem) if (JPL) then begin if (foc) then begin tdat[t+fbi:t+fbi+fei,*,0]=(ali) ? $ transpose(alog10(abs(complex(rdi[ch,*,0,0:fei,0],$ rdi[ch,*,1,0:fei,0]))+rds)):$ transpose(10.*alog10((abs(complex(rdi[ch,*,0,0:fei,0],$ rdi[ch,*,1,0:fei,0])))^2)) endif else begin tdat[t+fbi:t+fbi+fei,*,0]=(ali) ? $ transpose(alog10(abs(rdi[ch,*,0:fei,0])+rds)) : $ ; Ali scaling transpose(10.*alog10((abs(rdi[ch,*,0:fei,0]))^2)) ; power scaling ; 2009/11/18 NEP: not sure if squaring on previous line should be done endelse endif else begin tdat[t+fbi:t+fbi+fei,*,0]=(ali) ? $ alog10(((RPB)?transpose(rdi[ch,*,fbi:fbi+fei,0]):$ rdi[ch,fbi:fbi+fei,*,0])+rds):$ ; Ali scaling ( (pow) ? 10.*alog10(abs(((RPB)?transpose(rdi[ch,*,fbi:fbi+fei,0]): $ rdi[ch,fbi:fbi+fei,*,0]))) : $ ; power scaling ((UTS)?reverse(((RPB)?transpose(rdi[ch,*,fbi:fbi+fei,0]):$ rdi[ch,fbi:fbi+fei,*,0]),3) : $ ; no scaling ((RPB)?transpose(rdi[ch,*,fbi:fbi+fei,0]):$ rdi[ch,fbi:fbi+fei,*,0]))) endelse endif else if diag then $ print, " skipping frames ", strcompress(t+fbi+1,/rem), " to ", $ strcompress(t+fbi+fei+1,/rem) endfor rdat=tdat[t:t+nf-1,*,0] ; Check for non-finite values and warn if found, then null them out wnfn=where(finite(rdat) eq 0,wcount) if wcount gt 0 then begin scount=strcompress(wcount,/rem) print,"WARNING: found "+scount+" non-finite values, set to null value" rdat[wnfn]=null endif if diag then begin wrnn=where(rdat ne null) scaled=(pow or ali)?' scaled':'' print," min"+scaled+" value: ",min(rdat[wrnn]) print," max"+scaled+" value: ",max(rdat[wrnn]) print," mean"+scaled+" value: ",mean(rdat[wrnn]) endif ftn=t+1 ; first sequential trace number in file if nowrite eq 0 and diag then print,"Writing traces to "+ofile+":" for f=0l,nf-1 do begin ; loop over frames if nowrite eq 0 and t ge tro then begin ; CREATE SEGY trace header fields tt=t-tro c[tt]=intarr(ltr/2) c[0:lth/2-1,tt]=thead ; write long-word trace headers ; byte offset description ;d[0,tt]=[tt+1] ; 1-4 trace sequence number ;d[2,tt]=[tt+1] ; 9-12 field record number ;d[3,tt]=[tt+1] ; 13-16 field trace number ;d[5,tt]=[tt+1] ; 21-24 common depth point ensemble number ; 2014 Jan 14: use actual observation sequence numbers d[0,tt]=[t+1] ; 1-4 trace sequence number d[2,tt]=[t+1] ; 9-12 field record number d[3,tt]=[t+1] ; 13-16 field trace number d[5,tt]=[t+1] ; 21-24 common depth point ensemble number e[lth/4:ltr/4-1,tt]=rdat[f,*] ; write frame to trace file ; show output file size every 2000 frames if diag and $ (tt eq 0 or tt eq ntr-1 or (tt+1) mod 2000 eq 0) then begin fio=file_info(ofile) csiz=string(fio.size,format='(I12)') ;print,string(f+1-tro,format='(I7)')+" of "+$ print,string(tt+1,format='(I7)')+" of "+$ string(ntr,format='(I7)')+" traces, "+csiz+oss endif endif t=t+1l if t eq tro+ntr then break endfor ; end frame loop free_lun,ilun ; code for creating latlon file ; convert location file to sequential trace number file (calls orb2stn.pro) if nostn eq 0 then $ orb2stn,ob[n],nf*((UTS)?8:1),ostfn,ftn=ftn,wti=wti,lineid=lineid,$ diag=odiag,xy=xy,proj=proj,clat=clat,clon=clon,radii=radii if t eq tro+ntr then break endfor ; end input files loop rdat=0 ; remove rdat array from memory if nowrite eq 0 then free_lun,olun ; close output file if n_params() eq 2 then begin ; if called for, populate odat odat=tdat[tro:ntr+tro-1,*,0] ; subset for return ; Check for non-finite values and null them out wnfn=where(finite(odat) eq 0,wcount) if wcount gt 0 then begin scount=strcompress(wcount,/rem) odat[wnfn]=null endif if diag then begin wonn=where(odat ne null) print,"min value returned: ",min(odat[wonn]) print,"max value returned: ",max(odat[wonn]) print,"mean value returned: ",mean(odat[wonn]) endif endif if display then begin ddat=tdat[tro:ntr+tro-1,dso:dso+dns-1,0] ; subset for display ; Check for non-finite values and null them out wnfn=where(finite(ddat) eq 0,wcount) wafn=where(finite(ddat)) if wcount gt 0 then begin scount=strcompress(wcount,/rem) ddat[wnfn]=min(ddat[wafn]) endif ddat=rotate(ddat,7) ; flip time axis for display if !d.name eq 'X' then begin device,decomposed=0 ; allow proper colors win=(!d.window ge 0) ? !d.window : 0 window,win,xs=ntr/dsf,ys=dns/dst,title=oname ; set up window size endif loadct,ct,/silent ; color table (3 = RED TEMPERATURE) if dsf ne 1 or dst ne 1 then ddat=congrid(ddat,ntr/dsf,dns/dst) ; regrid ;if png then write_png,oname+'.png',rotate(bytscl(ddat),7) if png then write_png,oname+'.png',bytscl(ddat) if png and diag then print, 'Wrote '+oname+'.png' tvscl,ddat ; plot radargram if diag and (dns ne ns or dst ne 1 or dsf ne 1) then begin print,"min value for display: ",min(ddat) print,"max value for display: ",max(ddat) print,"mean value for display: ",mean(ddat) endif if debug then begin print,"min value in ",tfile,": ",min(tdat[0]) print,"max value in ",tfile,": ",max(tdat[0]) print,"mean value in ",tfile,": ",mean(tdat[0]) help endif endif free_lun,tlun ; close temporary file if debug eq 0 then file_delete, tfile if debug eq 0 and simpng eq 1 then file_delete,rf end