This repository has been archived by the owner on Nov 7, 2020. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 9
/
plot_results.py
166 lines (156 loc) · 6.13 KB
/
plot_results.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
# very simple code to pick up all the RH results and make a plot.
# only getting rid of the biggest outliers using a median filter
# Kristine Larson May 2019
# June 30, 2019 added extension to results directory
# require you to have a minimum number of points to use for a daily average
# was previously using a default value
import argparse
import datetime
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
from datetime import date
#
# where the results are stored
xdir = os.environ['REFL_CODE']
# must input start and end year
parser = argparse.ArgumentParser()
parser.add_argument("station", help="station name", type=str)
parser.add_argument("year1", help="first year", type=int)
parser.add_argument("year2", help="end year", type=int)
parser.add_argument("medfilter", help="median filter", type=float)
parser.add_argument("ReqTracks", help="required number of tracks", type=int)
# optional inputs: filename to output daily RH results
parser.add_argument("-txtfile", "--txtfile", default='None', type=str, help="txtfile for output")
parser.add_argument("-noscreen", "--noscreen", default='None', type=str, help="toggle to not plot to screen")
parser.add_argument("-extension", "--extension", default='None', type=str, help="extension for solution names")
args = parser.parse_args()
station = args.station
year1= args.year1
year2= args.year2
medfilter= args.medfilter
txtfile = args.txtfile
ReqTracks = args.ReqTracks
noscreen = args.noscreen
if args.extension == 'None':
extension = ''
else:
extension = args.extension
# where the summary files will be written to
txtdir = xdir + '/Files'
if not os.path.exists(txtdir):
print('make an output directory', txtdir)
os.makedirs(txtdir)
# outliers limit, defined in meters
howBig = medfilter;
k=0
# added standard deviation 2020 feb 14, changed n=6
n=7
# now require it as an input
# you can change this - trying out 80 for now
#ReqTracks = 80
# putting the results in a np.array, year, doy, RH, Nvalues, month, day
tv = np.empty(shape=[0, n])
obstimes = []
medRH = []
meanRH = []
plt.figure()
#plt.subplot(211)
yearEnd = year2 + 1
year_list = np.arange(year1,yearEnd,1)
print('Years to examine: ',year_list)
for yr in year_list:
if extension == '':
direc = xdir + '/' + str(yr) + '/results/' + station + '/'
else:
direc = xdir + '/' + str(yr) + '/results/' + station + '/' + extension + '/'
print('looking at ', yr, direc)
# check to make sure the directory exists
if os.path.isdir(direc):
all_files = os.listdir(direc)
for f in all_files:
fname = direc + f
L = len(f)
# file names have 7 characters in them ...
if (L == 7):
# check that it is a file and not a directory and that it has something/anything in it
try:
a = np.loadtxt(fname,skiprows=3,comments='%').T
numlines = len(a)
if (len(a) > 0):
y = a[0] +a[1]/365.25
rh = a[2]
doy = int(np.mean(a[1]))
# change from doy to month and day in datetime
d = datetime.date(yr,1,1) + datetime.timedelta(doy-1)
medv = np.median(rh)
cc = (rh < (medv+howBig)) & (rh > (medv-howBig))
good =rh[cc]; goodT =y[cc]
# only save if there are some minimal number of values
if (len(good) > ReqTracks):
rh = good
obstimes.append(datetime.datetime(year=yr, month=d.month, day=d.day, hour=12, minute=0, second=0))
medRH =np.append(medRH, medv)
plt.plot(goodT, good,'.')
# store the meanRH after the outliers are removed using simple median filter
meanRHtoday = np.mean(good)
stdRHtoday = np.std(good)
meanRH =np.append(meanRH, meanRHtoday)
# add month and day just cause some people like that instead of doy
# added standard deviation feb14, 2020
newl = [yr, doy, meanRHtoday, len(rh), d.month, d.day, stdRHtoday]
tv = np.append(tv, [newl],axis=0)
k += 1
else:
print('not enough retrievals on ', yr, d.month, d.day, len(good))
except:
print('problem reading ', fname, ' so skipping it')
else:
print('that directory does not exist - so skipping')
plt.ylabel('Reflector Height (m)')
plt.title('GNSS station: ' + station)
plt.gca().invert_yaxis()
plt.grid()
if (noscreen == 'None'):
plt.show()
#
#plt.subplot(211)
#plt.figure()
fig,ax=plt.subplots()
ax.plot(obstimes,meanRH,'.')
fig.autofmt_xdate()
#plt.xlabel('date')
plt.ylabel('Reflector Height (m)')
# so you know when it was computed
today = str(date.today())
plt.title(station.upper() + ': Daily Mean Reflector Height, Computed ' + today)
plt.grid()
plt.gca().invert_yaxis()
# save the second plot to a png file
pltname = txtdir + '/' + station + '_RH.png'
plt.savefig(pltname)
print('png file goes to this ', pltname)
if (noscreen == 'None'):
plt.show()
if txtfile == 'None':
print('no txt output')
else:
print('yes to making an output file')
# sort the time tags
ii = np.argsort(obstimes)
# apply time tags to a new variable
ntv = tv[ii,:]
N,M = np.shape(ntv)
outfile = txtdir + '/' + txtfile
xxx = str(datetime.datetime.now().strftime("%Y-%m-%d %H:%M"))
print('output file goes to this ', outfile)
fout = open(outfile, 'w+')
# change comment value from # to %
fout.write("{0:28s} \n".format( '% calculated on ' + xxx ))
fout.write("% year doy RH numval month day RH-sigma\n")
fout.write("% year doy (m) (m)\n")
fout.write("% (1) (2) (3) (4) (5) (6) (7)\n")
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} \n".format(ntv[i,0], ntv[i,1], ntv[i,2],ntv[i,3],ntv[i,4],ntv[i,5],ntv[i,6]))
fout.close()