In [1]:
%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 [2]:
import numpy as np
import pandas as pd
from datetime import datetime
from calendar import timegm
from pyfbf import Workspace
In [3]:
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 [4]:
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 [5]:
df = pd.read_csv(fn_iwg1, names=columns, header=None, converters=converters) # parse_dates=True,
In [6]:
df[1000:1005]
Out[6]:
In [7]:
W = Workspace('rdr20170411T161346end20170411T220220sdr20170414T160410')
In [8]:
acT, sT = W.aircraftTime[:], W.systemTimer[:]
r = np.argwhere((acT > 0) & (sT > 0)).squeeze()
len(r), r[:3]
Out[8]:
In [9]:
# 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[9]:
In [10]:
sT[31985:31999]
Out[10]:
In [11]:
# choose subset for polyfit, given descent data has one flaky point in systemTimer
rsub = r[:8000]
np.max(rsub)
Out[11]:
In [12]:
m,b = np.polyfit(acT[rsub], sT[rsub], 1)
m,b
Out[12]:
In [13]:
fig,ax = fax()
rss = rsub[:400]
ax.plot(acT[rss], sT[rss], '.')
ax.plot(acT[rss], m*acT[rss] + b)
Out[13]:
In [14]:
# estimate systemTimer based on IWG1 aircraftTime, noting that IWG1 starts before SHIS
esT = df['aircraftTime'] * m + b
esT[6000:6020]
Out[14]:
In [15]:
# 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 [16]:
fig,ax = fax()
ax.plot(D['instrumentLongitude'], D['instrumentLatitude'], '.')
ax.plot(W['instrumentLongitude'][:][r], W['instrumentLatitude'][:][r], '.')
Out[16]:
In [17]:
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[17]:
In [18]:
fig,ax = fax()
ax.plot(D['Roll'], '.')
ax.plot(D['Pitch'], '.')
ax.legend(['Roll', 'Pitch'])
Out[18]:
In [21]:
print("mkdir orig-nav")
for c in D.keys():
print("cp -n {}.real4 orig-nav/".format(c))
In [23]:
end=21753
acT[end-1], acT[end]
Out[23]:
In [24]:
# patch nav in-place
import os
for k,v in D.items():
path = os.path.join(W.path, k+'.real4')
print('patching {}'.format(path))
mm = np.memmap(path, dtype=np.float32, mode='r+', shape=v.shape)
mm[:end] = v[:end]
mm.flush()
In [25]:
# correct aircraftType
fn_type = os.path.join(W.path, 'aircraftType.real4')
aircraftType = np.memmap(fn_type, mode='r+', dtype=np.float32)
aircraftType[:] = 871
aircraftType.flush()
In [ ]: