source: ntrip/trunk/BNC/scripts/plot_xyz.py

Last change on this file was 10791, checked in by mervart, 3 weeks ago

BNC Multifrequency and PPPAR Client (initial version)
  • Property svn:executable set to *
File size: 8.1 KB
Line 
1#!/usr/bin/python3
2
3import sys
4import datetime
5import matplotlib.pyplot as plt
6import matplotlib.dates as mdates
7import matplotlib.ticker as mticker
8import numpy as np
9import pathlib as pth
10
11import utils
12
13sys.dont_write_bytecode = True
14
15# Data Class
16##########################################################################################
17class Data:
18 def __init__(self):
19 self.epo0 = None
20 self.epo = []
21 self.dNfix = []
22 self.dEfix = []
23 self.dUfix = []
24 self.dNflt = []
25 self.dEflt = []
26 self.dUflt = []
27 self.nFix = []
28 self.conv2D = {}
29
30 def append(self, epo, dN, dE, dU, isFix, nFix, flgReset):
31 self.epo.append(epo)
32 self.nFix.append(nFix)
33 if isFix:
34 self.dNfix.append(dN)
35 self.dEfix.append(dE)
36 self.dUfix.append(dU)
37 self.dNflt.append(np.nan)
38 self.dEflt.append(np.nan)
39 self.dUflt.append(np.nan)
40 else:
41 self.dNflt.append(dN)
42 self.dEflt.append(dE)
43 self.dUflt.append(dU)
44 self.dNfix.append(np.nan)
45 self.dEfix.append(np.nan)
46 self.dUfix.append(np.nan)
47
48 if flgReset or self.epo0 is None:
49 self.epo0 = epo
50
51 nSec = int((epo - self.epo0).total_seconds())
52 if nSec not in self.conv2D:
53 self.conv2D[nSec] = []
54 self.conv2D[nSec].append(np.sqrt(dN*dN + dE*dE))
55
56# Function to generate the plot
57##########################################################################################
58def run_plot_client(staName, data, title, timeStr, pngFile, statFile):
59 cm = 1/2.54
60 fig, (ax1, ax2, ax3, dummy, ax4) = plt.subplots(5, 1, figsize=(20*cm,20*cm),
61 gridspec_kw={'height_ratios': [1, 1, 1, 0.3, 1.5]})
62 dummy.axis("off")
63
64 ax1.set_title(title)
65 ax1.set_ylabel("North [m]")
66 ax1.tick_params(labelbottom=False)
67 ax1.plot(data.epo, data.dNflt, color = "crimson", linewidth = 1.0)
68 ax1.plot(data.epo, data.dNfix, color = "forestgreen", label = "fixed", linewidth = 1.0)
69
70 ax2.set_ylabel("East [m]")
71 ax2.tick_params(labelbottom=False)
72 ax2.plot(data.epo, data.dEflt, color = "crimson", label = "float", linewidth = 1.0)
73 ax2.plot(data.epo, data.dEfix, color = "forestgreen", linewidth = 1.0)
74
75 ax3.set_xlabel(timeStr)
76 ax3.set_ylabel("Height [m]")
77 ax3.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M'))
78 ax3.plot(data.epo, data.dUflt, color = "crimson", linewidth = 1.0)
79 ax3.plot(data.epo, data.dUfix, color = "forestgreen", linewidth = 1.0)
80
81 ax32 = ax3.twinx()
82 ax32.plot(data.epo, data.nFix, color = "lightskyblue", label = "% fix", linewidth = 0.4)
83 ax32.set_ylim(0,300)
84 ax32.set_yticks([0,50,100])
85 ax32.legend(loc="upper right")
86
87 for ax in [ax1, ax2, ax3]:
88 ax.set_ylim(-0.4, 0.4)
89 hh, ll = ax.get_legend_handles_labels()
90 if len(hh) > 0:
91 ax.legend(loc="upper right")
92 ax.grid(axis = 'y', which = 'both', color = 'grey', linestyle = '--', linewidth = 0.3)
93 ax.set_yticks([x / 10 for x in range(-4, 5, 2)])
94 ax.yaxis.set_minor_locator(mticker.AutoMinorLocator(2))
95
96 maxLen = 0
97 perc = {68: [], 95: []}
98 percMin = {10: {68 : 0.0, 95 : 0.0}, 20: {68 : 0.0, 95 : 0.0}, 30: {68 : 0.0, 95 : 0.0}}
99 keyMin = {10: {68 : 0.0, 95 : 0.0}, 20: {68 : 0.0, 95 : 0.0}, 30: {68 : 0.0, 95 : 0.0}}
100
101 for key in data.conv2D.keys():
102 if len(data.conv2D[key]) > maxLen:
103 maxLen = len(data.conv2D[key])
104
105 for pKey in (68, 95):
106 pp = np.percentile(data.conv2D[key], pKey)
107 perc[pKey].append(pp)
108 for minute in (10, 20, 30):
109 tMin = 60*minute
110 if key in range(60*minute-10, 60*minute+1):
111 if abs(key - tMin) < abs(keyMin[minute][pKey] - tMin):
112 keyMin[minute][pKey] = key
113 percMin[minute][pKey] = pp
114
115 print("%8s 68: %6.3f %6.3f %6.3f 95: %6.3f %6.3f %6.3f" %
116 (staName, percMin[10][68], percMin[20][68], percMin[30][68], percMin[10][95], percMin[20][95], percMin[30][95]))
117
118 for ii in range(0, maxLen):
119 dPos = []
120 for key in data.conv2D.keys():
121 if len(data.conv2D[key]) > ii:
122 dPos.append(data.conv2D[key][ii])
123 else:
124 dPos.append(np.nan)
125 ax4.plot(data.conv2D.keys(), dPos, linewidth = 0.5)
126
127 ax4.set_title("Convergence")
128 ax4.plot(data.conv2D.keys(), perc[68], linewidth = 1.0, color = "black", label = "68 %")
129 ax4.plot(data.conv2D.keys(), perc[95], linewidth = 2.0, color = "black", label = "95 %")
130 ax4.legend(loc="upper right")
131 ax4.set_ylim(0, 1.0)
132 ax4.set_xlabel("Seconds after Reset")
133 ax4.set_ylabel("Meters")
134 ax4.text( 0.1, 0.7,
135 "percentile after 10 min:\n68%% : %8.3f m\n95%% : %8.3f m" % (percMin[10][68], percMin[10][95]),
136 transform=ax4.transAxes)
137 ax4.text( 0.4, 0.5,
138 "percentile after 20 min:\n68%% : %8.3f m\n95%% : %8.3f m" % (percMin[20][68], percMin[20][95]),
139 transform=ax4.transAxes)
140 ax4.text( 0.7, 0.3,
141 "percentile after 30 min:\n68%% : %8.3f m\n95%% : %8.3f m" % (percMin[30][68], percMin[30][95]),
142 transform=ax4.transAxes)
143
144 ### plt.show()
145 plt.savefig(pngFile)
146 plt.close()
147
148 # Output Statistics
149 # -----------------
150 if statFile:
151 with open(statFile, 'w') as outStat:
152 print("%s %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f" %
153 (fileName,
154 percMin[10][68], percMin[10][95],
155 percMin[20][68], percMin[20][95],
156 percMin[30][68], percMin[30][95]), file = outStat)
157
158
159# Main Program
160##########################################################################################
161if __name__ == '__main__':
162 import sys
163 import argparse
164
165 parser = argparse.ArgumentParser()
166 parser.add_argument("crdFile")
167 parser.add_argument("fileName")
168 parser.add_argument("--pngFile", type=str)
169 parser.add_argument("--statFile", type=str)
170 parser.add_argument("--title", type=str)
171 args = parser.parse_args()
172
173 crdFile = args.crdFile
174 fileName = args.fileName
175 if args.pngFile is None:
176 pngFile = pth.Path(fileName).with_suffix(".png")
177 else :
178 pngFile = args.pngFile
179 if args.title is None:
180 title = fileName
181 else:
182 title = args.title
183 statFile = args.statFile
184
185 dateFmt = "%Y-%m-%d_%H:%M:%S.%f"
186 dfCrd = utils.readcrdfile(crdFile)
187 station = None
188
189 # Read Data
190 # ---------
191 data = Data()
192 flgReset = True
193 timeStr = None
194 with open(fileName, 'r') as inFile:
195 for line in inFile:
196
197 if line.find("RESET FILTER") == 0:
198 flgReset = True
199
200 elif line.find("X =") >= 0:
201 fields = line.split()
202
203 dateStr = fields[0]
204 staName = fields[1]
205 xx = float(fields[4])
206 yy = float(fields[9])
207 zz = float(fields[14])
208 isFix = (fields[32] == "fix")
209 if isFix:
210 numFix = int(fields[33])
211 else:
212 numFix = 0
213
214 if station is None:
215 station = utils.getsta(dfCrd, staName)
216 if station is None:
217 raise Exception("Station %s not found" % staName)
218 ell = utils.xyz2ell((xx, yy, zz))
219
220 xyz = (xx - station.xx,
221 yy - station.yy,
222 zz - station.zz)
223
224 neu = utils.xyz2neu(ell, xyz)
225
226 time = datetime.datetime.strptime(dateStr, dateFmt)
227 data.append(time, neu[0], neu[1], neu[2], isFix, numFix, flgReset)
228 flgReset = False
229
230 if timeStr is None:
231 timeStr = time.strftime("%Y-%m-%d")
232
233 run_plot_client(staName, data, title, timeStr, pngFile, statFile)
Note: See TracBrowser for help on using the repository browser.