Skip to content

Commit

Permalink
updates to daily_avg (Hortho), quickplt (ydoy)
Browse files Browse the repository at this point in the history
  • Loading branch information
@kristinemlarson
kristinemlarson committed Mar 10, 2024
1 parent f7a27e4 commit 7885fd3
Show file tree
Hide file tree
Showing 6 changed files with 48 additions and 23 deletions.
34 changes: 22 additions & 12 deletions gnssrefl/daily_avg.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@

# my code
import gnssrefl.gps as g
import gnssrefl.sd_libs as sd
#

def fbias_daily_avg(station):
Expand Down Expand Up @@ -524,24 +525,28 @@ def daily_avg_stat_plots(obstimes,meanRH,meanAmp, station,txtdir,tv,ngps,nglo,ng
print('Number of values used in average RH file saved as: ', pltname)


def write_out_RH_file(obstimes,tv,outfile,csvformat):
def write_out_RH_file(obstimes,tv,outfile,csvformat,station,extension):
"""
write out the daily average RH values
Parameters
----------
obstimes : datetime object
time of observation
tv : numpy array
content of a LSP results file
outfile : string
full name of output file
csvformat : boolean
true if you want csv format output
station : str
4 ch station name
extension : str, optional
analysis extension name
"""
# ignore extension for now
Hortho = sd.find_ortho_height(station,extension)
# sort the time tags
nr,nc = tv.shape
# nothing to write
Expand All @@ -557,20 +562,25 @@ def write_out_RH_file(obstimes,tv,outfile,csvformat):
# header of a sorts
# change comment value from # to %
# 2021 november 8 added another column that has mean amplitude
fout.write("{0:28s} \n".format( '% calculated on ' + xxx ))
# 2024 march 10 added a column with orthometric height
fout.write("{0:48s} \n".format( '% station:' + station + ' calculated on ' + xxx ))
fout.write("% Daily average reflector heights (RH) calculated using the gnssrefl software \n")
fout.write("% Nominally you should assume this average is associated with 12:00 UTC hours \n")
fout.write("% year doy RH numval month day RH-sigma RH-amp\n")
fout.write("% (m) (m) (v/v)\n")
fout.write("% (1) (2) (3) (4) (5) (6) (7) (8) \n")
fout.write("% year doy RH numval month day RH-sigma RH-amp Hortho-RH\n")
fout.write("% (m) (m) (v/v) (m) \n")
fout.write("% (1) (2) (3) (4) (5) (6) (7) (8) (9) \n")

# these print statements could be much smarter
if csvformat:
for i in np.arange(0,N,1):
fout.write(" {0:4.0f}, {1:3.0f},{2:7.3f},{3:3.0f},{4:4.0f},{5:4.0f},{6:7.3f},{7:6.2f} \n".format(ntv[i,0],
ntv[i,1], ntv[i,2],ntv[i,3],ntv[i,4],ntv[i,5],ntv[i,6],ntv[i,7]))
dv = Hortho - ntv[i,2]
fout.write(" {0:4.0f}, {1:3.0f},{2:7.3f},{3:3.0f},{4:4.0f},{5:4.0f},{6:7.3f},{7:6.2f},{8:7.3f} \n".format(ntv[i,0],
ntv[i,1], ntv[i,2],ntv[i,3],ntv[i,4],ntv[i,5],ntv[i,6],ntv[i,7],dv ))
else:
for i in np.arange(0,N,1):
fout.write(" {0:4.0f} {1:3.0f} {2:7.3f} {3:3.0f} {4:4.0f} {5:4.0f} {6:7.3f} {7:6.2f} \n".format(ntv[i,0],
ntv[i,1], ntv[i,2],ntv[i,3],ntv[i,4],ntv[i,5],ntv[i,6], ntv[i,7]))
dv = Hortho - ntv[i,2]
fout.write(" {0:4.0f} {1:3.0f} {2:7.3f} {3:3.0f} {4:4.0f} {5:4.0f} {6:7.3f} {7:6.2f} {8:7.3f} \n".format(ntv[i,0],
ntv[i,1], ntv[i,2],ntv[i,3],ntv[i,4],ntv[i,5],ntv[i,6], ntv[i,7],dv))
fout.close()


Expand Down
2 changes: 1 addition & 1 deletion gnssrefl/daily_avg_cl.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,7 @@ def daily_avg(station: str , medfilter: float, ReqTracks: int, txtfile: str = No
outfile = txtdir + '/' + txtfile

if (nr > 0):
da.write_out_RH_file(obstimes, tv, outfile, csv)
da.write_out_RH_file(obstimes, tv, outfile, csv,station,extension)


def main():
Expand Down
9 changes: 5 additions & 4 deletions gnssrefl/gps.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1982,8 +1982,9 @@ def glonass_channels(f,prn):

def open_outputfile(station,year,doy,extension):
"""
opens output file in
REFL_CODE/year/results/station/extension directory
opens an output file in
$REFL_CODE/year/results/station/extension directory
for lomb scargle periodogram results
Parameters
----------
Expand Down Expand Up @@ -2018,9 +2019,9 @@ def open_outputfile(station,year,doy,extension):
# changed to a function
filepath1,fexit = LSPresult_name(station,year,doy,extension)
#print('Output will go to:', filepath1)
versionNumber = version('gnssrefl')
versionNumber = 'v' + str(version('gnssrefl'))
#versionNumber = 'working-on-it'
tem = '% gnssrefl, https://github.com/kristinemlarson, ' + str(versionNumber) + ' \n'
tem = '% station ' + station + ' https://github.com/kristinemlarson/gnssrefl ' + versionNumber + '\n'
try:
fout=open(filepath1,'w+')
# put a header in the output file
Expand Down
2 changes: 1 addition & 1 deletion gnssrefl/quicklib.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ def trans_time(tvd, ymd, ymdhm, convert_mjd, ydoy ,xcol,ycol,utc_offset):
yval = []

nr,nc = tvd.shape
print('rows and columns ', nr,nc)
#print('rows and columns ', nr,nc)
if (ycol+1 > nc):
print('You asked to plot column', ycol+1, ' and that column does not exist in the file')
sys.exit()
Expand Down
3 changes: 3 additions & 0 deletions gnssrefl/quickplt.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -223,6 +223,9 @@ def run_quickplt (filename: str, xcol: int, ycol: int, errorcol: int=None, mjd:
if ydoy:
print('Making obstimes for ydoy x-axis')
tval = g.ydoy2datetime(tvd[:,0], tvd[:,1])
if ydoy & secondFile:
print('Making obstimes for second ydoy x-axis')
tval2 = g.ydoy2datetime(tvd2[:,0], tvd2[:,1])

fig,ax=plt.subplots()

Expand Down
21 changes: 16 additions & 5 deletions gnssrefl/sd_libs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,9 @@ def write_spline_output(iyear, th, spline, delta_out, station, txtdir,Hortho):
$REFL_CODE/Files/ssss/ssss_2023_spline_out.txt
I do not think this is used anymore. It has been consolidated with the plot code.
Parameters
----------
iyear : int
Expand Down Expand Up @@ -72,13 +75,17 @@ def write_spline_output(iyear, th, spline, delta_out, station, txtdir,Hortho):
# but you only want those values when we have data ....
splinefileout = txtdir + '/' + station + '_' + str(iyear) + '_spline_out_orig.txt'
print('Writing evenly sampled file to: ', splinefileout)
vn = ' gnssrefl v' + str(g.version('gnssrefl'))
fout = open(splinefileout,'w+')
fout.write('{0:1s} {1:30s} \n'.format('%','station: ' + station + vn))
fout.write('{0:1s} {1:30s} \n'.format('%','This is NOT observational data - be careful when interpreting it.'))
fout.write('{0:1s} {1:30s} \n'.format('%','If the data are not well represented by the spline functions, you will '))
fout.write('{0:1s} {1:30s} \n'.format('%','have a very poor representation of the data. I am also writing out station '))
fout.write('{0:1s} {1:30s} {2:8.3f} \n'.format('%','orthometric height minus RH, where Hortho (m) is ', Hortho ))
fout.write('{0:1s} {1:30s} \n'.format('%','This assumes RH is measured relative to the L1 phase center. '))
fout.write('{0:1s} {1:30s} \n'.format('%','MJD, RH(m), YY,MM,DD,HH,MM,SS, quasi-sea-level(m)'))
fout.write('{0:1s} {1:30s} \n'.format('%','MJD RH(m) YYYY MM DD HH MM SS quasi-sea-level(m)'))
fout.write('{0:1s} {1:30s} \n'.format('%','(1) (2) (3) (4) (5) (6) (7) (8) (9)'))


dtime = False
for i in range(0,N):
Expand All @@ -87,7 +94,7 @@ def write_spline_output(iyear, th, spline, delta_out, station, txtdir,Hortho):
utc= 24*(tplot[i] - doy)
bigt,yy,mm,dd,hh,mi,ss = g.ymd_hhmmss(iyear,doy,utc,dtime)
if (tplot[i] >= firstpoint) & (tplot[i] <= lastpoint):
fout.write('{0:15.7f} {1:10.3f} {2:4.0f} {3:2.0f} {4:2.0f} {5:2.0f} {6:2.0f} {7:2.0f} {8:10.3f} \n'.format(
fout.write('{0:15.7f} {1:10.3f} {2:4.0f} {3:3.0f} {4:3.0f} {5:3.0f} {6:3.0f} {7:3.0f} {8:10.3f} \n'.format(
modjul, spline_even[i], yy,mm,dd,hh,mi,ss, Hortho-spline_even[i]))
fout.close()

Expand Down Expand Up @@ -991,12 +998,16 @@ def RH_ortho_plot2( station, H0, year, txtdir, fs, time_rh, rh, gap_min_val,th,
splinefileout = txtdir + '/' + station + '_spline_out.txt'
print('Writing evenly sampled file to: ', splinefileout)
fout = open(splinefileout,'w+')
vn = station + ' gnssrefl v' + str(g.version('gnssrefl'))
fout.write('{0:1s} {1:30s} \n'.format('%','station ' + vn))
fout.write('{0:1s} {1:30s} \n'.format('%','This is NOT observational data - be careful when interpreting it.'))
fout.write('{0:1s} {1:30s} \n'.format('%','If the data are not well represented by the spline functions, you will '))
fout.write('{0:1s} {1:30s} \n'.format('%','have a very poor representation of the data. I am also writing out station '))
fout.write('{0:1s} {1:30s} {2:8.3f} \n'.format('%','orthometric height minus RH, where Hortho (m) is ', H0 ))
fout.write('{0:1s} {1:30s} \n'.format('%','This assumes RH is measured relative to the L1 phase center. '))
fout.write('{0:1s} {1:30s} \n'.format('%','MJD, RH(m), YY,MM,DD,HH,MM,SS, quasi-sea-level(m)'))
#fout.write('{0:1s} {1:30s} \n'.format('%','This assumes RH is measured relative to the L1 phase center. '))
#fout.write('{0:1s} {1:30s} \n'.format('%','MJD, RH(m), YY,MM,DD,HH,MM,SS, quasi-sea-level(m)'))
fout.write('{0:1s} {1:30s} \n'.format('%','MJD RH(m) YYYY MM DD HH MM SS quasi-sea-level(m)'))
fout.write('{0:1s} {1:30s} \n'.format('%','(1) (2) (3) (4) (5) (6) (7) (8) (9)'))


# difference function to find time between all rh measurements
Expand Down Expand Up @@ -1039,7 +1050,7 @@ def RH_ortho_plot2( station, H0, year, txtdir, fs, time_rh, rh, gap_min_val,th,
for i in range(0,N_new):
if not np.isnan(spline_new[i]):
rhout = spline_new[i]
fout.write('{0:15.7f} {1:10.3f} {2:4.0f} {3:2.0f} {4:2.0f} {5:2.0f} {6:2.0f} {7:2.0f} {8:10.3f} \n'.format(
fout.write('{0:15.7f} {1:10.3f} {2:4.0f} {3:3.0f} {4:3.0f} {5:3.0f} {6:3.0f} {7:3.0f} {8:10.3f} \n'.format(
mjd_new[i], rhout, theyear[i], xm[i], xd[i], xh[i], xmin[i], xs[i], H0-rhout))


Expand Down

0 comments on commit 7885fd3

Please sign in to comment.