; Procedure to extract a profile from CO-SHARPS MOLA topography or radii files. ; Useful for importing as a horizon into SeisWare. ; Created 2007 Mar 06 Than Putzig ; Updated 2009 Apr 01 NEP for UPB and directory structure ; Updated 2011 Aug 01 NEP for optionally using mola radii file as input ; Updated 2011 Aug 02 NEP added decimate switch ; Updated 2011 Aug 18 NEP to output prod_id, counter p+1 (as "CDP") for SeisWare ; Updated 2011 Oct 15 NEP to parameterize prodir ; Updated 2012 Jan 20 NEP to parameterize processor type PROC (UPB, FPB, etc.) ; Updated 2012 May 24 NEP v02 to use /data/*/WUSHARPS/Archive directories ; and to prompt for missing input file ; Updated 2018 May 05 NEP v03 to use /data/f for UPB pro mprofile,product,profile,prodir=prodir,proc=proc,modes=modes,polar=polar,$ radii=radii,interpolate=interpolate,decimate=decimate,plot=plot,help=help,$ debug=debug pron='MPROFILE' prom='['+pron+']: ' version='03' IF N_PARAMS() EQ 0 OR KEYWORD_SET(HELP) THEN BEGIN PRINT,"" PRINT,"Version "+version PRINT,pron+",PRODUCT{,PROFILE,PROC=PROC,PRODIR=PRODIR,MODES=MODES,POLAR=POLAR,/INTERPOLATE,/DECIMATE,/RADII,/PLOT,/DEBUG,/HELP}" IF KEYWORD_SET(HELP) THEN BEGIN PRINT," Parameter Description {default value}" PRINT," --------------------------------------------------------" PRINT," PRODUCT binary map or COSHARPS product (REQUIRED) {}" PRINT," e.g., 'my_obs_topo.dat' or '219001000_1_01'" PRINT," PROFILE extracted profile {not returned}" PRINT," PROC processor (FPB,UPB) {UPB}" PRINT," PRODIR product dir {'/data/f/WUSHARPS/Archive/UPB_PROD'}" PRINT," MODES e.g., [5,6,7] {ALL}" PRINT," POLAR e.g., [0,1,1] {NONE}" PRINT," /INTERPOLATE 4x interp to match equatorial to polar {}" PRINT," defaulted to ON if POLAR has nonzero elements. " PRINT," /DECIMATE x/4 rebin to match polar to equatorial {}" PRINT," overrides /INTERPOLATE" PRINT," /RADII extract MOLA radii {MOLA topography}" PRINT," /PLOT displays profile {}" PRINT," /DEBUG print parameters for debugging {}" ENDIF PRINT, " /HELP prints parameter details" PRINT,"" RETURN STOP ENDIF ; example input/output file names ;ifn='UPB_176902000_1_02_topo_002.dat' ; input topographic file name ;ifn='UPB_176902000_1_02_mola_002.dat' ; input radii file name ;ofn='UPB_176902000_1_02_topo_002.txt' ; output topographic file name ;ofn='UPB_176902000_1_02_mola_002.txt' ; output radii file name ; set debug keyword if n_elements(debug) eq 0 then debug=0b ; check RADII keyword and set MTYPE to get either radii or topo input file RADII=(keyword_set(radii)) MTYPE=(RADII) ? '_mola' : '_topo' MTYPD=(RADII) ? 'radii' : 'topographic' ; set default products directory if n_elements(proc) eq 0 then proc='UPB' ; if n_elements(prodir) eq 0 then prodir='/data/t/'+proc+'_PROD' if n_elements(prodir) eq 0 then begin arcv=(proc eq 'UPB')?'f':((proc eq 'FPB')?'g':((proc eq 'QDA')?'b1':'e')) prodir='/data/'+arcv+'/WUSHARPS/Archive/'+proc+'_PROD/' endif ; process PRODUCT parameter sproduct=strcompress(product,/rem) ; Test whether PRODUCT is a file name. PIF=file_test(sproduct,/regular) if PIF then begin ifn=sproduct fbn=file_basename(sproduct,'.dat') ofn=fbn+'.txt' dirp=strsplit(fbn,"_",/EXTRACT) pid=dirp[1] endif else begin dirp=strsplit(sproduct,"_",/EXTRACT) pid=dirp[0] pdir=proc+"_PROD_"+string(format='(I05)',long64(pid)/10000000*100) idir=prodir+'/'+pdir+'/'+proc+'_'+sproduct+'/' ; input file directory prefix=proc+'_'+sproduct+MTYPE ; input file prefix nm=n_elements(modes) ; count modes if nm ne 0 then begin smodes=string(modes,format='(I03)') ; convert modes to string ifn=idir+prefix+'_'+smodes+'.dat' ; input topo or radii files ofn=prefix+strjoin('_'+smodes)+'.txt' ; output profile file endif else begin ifn=file_search(idir+prefix+'_???.dat') ; all-mode input files ofn=prefix+'_A.txt' ; all-mode output profile file endelse endelse nf=n_elements(ifn) ; number of input files for n=0l,nf-1 do if file_test(ifn[n]) eq 0 then begin PRINT,"No input products "+idir+prefix+"_???.dat found." title='Please select a MAPTRACK '+MTYPD+' file:' ifn[n]=dialog_pickfile(title=title,filter='*_'+MTYPE+'_*.dat') if ifn[0] eq '' then begin print, 'No '+MTYPD+' file selected; exiting.' RETURN STOP endif endif ; set interpolate flag if given by user or polar array has non-zero elements if n_elements(polar) gt 0 then tpolar=total(polar) else tpolar=0 interpolate=((keyword_set(interpolate)) or (tpolar gt 0)) ; if polar not specified per-mode, set it uniformly (to 1 if /polar or /interp) if n_elements(polar) ne nf then $ polar=bytarr(nf)+((keyword_set(polar)) or interpolate) decimate=(keyword_set(decimate)) if decimate then interpolate=0 ; set constants mrp=0.001953125d ; map res (polar, dpp) mre=0.0078125 ; map res (equatorial, dpp) mw=3.5d ; maximum image width (deg) openw,ou,ofn,/get_lun ; open output file for n=0l,nf-1 do begin ; loop over input files openr,iu,ifn[n],/get_lun ; open input file nx=long(mw/(polar[n]?mrp:mre)) ; cross-track width in pixels spawn,"ls -l "+ifn[n]+" | awk '{ print $5 }'",fs ; find file size in bytes if RADII then begin na=long(fs[0]/nx/4) ; along-track length in pixels md=fltarr(nx,na) ; initialize MOLA data array ctp=fltarr(na) ; initialize profile array endif else begin fs=double(fs)/2 ; 16-bit-integer file size na=long(fs[0]/nx) ; along-track length in pixels md=intarr(nx,na) ; initialize MOLA data array ctp=intarr(na) ; initialize profile array endelse a=assoc(iu,md) ; associate data to variable md=a[0] ; copy data to array for m=0l,na-1 do ctp[m]=md[nx/2,m] ; extract center-track profile ; interpolate equatorial x4 to match polar if interpolate and polar[n] eq 0 then $ ctp=[long(interpol(ctp,4*(na-1)-1)),intarr(3)+ctp[na-1]] ; rebin polar /4 to match equatorial if decimate and polar[n] eq 1 then ctp=rebin(ctp[0:na/4*4-1],na/4) if debug then help,ctp profile=(n eq 0) ? ctp : [profile,ctp] ; append to profile free_lun, iu endfor ; print counter+1 and profile values to file for p=0l,n_elements(profile)-1 do printf,ou,pid,p+1,profile[p] if debug or keyword_set(plot) then begin ; diagnostic plot ytitle=(radii)?'radius (km)':'elevation (m)' scale=(radii)?1000:1 plot,profile/scale,xtitle='frame',ytitle=ytitle endif if debug then help ; diagnostic message if debug then print,"POLAR is ",polar ; diagnostic message free_lun, ou end