from glob import glob
import cartopy.crs as ccrs
import h5py
import matplotlib.pyplot as plt
import netCDF4
import numpy as np
import pandas as pd
import pyproj
import seaborn as sns
from astropy.time import Time, TimeDelta
sns.set(style='darkgrid', context='talk')
def read_l1b(l1b_file):
with netCDF4.Dataset(l1b_file) as d:
frame = pd.DataFrame({
'obs_time': np.repeat(d['obs_time_tai93'][:][:,:,None], 9, axis=2).ravel(),
'lat': d['lat'][:].ravel(),
'lon': d['lon'][:].ravel(),
'lat_geoid': d['lat_geoid'][:].ravel(),
'lon_geoid': d['lon_geoid'][:].ravel(),
'sat_zen': d['sat_zen'][:].ravel(),
'sat_azi': d['sat_azi'][:].ravel(),
'sol_zen': d['sol_zen'][:].ravel(),
'sol_azi': d['sol_azi'][:].ravel(),
'surf_alt': d['surf_alt'][:].ravel()})
frame['obs_time'] = (Time('1993-01-01') + TimeDelta(frame.obs_time, format='sec')).datetime
frame['fov'] = np.arange(len(frame)) % 9 + 1
frame = frame.set_index(['obs_time', 'fov'])
return frame
def read_sdr(gcrso_file):
with h5py.File(gcrso_file, 'r') as f:
g = f['All_Data/CrIS-SDR-GEO_All']
frame = pd.DataFrame({
'obs_time': np.repeat(g['FORTime'][:][:,:,None], 9, axis=2).ravel(),
'lat': g['Latitude'][:].ravel(),
'lon': g['Longitude'][:].ravel(),
'sat_zen': g['SatelliteZenithAngle'][:].ravel() % 360,
'sat_azi': g['SatelliteAzimuthAngle'][:].ravel() % 360,
'sol_zen': g['SolarZenithAngle'][:].ravel(),
'sol_azi': g['SolarAzimuthAngle'][:].ravel()})
frame['obs_time'] = (
Time('1958-01-01', scale='tai')
+ TimeDelta(frame.obs_time * 1e-6, format='sec')
).utc.datetime
frame['fov'] = np.arange(len(frame)) % 9 + 1
frame = frame.set_index(['obs_time', 'fov'])
return frame
def show_diff(data):
sns.distplot(data, kde=False)
return data.describe()
def diff_azi(x, y):
"""Calculate difference between azimuth values; assumes inputs are in [0,360)"""
diff = np.abs(x - y)
return np.minimum(diff, 360 - diff)
l1b = pd.concat([read_l1b(f) for f in sorted(glob('data/SNDR.*'))])
sdr = pd.concat([read_sdr(f) for f in sorted(glob('data/GCRSO_*'))])
l1b, sdr = l1b.align(sdr, 'inner', axis=0)
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
ax.scatter(l1b.lon, l1b.lat, transform=ccrs.Geodetic())
print('Total FOVs:', len(l1b))
Total FOVs: 214110
l1b.iloc[[0, -1]]
| lat | lon | lat_geoid | lon_geoid | sat_zen | sat_azi | sol_zen | sol_azi | surf_alt | ||
|---|---|---|---|---|---|---|---|---|---|---|
| obs_time | fov | |||||||||
| 2021-08-24 13:14:47.983970 | 1 | -11.103292 | -6.045043 | -11.103292 | -6.045043 | 59.702019 | 83.279739 | 25.074383 | 330.984680 | 0.0 |
| 2021-08-24 14:59:57.783950 | 9 | 4.887642 | -15.219758 | 4.887642 | -15.219758 | 57.631771 | 259.607727 | 29.524368 | 283.504181 | 0.0 |
geod = pyproj.Geod(ellps='WGS84')
dist = pd.Series(
geod.inv(l1b.lon_geoid.values, l1b.lat_geoid.values, sdr.lon.values, sdr.lat.values)[2],
l1b.index,
name='dist_vs_sdr')
show_diff(dist)
count 164310.000000 mean 15.430414 std 17.888645 min 0.000000 25% 3.614095 50% 8.214800 75% 20.102128 max 129.629733 Name: dist_vs_sdr, dtype: float64
no_tc = l1b.surf_alt < 100 # exclude cases where terrain correction is significant
l1b = l1b[no_tc]
sdr = sdr[no_tc]
len(l1b)
164310
show_diff(l1b.sat_zen - sdr.sat_zen)
count 164310.000000 mean -0.000103 std 0.000270 min -0.006184 25% -0.000240 50% -0.000078 75% 0.000065 max 0.001240 Name: sat_zen, dtype: float64
show_diff(diff_azi(l1b.sat_azi, sdr.sat_azi))
count 164310.000000 mean 0.000699 std 0.001785 min 0.000000 25% 0.000137 50% 0.000305 75% 0.000641 max 0.352631 Name: sat_azi, dtype: float64
show_diff(l1b.sol_zen - sdr.sol_zen)
count 164310.000000 mean 0.000024 std 0.000151 min -0.004028 25% -0.000015 50% 0.000023 75% 0.000061 max 0.002097 Name: sol_zen, dtype: float64
show_diff(diff_azi(l1b.sol_azi, sdr.sol_azi))
count 164310.000000 mean 0.000189 std 0.001057 min -0.010834 25% 0.000027 50% 0.000103 75% 0.000244 max 0.352783 Name: sol_azi, dtype: float64
