7 from os.path
import expanduser
8 from scipy.interpolate
import interp1d
14 from log_reader
import LogReader
17 from pylib.data_sets
import *
18 from pylib.pose
import *
20 from pylib.ISToolsDataSorted
import refLla, getTimeFromTowMs, getTimeFromTow, setGpsWeek, getTimeFromGTime
22 RAD2DEG = 180.0 / np.pi
23 DEG2RAD = np.pi / 180.0
35 def load(self, directory, serials=['ALL']):
44 raise ValueError(
"No devices found in log")
48 refIdx = self.
serials.index(10101)
54 self.
rtk =
'Rov' in str(self.
data[0, DID_DEV_INFO][
'addInfo'][-1])
55 self.
navMode = (self.
data[0, DID_INS_2][
'insStatus'][-1] & 0x1000) == 0x1000
65 while dev_id >= len(self.
data):
66 self.
data.append([[]
for i
in range(NUM_DIDS+1)])
67 self.
data[dev_id][did] = arr
72 if self.
data[dev_id][did] == []:
73 self.
data[dev_id][did] = [[]
for i
in range(6)]
74 if dev_id < len(self.
data)
and did < len(self.
data[dev_id])
and 6 <= len(self.
data[dev_id][did]):
75 if msg_type == eRawDataType.raw_data_type_observation.value:
76 self.
data[dev_id][did][0].append(arr)
77 elif msg_type == eRawDataType.raw_data_type_ephemeris.value:
78 self.
data[dev_id][did][1] = arr
79 elif msg_type == eRawDataType.raw_data_type_glonass_ephemeris.value:
80 self.
data[dev_id][did][2] = arr
81 elif msg_type == eRawDataType.raw_data_type_sbas.value:
82 self.
data[dev_id][did][3] = arr
83 elif msg_type == eRawDataType.raw_data_type_ionosphere_model_utc_alm.value:
84 self.
data[dev_id][did][4] = arr
85 elif msg_type == eRawDataType.raw_data_type_base_station_antenna_position.value:
86 self.
data[dev_id][did][5] = arr
91 GPS_start_Time = datetime.datetime.strptime(
'6/Jan/1980',
"%d/%b/%Y")
94 if np.size(self.
data[0][DID_INS_2])==0
and np.size(self.
data[0][DID_INS_1])!=0:
95 self.
data[0][DID_INS_2].resize(np.size(self.
data[0][DID_INS_1]))
96 self.
data[0][DID_INS_2][
'timeOfWeek'] = self.
data[0][DID_INS_1][
'timeOfWeek']
97 self.
data[0][DID_INS_2][
'week'] = self.
data[0][DID_INS_1][
'week']
98 self.
data[0][DID_INS_2][
'insStatus'] = self.
data[0][DID_INS_1][
'insStatus']
99 self.
data[0][DID_INS_2][
'hdwStatus'] = self.
data[0][DID_INS_1][
'hdwStatus']
101 self.
data[0][DID_INS_2][
'uvw'] = self.
data[0][DID_INS_1][
'uvw']
102 self.
data[0][DID_INS_2][
'lla'] = self.
data[0][DID_INS_1][
'lla']
104 week_time = GPS_start_Time + (datetime.timedelta(weeks=int(self.
data[0][DID_INS_2][
'week'][-1])))
106 for d, dev
in enumerate(self.
data):
107 if len(dev[DID_INS_2]) == 0:
108 print(
"\033[93m" +
"missing DID_INS_2 data: removing device\033[0m")
112 if len(dev[DID_DEV_INFO]) == 0:
113 print(
"\033[93m" +
"missing DID_DEV_INFO data: making some up\033[0m")
114 self.
data[d][DID_DEV_INFO] = np.resize(self.
data[d][DID_DEV_INFO], 1)
115 self.
data[d][DID_DEV_INFO][
'serialNumber'][0] = d
117 for d, dev
in enumerate(self.
data):
119 if isinstance(dev[did], list):
121 for field
in [
'towMs',
'timeOfWeekMs',
'tow',
'timeOfWeek']:
122 if field
in dev[did].dtype.names:
123 if (np.diff(dev[did][field].astype(np.int64)) < 0).any():
124 idx = np.argmin(np.diff(dev[did][field].astype(np.int64)))
126 t1 = week_time + datetime.timedelta(milliseconds=int(dev[did][field][idx]))
127 t2 = week_time + datetime.timedelta(milliseconds=int(dev[did][field][idx+1]))
129 t1 = week_time + datetime.timedelta(seconds=int(dev[did][field][idx]))
130 t2 = week_time + datetime.timedelta(seconds=int(dev[did][field][idx+1]))
131 print(
"\033[93m" +
"Time went backwards in " + did_name_lookup[did] +
r"!!!, " +
132 " Time went from " + str(t1) +
" to " + str(t2) +
". removing all data " 133 + (
"before" if idx < len(dev[did][field])/2.0
else "after") +
"\033[0m")
134 if idx < len(dev[did][field])/2.0:
135 self.
data[d][did] = dev[did][idx+1:]
137 self.
data[d][did] = dev[did][:idx]
138 ms_multiplier = 1000.0
if 'Ms' in field
else 1.0
139 if (np.diff(dev[did][field]) > 3600 * ms_multiplier).any():
140 print(
"\033[93m" +
"greater than 1 minute gap in " + did_name_lookup[did]
141 +
" data, assuming GPS fix was acquired during data set, and chopping data"+
"\033[0m")
142 idx = np.argmax(np.diff(dev[did][field])) + did
143 self.
data[d][did] = dev[did][idx:]
147 print(
"\nComputing RMS Accuracies: (%d devices)" % (self.
numDev))
150 data = [np.hstack((self.
data[i, DID_INS_2][
'timeOfWeek'][:,
None],
151 self.
data[i, DID_INS_2][
'lla'],
152 self.
data[i, DID_INS_2][
'uvw'],
157 if (np.diff(data[dev][:, 0]) > 10.0).any():
158 print(
"\033[93m" +
"large gaps in data for dev" + str(dev)
159 +
"chopping off data before gap".format(dev) +
"\033[0m")
160 idx = np.argmax(np.diff(data[dev][:, 0])) + 1
161 data[dev] = data[dev][idx:, :]
169 time_of_fix_ms = [self.
data[dev, DID_GPS1_POS][
'timeOfWeekMs'][np.argmax(self.
data[dev, DID_GPS1_POS][
'status'] & 0x08000000)] / 1000.0
for dev
in range(self.
numDev)]
182 data[i] = data[i][data[i][:, 0] > self.
min_time]
183 data[i] = data[i][data[i][:, 0] < self.
max_time]
189 fi = interp1d(data[i][:, 0], data[i][:, 1:].T, kind=
'cubic', fill_value=
'extrapolate',
191 data[i] = np.hstack((t[:,
None], fi(t).T))
194 data[i][:, 7:] /=
norm(data[i][:, 7:], axis=1)[:,
None]
197 data = np.array(data)
200 refLla = data[0, int(round(len(t) / 2.0)), 1:4].copy()
202 data[i, :, 1:4] =
lla2ned(refLla, data[i, :, 1:4])
209 means = np.empty((len(self.
stateArray[0]), 10))
210 means[:, :6] = np.mean(self.
stateArray[:, :, 1:7], axis=0)
240 self.
RMS[:,:6] = np.sqrt(np.mean(np.square(self.
stateArray[:, :, 1:7] - self.
truth[:,0:6]), axis=1))
241 self.
RMS[:,6:] = np.sqrt(np.mean(np.square(self.
att_error[:, :, :]), axis=1))
259 filename = os.path.join(self.
directory,
'RMS_report_new_logger.txt')
260 thresholds = np.array([0.35, 0.35, 0.8,
266 thresholds[:6] = np.inf
273 thresholds[6:] *= DEG2RAD
279 f = open(filename,
'w')
280 f.write(
'***** Performance Analysis Report - %s *****\n' % (self.
directory))
282 f.write(
'Directory: %s\n' % (self.
directory))
286 if self.
refINS: mode +=
" With NovaTel Reference" 293 '--------------------------------------------------- RMS Accuracy -------------------------------------------\n')
294 line = line +
'UVW[ (m/s) (m/s) (m/s) ], NED[ (m) (m) (m) ],' 296 f.write(
'-------------- RMS Accuracy --------------\n')
297 line = line +
' Att [ (deg) (deg) (deg) ]\n' 300 for n, dev
in enumerate(uINS_device_idx):
301 devInfo = self.
data[dev,DID_DEV_INFO][0]
302 line =
'%2d SN%d ' % (n, devInfo[
'serialNumber'])
304 line = line +
'[ %6.4f %6.4f %6.4f ], ' % (self.
RMS[n, 3], self.
RMS[n, 4], self.
RMS[n, 5])
305 line = line +
'[ %6.4f %6.4f %6.4f ], ' % (self.
RMS[n, 0], self.
RMS[n, 1], self.
RMS[n, 2])
306 line = line +
'[ %6.4f %6.4f %6.4f ]\n' % (
313 '------------------------------------------------------------------------------------------------------------\n')
317 f.write(
'------------------------------------------\n')
318 line = line +
'[%7.4f %7.4f %7.4f ]\n' % (
324 line = line +
'[%7.4f %7.4f %7.4f ], ' % (thresholds[3], thresholds[4], thresholds[5])
325 line = line +
'[%7.4f %7.4f %7.4f ], ' % (thresholds[0], thresholds[1], thresholds[2])
326 line = line +
'[%7.4f %7.4f %7.4f ]\n' % (
327 thresholds[6] * RAD2DEG, thresholds[7] * RAD2DEG, thresholds[8] * RAD2DEG)
333 '------------------------------------------------------------------------------------------------------------\n')
337 f.write(
'------------------------------------------\n')
343 line = line +
'[ %s %s %s ], ' % (
345 line = line +
'[ %s %s %s ], ' % (
347 line = line +
'[ %s %s %s ]\n' % (
355 f.write(
'(' + mode +
')\n\n')
358 f.write(
'--------------- Angular Mounting Biases ----------------\n')
359 f.write(
'Device Euler Biases[ (deg) (deg) (deg) ]\n')
360 for n, dev
in enumerate(uINS_device_idx):
361 devInfo = self.
data[dev, DID_DEV_INFO][0]
362 f.write(
'%2d SN%d [ %7.4f %7.4f %7.4f ]\n' % (
367 f.write(
"----------------- Average Attitude ---------------------\n")
368 f.write(
"Dev: \t[ Roll\t\tPitch\t\tYaw ]\n")
372 f.write(
"%d\t%f\t%f\t%f\n" % (self.
serials[i], euler[0], euler[1], euler[2]))
376 '\n\n------------------------------------------- Device Info -------------------------------------------------\n')
377 for n, dev
in enumerate(uINS_device_idx):
378 devInfo = self.
data[dev, DID_DEV_INFO][0]
379 hver = devInfo[
'hardwareVer']
380 cver = devInfo[
'protocolVer']
381 fver = devInfo[
'firmwareVer']
382 buld = devInfo[
'buildNumber']
383 repo = devInfo[
'repoRevision']
384 date = devInfo[
'buildDate']
385 time = devInfo[
'buildTime']
386 addi = devInfo[
'addInfo']
388 '%2d SN%d HW: %d.%d.%d.%d FW: %d.%d.%d.%d build %d repo %d Proto: %d.%d.%d.%d Date: %04d-%02d-%02d %02d:%02d:%02d %s\n' % (
389 n, devInfo[
'serialNumber'],
390 hver[3], hver[2], hver[1], hver[0],
391 fver[3], fver[2], fver[1], fver[0], buld, repo,
392 cver[3], cver[2], cver[1], cver[0],
393 2000 + date[2], date[1], date[0],
394 time[3], time[2], time[1],
405 import matplotlib.pyplot
as plt
406 colors = [
'r', 'g', 'b', 'm']
409 plt.title(
"position error")
413 plt.plot(self.
stateArray[0,:,0], self.
truth[:, m], linewidth=2, color = colors[m])
415 plt.title(
"velocity error")
419 plt.plot(self.
stateArray[0,:,0], self.
truth[:, m+3], linewidth=2, color = colors[m])
421 plt.title(
"attitude")
425 plt.plot(self.
stateArray[0,:,0], self.
truth[:, m+6], linewidth=2, color = colors[m])
429 plt.subplot(3, 1, m +1)
435 filename = os.path.join(self.
directory,
'RMS_report_new_logger.txt')
436 if 'win' in sys.platform:
437 subprocess.Popen([
"notepad.exe", filename])
438 if 'linux' in sys.platform:
439 subprocess.Popen([
'gedit', filename])
442 if __name__ ==
'__main__':
443 np.set_printoptions(linewidth=200)
445 home = expanduser(
"~")
448 if len(sys.argv) >= 2:
449 directory = sys.argv[1]
452 if 'directory' not in locals():
453 print(
"First parameter must be directory!")
457 file = open(home +
"/Documents/Inertial_Sense/config.yaml",
'r') 458 config = yaml.load(file) 459 directory = config["directory"]
def did_callback(self, did, arr, dev_id)
bool init(test_data_t &t)
GeneratorWrapper< T > range(T const &start, T const &end, T const &step)
#define max(a, b)
Takes the maximal value of a and b.
def pass_fail(self, ratio)
def gps_raw_data_callback(self, did, arr, dev_id, msg_type)
def load(self, directory, serials=['ALL'])
def calcAttitudeError(self)
def euler2quatArray(euler)