import os
import pandas as pd
import xarray as xr
import numpy as np
import sasktran as sk
from skretrieval.util import rotation_matrix
from skretrieval.core import OpticalGeometry
[docs]
def load_er2_data(flight, folder=r'C:\Users\lar555\Documents\ACCP\ER2\telemetry', resample='400ms', load=True):
"""
Load the IWG1 data for a flight into a xarray.Dataset.
Parameters
----------
flight: Name of the flight used in the IWG1 file
folder: directory where the file is stored
resample: If set data will be resampled to this resolution
load: If true the resampled file will be loaded if it exists
Returns
-------
xarray.Dataset containing the IWG1 data
"""
if load:
try:
data_lores = xr.open_dataset(os.path.join(folder, f'IWG1.{flight}_{resample}.nc'))
return data_lores
except FileNotFoundError:
load = False
data = pd.read_csv(os.path.join(folder, f'IWG1.{flight}.txt'),
index_col=0,
usecols=[1, 2, 3, 4, 8, 13, 14, 16, 17],
names=['time', 'latitude', 'longitude', 'altitude', 'speed', 'heading', 'track',
'pitch_angle', 'roll_angle'],
parse_dates=True, skiprows=11863).to_xarray().sortby('time')
data = data.where(np.isfinite(data.latitude) & np.isfinite(data.roll_angle), drop=True)
data_lores = data.resample(time=resample).mean()
data_lores.to_netcdf(os.path.join(folder, f'IWG1.{flight}_{resample}.nc'))
return data_lores
[docs]
def er2_orientation(er2, instrument_pitch_deg=0.0, instrument_yaw_deg=0.0, instrument_roll_deg=0.0):
"""
Calculate the `OpticalGeometry` used in Sasktran from an IWG1 point
Parameters
----------
er2: xr.Dataset
Containing a single time and at least attributes: ['longitude', 'latitude', 'altitude', 'heading', 'pitch', 'roll', 'time']
instrument_pitch_deg: [Optional]
Pitch angle in degrees of the instrument relative to teh aircraft
instrument_yaw_deg: [Optional]
Yaw angle in degrees of the instrument relative to teh aircraft
instrument_roll_deg: [Optional]
Roll angle in degrees of the instrument relative to teh aircraft
Returns
-------
LookGeometry
"""
try:
geo = sk.Geodetic()
geo.from_lat_lon_alt(float(er2.latitude.values), float(er2.longitude.values), float(er2.altitude.values))
pos = geo.location
north = -geo.local_south
# create frame along heading direction
up = geo.local_up
look = north @ rotation_matrix(up, (float(er2.heading.values) + instrument_yaw_deg) * np.pi / 180)
hor = np.cross(look, up)
# rotate frame by pitch angle
look = look @ rotation_matrix(hor, -(float(er2.pitch_angle.values) + instrument_pitch_deg) * np.pi / 180)
up = up @ rotation_matrix(hor, -(float(er2.pitch_angle.values) + instrument_pitch_deg) * np.pi / 180)
# rotate frame by roll angle
up = up @ rotation_matrix(look, -(float(er2.roll_angle.values) + instrument_roll_deg) * np.pi / 180)
return OpticalGeometry(observer=pos,
look_vector=look,
local_up=up,
mjd=(er2.time.values - np.datetime64('1858-11-17')) / np.timedelta64(1, 'D'))
except Exception as e:
print(e)
return None