; procedure to convert CO-SHARPS orbit files to sequential trace number files ; for use in associating location to SeisWare SEGY radargrams. ; ; Written 2007 Feb 08 v00 Than Putzig ; Updated 2007 Oct 03 v00 Than Putzig, to handle JPL UPA/FPA format (*.etm) ; Updated 2007 Oct 03 v00 Than Putzig, to add wti (write trace interval) ; Updated 2009 Jul 30 v01 Than Putzig, to add /xy keyword (projected XYs) ; Updated 2011 Jul 21 v02 Than Putzig, to treat LINEID as long64 ; Updated 2014 Jan 10 v04 Than Putzig, to allow Smithsonian FPB PDS geom files ; Updated 2014 Jan 27 v05 Than Putzig, to allow Smithsonian FPB pixlatlon files ; and to allow user-specified projections ; Updated 2014 Nov 13 v05 Than Putzig, to add RADII keyword for PDS ; Updated 2015 Feb 18 v06 Than Putzig, to bypass time calcs for PDS files ; Updated 2015 Feb 27 v07 Than Putzig, to add /INTERP keyword ; Updated 2015 Mar 02 v07 Than Putzig, to strip lines beginning with # or $ ; (i.e., allow decoder orbit files). ; Updated 2015 Jun 24 v08 Than Putzig, fix for FPB v9.1 naming ; ; Orbit file (decoder, UPB,FPB) fields: ; UTC (from ET) Lat Lon Radius Vtang Vradial Pos X Pos Y Pos Z ; Vel X Vel Y Vel Z Roll Pitch Yaw HGAin HGAout SAPXin SAPXout ; SAMXin SAMXout SZA Mag_field Sun_dist ; ; Example geom (FPB PDS) format: ; Frame,UT_obs_time,Lat,Lon,Datum_rad,MRO_rad,Vradial,Vtang,SZA,Phase/1e16 ; 1,2006-12-09T01:00:56.075, 87.2214,346.8616,3378.220,3689.077, 1.9152,3398.7878, 74.99, 1.492 ; ; Example pixlatlon format: ; 72.915182, 232.98771 ; 72.922865, 232.98290 ; 72.930549, 232.97809 ; ; Example etm (UPA, FPA, QDA) format (field numbers): ; 1 2 3 4 5 6 7 8 9 10 11 12 13 14 ; 0 241586138.3946 316.7760 116.9090 79.9063 81.9087 3374.3550 -0.0045 3.3845 0.3093 2045.8018 319.0305 2007-08-28-T15:14:33.2118 2007-08-28-T15:14:33.2118 ; ; According to Ali, etm fields are (* = used here): ; 1. Frame Index ; 2. Epoch (J2000 Seconds past 2000) ; 3. Altitude above the ellipsoid (r1=3396, r2=3376 km) ; 4. Solar zenith angle ; 5. Latitude (deg) * ; 6. Longitude (deg) * ; 7. Nadir MOLA radius ; 8. vrad ; 9. vtan ; 10. along_track_travel (km) ; 11. tshift*1.0e6 (this is the timeshift applied to the frame to align with ; the ellipsoid. This includes a shift by half of the 3600 sample window ; sampled at 80/3 MHz. ; 12. S/C height ; 13. UTC time * ; 14. Duplicate of 13 (this is a spare field) ; pro orb2stn,infile,nframes,fileout,outfile=outfile,ftn=ftn,stn=stn,wti=wti,$ lineid=lineid,proj=proj,clat=clat,clon=clon,xy=xy,radii=radii,diag=diag,$ interp=interp,help=help version='08' pron='ORB2STN' prom='['+pron+']: ' debug=0 IF N_PARAMS() LT 1 OR KEYWORD_SET(HELP) THEN BEGIN PRINT,"" PRINT,"Version ", version PRINT,pron +",INFILE{,NFRAMES,OUTFILE,OUTFILE=OUTFILE,FTN=FTN,WTI=WTI,"+$ "STN=STN,LINEID=LINEID,PROJ=PROJ,CLAT=CLAT,CLON=CLON,/XY,/INTERP,"+$ "/RADII,/DIAG,/HELP}" IF KEYWORD_SET(HELP) THEN BEGIN PRINT," Parameter Description {default value}" PRINT," --------------------------------------------------------" PRINT," INFILE CO-SHARPS Orbit or PDS geom file (REQUIRED) {}" PRINT," e.g., 'OBS_219001000_1_Orbit_1.txt' (decoder)" PRINT," or 'UPB_219001000_1_01_Orbit_001.txt' (SI UPB)" PRINT," or 'FPB_219001000_1_01_pixlatlon_001.txt' (SI FPB)" PRINT," or 'QDA_219001000_1_Mode_001.dat.etm' (JPL)" PRINT," or 's_00219001_geom.tab' (SI PDS)" PRINT," or 'FPB_219001000_1_01_geom_001.txt' (SI FPB 9.1+)" PRINT," NFRAMES Number of frames in mode. IGNORED if INFILE " PRINT," is PDS or pixlatlon; otherwise, REQUIRED. {}" PRINT," OUTFILE Output file for sequential trace information. " PRINT," If not given, output is printed to screen. May " PRINT," be given as a parameter following NFRAMES. {}" PRINT," FTN First sequential trace number this mode. {1}" PRINT," WTI Write trace interval relative to input file.{1}" PRINT," STN Returned array of seq trace numbers. {}" PRINT," LINEID Specified or returned line id; default is to" PRINT," add mode number to product id {e.g., 219001001}" PRINT," /XY Output projected X and Y in meters. {}" PRINT," Associated parameters:" PRINT," PROJ Projection. {106 -> polar stereographic coords}" PRINT," CLAT Center latitude for projection. {90}" PRINT," CLON Center longitude for projection. {0}" PRINT," /INTERP Interpolate to unitary seq trace interval. " PRINT," Occurs in projection only. Ignored for unitary" PRINT," inputs (geom, pixlatlon). {}" PRINT," /RADII Output datum and MRO radii in km (PDS). {}" PRINT," /DIAG Print diagnostic messages." ENDIF PRINT, " /HELP Print parameter details." PRINT,"" RETURN STOP ENDIF fext=strmid(infile,2,3,/rev) ; input file extension JPL=(fext eq 'etm')?1:0 ; if true, input file is in JPL format gtyp=strmid(infile,7,4,/rev) ; extract 1st 4 of last 8 characters PDS=(strlowcase(gtyp) eq 'geom')?1:0 ; if geom, input file is SI PDS format gtyp2=strmid(infile,11,4,/rev) ; extract 1st 4 of last 12 characters PDS2=(strlowcase(gtyp2) eq 'geom')?1:0 ; if geom, input file is SI FPB v9.1+ PDF=(PDS or PDS2) xy=(keyword_set(xy)) ; projection keyword interp=(keyword_set(interp)) ; interpolation keyword radii=(keyword_set(radii)) ; radii keyword diag=(keyword_set(diag)) ; print diagnostics outfile=(n_elements(fileout) ne 0)?fileout:outfile prin=(n_elements(outfile) eq 0 and n_elements(fileout) eq 0); print output if n_elements(ftn) eq 0 then ftn=1 ; first sequential trace number if n_elements(wti) eq 0 then wti=1 ; write trace interval part=strsplit(file_basename(infile),'_.',/extract) ; parse input orbit file name npart=n_elements(part) if n_elements(lineid) eq 0 then begin ; output lineid obsn=long64(part[1])*((PDS)?1000:1) ; observation number from file name mode=(PDS or npart lt 7)?'A':part[5-JPL] ; mode from file name lineid=obsn+((mode ne 'A') ? fix(mode) : 0) endif PIX=(part[n_elements(part)-3] eq 'pixlatlon')?1:0 if diag then print,'Number of frames: ',nframes if PIX eq 0 and PDF eq 0 and n_elements(nframes) eq 0 then begin PRINT,"ERROR: NFRAMES is required if INFILE is neither PDS nor pixlatlon." RETURN STOP endif fsep=(PDF)?' -F,':'' ; awk field separator flat=(JPL)? '5':((PDF)?'3':((PIX)?'1':'2')) ; field for latitude flon=(JPL)? '6':((PDF)?'4':((PIX)?'2':'3')) ; field for longitude cawk="cat "+infile+ " | grep -v \^# | grep -v \^$ | awk" if PIX eq 0 and PDF eq 0 then begin ;ftim=(JPL)?'13':((PDF)?'2':'1') ; field for timestamp ftim=(JPL)?'13':'1' ; field for timestamp tawk=cawk+fsep+" '{ print $"+ftim+" }' | awk -FT '{ print $2 }' | awk -F:" spawn,tawk+" '{ print $1 }'",hou ; extract hour spawn,tawk+" '{ print $2 }'",min ; extract minute spawn,tawk+" '{ print $3 }'",sec ; extract second ; check for empty string in sec (possible with truncated records) estr=where(sec eq '',count) extr=0 if count ne 0 then begin if count eq 1 and estr[0] eq file_lines(infile)-1 then begin print, 'WARNING: Last line of '+infile print, ' is apparently truncated; will extrapolate last frame.' extr=1 endif else begin print, 'ERROR: ',count,' bad time-stamp values in ',infile RETURN STOP endelse endif ; check for asterisk string in sec (occurs in some GEOM.TAB files) estr=where(strmatch(sec,'*\**') eq 1,count) extr=0 if count ne 0 then begin print, 'ERROR: ',count,' bad time-stamp values in ',infile RETURN STOP endif ; double input values hou=double(hou) min=double(min) sec=double(sec) endif spawn,cawk+fsep+" '{ print $"+flat+" }' ",lat ; extract latitude spawn,cawk+fsep+" '{ print $"+flon+" }' ",lon ; extract longitude if radii then begin spawn,cawk+fsep+" '{ print $5 }' ",rdat ; extract datum radii spawn,cawk+fsep+" '{ print $6 }' ",rmro ; extract MRO radii endif nlo=n_elements(lat) ; number of lines in input file if diag then print,'Number of lines in input file: ',nlo lat=double(lat) lon=double(lon) if PIX eq 0 and PDF eq 0 then begin tim=hou*3600+min*60+sec ; compute time in seconds if extr then begin tim[nlo-1]=tim[nlo-2]+(tim[nlo-2]-tim[nlo-3]) lat[nlo-1]=lat[nlo-2]+(lat[nlo-2]-lat[nlo-3]) lon[nlo-1]=lon[nlo-2]+(lon[nlo-2]-lon[nlo-3]) endif ; compute seconds per frame ;;sif=tim[nlo-1]-tim[0]+86400.*(tim[0] gt tim[nlo-1]) ; seconds in file ;;spf=sif/(nframes) ; seconds per frame nlu=nlo/wti*wti ; num lines to use for spf calculation sif=tim[nlu-1]-tim[wti-1]+86400.*(tim[wti-1] gt tim[nlu-1]) ; seconds in file spf=sif/(nframes-1) ; seconds per frame if diag then print,'seconds in file: ',sif if diag then print,'seconds per frame: ',spf endif else nframes=nlo if diag then begin xylbl = (xy) ? ' X Y' : '' radlbl = (radii) ? 'Rad datum Rad MRO' : '' print,' Product ID Seq Trace # Latitude Longitude'+$ xylbl+radlbl+' Time (s)' endif interpolated=0 ; project XY from lat/lon if called for by xy keyword if xy then begin if n_elements(proj) eq 0 then proj=106 ; GCTP polar stereographic projection smja=3396000 ; semimajor axis for projection smna=3376000 ; semiminor axis for projection feas=2000000 ; false easting fnor=2000000 ; false northing if n_elements(clon) eq 0 then clon=0 ; center longitude for projection if n_elements(clat) eq 0 then clat=90 ; center latitude for projection mstr=map_proj_init(proj,semimajor_axis=smja,semiminor_axis=smna,$ center_longitude=clon,center_latitude=clat,$ false_easting=feas,false_northing=fnor) mpxy=map_proj_forward(lon,lat,map_structure=mstr) mpx=mpxy[0,*] mpy=mpxy[1,*] if interp and nframes gt nlo then begin mpx=interpol(mpx,nframes) mpy=interpol(mpy,nframes) lola=map_proj_inverse(mpx,mpy,map_structure=mstr) lon=lola[0,*] lat=lola[1,*] wnl=where(lon lt 0.,count) if count gt 0 then lon[wnl]=lon[wnl]+360. if PIX eq 0 and PDF eq 0 then tim=interpol(tim,nframes) nlo=nframes interpolated=1 endif endif nls=nlo/wti ; number of lines in stn output file if diag then print,'Number of lines intended for output file: ',nls dtn=dblarr(nlo) ; sequential trace number (double) stn=lonarr(nls) ; sequential trace number (long) ; set format for output ofmt='F12.6,F12.6'+((xy)?',F12.1,F12.1':'')+((radii)?',F9.3,F9.3':'') dfmt='(I12,F12.3,'+ofmt+',F12.3)' ofmt='(I12,I12,'+ofmt+')' ; determine sequential trace number for each line in input file if prin eq 0 then openw,olun,outfile,/get_lun,/append for n=0L, nlo-1 do begin dtn[n]=(n eq 0 or PIX or PDF or interpolated) ? double(n) : $ dtn[n-1]+(tim[n]-tim[n-1]+86400.*(tim[n] lt tim[n-1]))/spf ;m=n/wti m=round(dtn[n])/wti if m*wti eq n and m lt nls then begin ;stn[m]=round(dtn[n])+ftn stn[m]=ftn+m prxy=(xy)?',mpx[n],mpy[n]':'' prr=(radii)?',rdat[n],rmro[n]':'' ptm=(diag)?',tim[n]':'' fmt=(diag)?dfmt:ofmt if diag or prin then $ x=execute('print,format=fmt,lineid,stn[m],lat[n],lon[n]'+$ prxy+prr+ptm) if prin eq 0 then $ x=execute('printf,format=ofmt,olun,lineid,stn[m],lat[n],lon[n]'+$ prxy+prr) endif endfor ;if diag then help if prin eq 0 then free_lun,olun end