In [35]:
%matplotlib inline
import matplotlib.pyplot as plt
def fax(w=15, h=10, n=111):
fig = plt.figure(figsize=(w,h))
ax = fig.add_subplot(n)
ax.grid()
return fig, ax
In [6]:
import numpy as np
import pandas as pd
from datetime import datetime
from calendar import timegm
from pyfbf import Workspace
In [51]:
kts2mps = lambda x: 0.514444 * float(x) if x!='' else np.NAN
ft2m = lambda x: 0.3048 * float(x) if x!='' else np.NAN
def iso2aircraftTime(s):
if s=='':
return np.NAN
dt = datetime.strptime(s, '%Y-%m-%dT%H:%M:%S.%f')
midnight = datetime(dt.year, dt.month, dt.day)
hours = (dt - midnight).total_seconds() / 3600.0
return hours
In [52]:
fn_iwg1 = '/data/rdr/shis/dock/sh170411/IWG1.csv'
columns = [ "iwg", "aircraftTime", "instrumentLatitude", "instrumentLongitude", "_",
"Altitude", "barometricAltitude", "radarAltitude", "groundSpeed", "_",
"aircraftIndicatedAirSpeed", "aircraftMachNumber", "aircraftVerticalSpeed",
"Heading", "aircraftTrackAngle", "driftAngle",
"Pitch", "Roll", "aircraftSideSlip", "aircraftAttackAngle",
"staticAirTemp", "dewpointTemp", "totalAirTemp", "aircraftExtPressure", "_", "_",
"windSpeed", "windDirection", "_", "solarZenithAngle",
"sunAircraftElevation", "sunGroundAzimuth", "sunAircraftNoseAzimuth" ]
converters = {'barometricAltitude': ft2m, 'radarAltitude': ft2m, 'aircraftIndicatedAirSpeed': kts2mps, 'aircraftTime': iso2aircraftTime}
In [53]:
df = pd.read_csv(fn_iwg1, names=columns, header=None, converters=converters) # parse_dates=True,
In [54]:
df[1000:1005]
Out[54]:
In [7]:
W = Workspace('rdr20170411T161346end20170411T220220sdr20170411T224336')
In [13]:
acT, sT = W.aircraftTime[:], W.systemTimer[:]
r = np.argwhere((acT > 0) & (sT > 0)).squeeze()
len(r), r[:3]
Out[13]:
In [28]:
# confirm that records of interest are monotonic
r[np.argwhere(np.diff(acT[r])<0).squeeze()], r[np.argwhere(np.diff(sT[r])<0).squeeze()]
Out[28]:
In [29]:
sT[31985:31999]
Out[29]:
In [32]:
# choose subset for polyfit, given descent data has one flaky point in systemTimer
rsub = r[:8000]
np.max(rsub)
Out[32]:
In [56]:
m,b = np.polyfit(acT[rsub], sT[rsub], 1)
m,b
Out[56]:
In [59]:
fig,ax = fax()
rss = rsub[:400]
ax.plot(acT[rss], sT[rss], '.')
ax.plot(acT[rss], m*acT[rss] + b)
Out[59]:
In [65]:
# estimate systemTimer based on IWG1 aircraftTime, noting that IWG1 starts before SHIS
esT = df['aircraftTime'] * m + b
esT[6000:6020]
Out[65]:
In [68]:
# piecewise linear interpolate IWG1 to actual systemTimer values
D = {}
rzero = np.argwhere(sT==0).squeeze()
for c in columns[1:]:
if '_' in c:
continue
print(c)
q = np.require(np.interp(sT, esT, np.array(df[c])), dtype=np.float32)
q[rzero] = np.NAN
D[c] = q
In [75]:
fig,ax = fax()
ax.plot(D['instrumentLongitude'], D['instrumentLatitude'], '.')
ax.plot(W['instrumentLongitude'][:][r], W['instrumentLatitude'][:][r], '.')
Out[75]:
In [79]:
fig,ax = fax()
rss = r[2000:2100]
ax.plot(sT[rss], D['instrumentLatitude'][rss], '.')
ax.plot(sT[rss], W['instrumentLatitude'][:][rss], '.')
ax.legend(['interp latitude', 'orig latitude'])
Out[79]:
In [80]:
fig,ax = fax()
ax.plot(D['Roll'], '.')
ax.plot(D['Pitch'], '.')
ax.legend(['Roll', 'Pitch'])
Out[80]:
In [83]:
print("mkdir orig-nav")
for c in D.keys():
print("cp -n {}.real4 orig-nav/".format(c))
In [84]:
import os
nn = os.path.join(W.path,'new-nav')
if not os.path.isdir(nn):
os.mkdir(nn)
for k,v in D.items():
path = os.path.join(nn, k+'.real4')
print('writing {}'.format(path))
mm = np.memmap(path, dtype=np.float32, mode='w+', shape=v.shape)
mm[:] = v
mm.flush()
In [ ]: