In [3]:
from datetime import datetime, timedelta

import matplotlib.pyplot as plt
import numpy as np

from geopak.astro import transform
from geopak.ellps import ecef2geod, ecef2geod_surf
from geopak.gransim import calc_surface_observation, get_nearest_tle
from geopak.leap import utc2tai
from geopak.sgp4 import sgp4
from geopak.tle import decode_tle
In [9]:
# simulated TLE; we'll make a quick plot showing location and altitude
# over a couple hours starting from the TLE epoch
tle = (
    "1 00001U          22043.00000000 -.00000735  00000-0 -79424-4 0 00002\n" +
    "2 00001 098.6022 182.1527 0010584 271.5423 088.4558 14.88397198000015\n"
)
satrec = sgp4.init_from_tle(tle)
t = utc2tai(datetime(2022, 2, 12)) + 60e6 * np.arange(120)
pos, vel = sgp4.run(satrec, t)
pos, vel = transform("teme", "ecef", t, pos, vel)
lat, lon, alt = ecef2geod(pos)
plt.scatter(np.degrees(lon), np.degrees(lat), c=alt)
plt.xlabel("Longitude")
plt.xlim(-180, 180)
plt.ylabel("Latitude")
plt.ylim(-90, 90)
plt.colorbar(label="Altitude (km)")
Out[9]:
<matplotlib.colorbar.Colorbar at 0x7f14961e8f40>
In [10]:
# here we define a bunch of parameters, mostly grabbed from the
# "MWinds Obervatory" scan geometry diagram I was given

# we'll make a sample 10-scan granule starting from the TLE epoch
granule_start = utc2tai(datetime(2022, 2, 12))
num_atrack = 10

# dimensions of each individual scan
num_xtrack = 1068
num_fov = 512

# it seems very likely these values need to be adjusted, espcially
# the along-track "pitch", given what the results below look like...
max_scan_angle = np.radians(52)
pitch_angle_step = 1.7e-3  # 1.7 milliradians
max_pitch_angle = pitch_angle_step * (num_xtrack - 1) / 2

scan_duration = 93e6  # earth view time for each scan
scan_period = 110e6  # interval from one scan to the next
In [12]:
# now produce the simulated geolocation
obs_time = (
    granule_start
    + scan_period * np.arange(num_atrack)[:, np.newaxis]
    + np.linspace(0, scan_duration, num_xtrack)[np.newaxis, :]
)

# scan is left-to-right, and we order FOVs at each scan position
# aft-to-fore
scan_angle = np.linspace(max_scan_angle, -max_scan_angle, num_xtrack)
pitch_angle = np.linspace(-max_pitch_angle, max_pitch_angle, num_fov)

obs = calc_surface_observation(
    satrec,
    obs_time[:, :, np.newaxis],
    scan_angle[np.newaxis, :, np.newaxis],
    pitch_angle[np.newaxis, np.newaxis, :],
)
lat, lon = ecef2geod_surf(obs)
In [19]:
# plot a subset of the obs, colored by obervation
plt.figure(figsize=[16,8])
idx = (slice(None), slice(None, None, 20), slice(None, None, 16))
x, y, c = [
    a.ravel() for a in np.broadcast_arrays(
        np.degrees(lon[idx]),
        np.degrees(lat[idx]),
        obs_time[idx[:2] + (np.newaxis,)]
    )
]
plt.scatter(x, y, c=c)
plt.xlabel("Longitude")
plt.ylabel("Latitude")
Out[19]:
Text(0, 0.5, 'Latitude')
In [20]:
# based on the observation spacing near the scan edges, it looks like
# we may be viewing off the limb. if that's happening we'll get NaN
# values for latitude and longitude...
nan_flag = np.isnan(lat)
print(f"{nan_flag.sum()} of {nan_flag.size} observations could not be geolocated")

14078 of 5468160 observations could not be geolocated