Satellite sweep plot

Here, we'll plot the data collected by a satellite as it orbits the Earth. Note that we've chosen to fake the data by simulating the satellite sweep over a high-level data product. See the satellite_sweep_plot.jl example in GeoMakie/examples/specialized/satellite/ for a more complete example using real data.

using LinearAlgebra: normalize
using Geodesy
using GeometryBasics: Point2d
using SatelliteToolbox, SatelliteAnalysis
using GLMakie, GeoMakie

"""
    sweep_points(position, velocity, distances)

Sweep points along a line of sight, given a position and velocity in the ECEF frame,
and distances in meters.

# Arguments
- `position`: The position of the satellite.
- `velocity`: The velocity of the satellite.
- `distances`: The distances to sweep along the line of sight, in meters.  Will be interpreted as meters in the ENU frame.

# Returns
A vector of 2-tuples of (longitude, latitude) in degrees, representing the sweep points according to `distances`.

# Extended help

# How does this work?

The algorithm is pretty simple.  Given a position and velocity in the earth-centered earth-fixed frame (ECEF),
we can calculate the position and velocity in the local frame of the satellite's position on the surface
of the ellipsoidal Earth.  This is the East-North-Up (ENU) frame.

We can then calculate the points along the line of sight by sweeping along a vector perpendicular to the velocity vector,
in this ENU space which is essentially the projection of the local tangent plane to the sphere.
This vector can be adjusted, rotated, etc. as necessary to change the angle of the sensor etc.

From ENU space offsets we go back to long-lat space, which is the standard way to represent positions on the Earth.

# Performance

This runs with 2 allocations (not sure why 2, probably one of the transforms)
and in ~6 μs for 50 points on my machine.  Scales linearly with `length(distances)`.
"""
function sweep_points(position, velocity, distances)
    # Input is in ECEF space.
    # Operate in ENU space and then convert to long lat
    enu_transf = ENUfromECEF(position, wgs84)
    lla_transf = LLAfromENU(position, wgs84)
    enu_position = ENU(0.0, 0.0, 0.0) # by definition of the enu frame
    enu_velocity = enu_transf(position +velocity) # should be enu(position + velocity) - enu(position) but we don't actually need to do that

    # Calculate the direction of the velocity
    enu_direction = normalize(Point2d(enu_velocity.e, enu_velocity.n))
    sweep_left_direction = Point2d(enu_direction[2], -enu_direction[1])

    # Calculate the points along the line of sight
    sweep_points = map(distances) do distance
        p = Point2d(enu_position.e, enu_position.n) + distance * sweep_left_direction
        lla = lla_transf(ENU(p[1], p[2], 0.0))
        (lla.lon, lla.lat)
    end
end
# Convenience method for propagating a satellite and then sweeping points along its line of sight
sweep_points(sv::SatelliteToolbox.OrbitStateVector, distances) = sweep_points(ECEF(sv.r), ECEF(sv.v), distances)
# Simple example usage
JFK = LLA(; lat = 40.6413, lon = -73.7781, alt = 0.0)
sp = sweep_points(ECEFfromLLA(wgs84)(JFK), ECEF(0, 0, 1), LinRange(0, 49, 50))

50-element Vector{Tuple{Float64, Float64}}:
 (-73.7781, 40.64130000000001)
 (-73.77808817822306, 40.6412999999994)
 (-73.77807635644612, 40.64129999999758)
 (-73.77806453466918, 40.64129999999455)
 (-73.77805271289223, 40.641299999990316)
 (-73.7780408911153, 40.64129999998488)
 (-73.77802906933834, 40.64129999997822)
 (-73.77801724756141, 40.64129999997035)
 (-73.77800542578446, 40.64129999996127)
 (-73.77799360400752, 40.64129999995098)
 ⋮
 (-73.77761530714538, 40.64129999898284)
 (-73.77760348536845, 40.64129999893263)
 (-73.7775916635915, 40.64129999888119)
 (-73.77757984181457, 40.64129999882855)
 (-73.77756802003762, 40.6412999987747)
 (-73.77755619826067, 40.64129999871964)
 (-73.77754437648373, 40.641299998663364)
 (-73.7775325547068, 40.64129999860587)
 (-73.77752073292984, 40.64129999854718)

A more complex example, simulating a satellite

amz1_tle = SatelliteToolbox.tle"""
       AMAZONIA 1
       1 47699U 21015A   25205.73369244  .00000207  00000+0  78058-4 0  9996
       2 47699  98.3576 279.7581 0001748  96.4737 263.6651 14.40869396231498
"""

prop = SatelliteToolbox.Propagators.init(Val(:SGP4), amz1_tle)

sv_teme = SatelliteToolbox.Propagators.propagate!(prop, 0:12:(3600*10), OrbitStateVector)
eop = SatelliteToolbox.fetch_iers_eop()
sv_itrf = SatelliteToolbox.sv_eci_to_ecef.(sv_teme, (TEME(),), (SatelliteToolbox.ITRF(),), (eop,))

sweep_points_timeseries = sweep_points.(sv_itrf, (LinRange(-125000, 125000, 5),))

3001-element Vector{Vector{Tuple{Float64, Float64}}}:
 [(71.95316953930954, -0.21383296238863847), (72.44341434106218, -0.10672778590591835), (72.93373095997808, 0.00040178989632720587), (73.4240475916748, 0.10753133607211557), (73.9142924317587, 0.2146364236951021)]
 [(71.79837594238462, 0.5028374123312584), (72.28863153943013, 0.6100252547664783), (72.77898176649852, 0.717184662124191), (73.2693548110442, 0.824291182879777), (73.75967884047105, 0.9313204005672475)]
 [(71.643397709899, 1.2195133618358747), (72.13373868014362, 1.326789453020063), (72.62419714562907, 1.433984263670049), (73.11470125604062, 1.5410733113401467), (73.60517912093486, 1.6480321837049312)]
 [(71.48820192960099, 1.9361896470430566), (71.97870290956583, 2.043559606898093), (72.4693443235305, 2.150795415248583), (72.96005425374915, 2.2578725503553723), (73.45076072219447, 2.364766595701047)]
 [(71.3327561694335, 2.652861002996834), (71.82349190842565, 2.7603304981594827), (72.31439111505352, 2.8676129308570704), (72.80538177437731, 2.9746837316413686), (73.29639179091184, 3.0815184714444985)]
 [(71.17702708859703, 3.3695221376646294), (71.66807250305132, 3.477096891141483), (72.15930453519523, 3.5844316166872057), (72.6506510429859, 3.691501688672333), (73.1420397834233, 3.7982826570808625)]
 [(71.02098157931484, 4.086167730682282), (71.51241180739876, 4.193853531585412), (72.00405194144341, 4.301246270015143), (72.49582968220183, 4.408321255618538), (72.98767260887288, 4.515054008981238)]
 [(70.86458621371953, 4.802792432031457), (71.35647667029443, 4.910595145410586), (71.8486004814296, 5.018051678055471), (72.34088516033474, 5.125137266263643), (72.83325807784513, 5.231827392706035)]
 [(70.7078067059093, 5.5193908606365945), (71.20023313794199, 5.627316437423581), (71.69291655594137, 5.734842616766145), (72.18578425494935, 5.841944552888433), (72.67876336655124, 5.9485976819504)]
 [(70.55060886825062, 6.235957602880702), (71.04364741077077, 6.344012089962797), (71.53696677637318, 6.4516138496050415), (72.03049401097479, 6.558737945122274), (72.52415597565175, 6.665359757473237)]
 ⋮
 [(-76.82197859240824, -5.897124971829325), (-76.32916534197449, -5.790479821915833), (-75.83646182682448, -5.683389917469842), (-75.34394084402452, -5.575879804543099), (-74.8516750288412, -5.467974308293604)]
 [(-76.97644969738342, -5.180675874667987), (-76.4842214410418, -5.073992138885552), (-75.99207977357447, -4.966916910666404), (-75.5000972620835, -4.859474691993342), (-75.00834633289188, -4.751690228483595)]
 [(-77.13084943480415, -4.464180422848143), (-76.639130130096, -4.357452476034279), (-76.14747438408355, -4.250386194181981), (-75.65595456467223, -4.1430060437771115), (-75.16464291975588, -4.035336699583248)]
 [(-77.28521008108578, -3.747643689892926), (-76.79392406136635, -3.6408659452754235), (-76.30267866460879, -3.533802933176803), (-75.81154608998298, -3.4264790930826345), (-75.32059843719996, -3.318919037496805)]
 [(-77.43956429933361, -3.0310707666360495), (-76.94863621767215, -2.9242376662264045), (-76.45772589615395, -2.817172290800063), (-75.96690539541504, -2.709899061277267), (-75.47624669700869, -2.602442536416911)]
 [(-77.5939439966591, -2.3144667622926036), (-77.10329877030925, -2.207572767333974), (-76.61264849287268, -2.100499429386744), (-76.12206511644085, -1.9932711591926184), (-75.63162053425889, -1.8859124702145547)]
 [(-77.74838171240559, -1.597836805514851), (-77.25794446774934, -1.490876386969428), (-76.76747939118918, -1.3837895116090078), (-76.27705835609476, -1.27660058834871), (-75.78675319712265, -1.169334093752442)]
 [(-77.90290948042906, -0.881186045439456), (-77.41260549829906, -0.7741536744998433), (-76.92225091278378, -0.6670477015880695), (-76.43191754823286, -0.5598925421229886), (-75.94167721034945, -0.4527126441274755)]
 [(-78.05755937064139, -0.16451965274303415), (-77.5673140319275, -0.057409791352632904), (-77.0769953066537, 0.04972083401577737), (-76.58667499977545, 0.15684779311777372), (-76.0964249176379, 0.26394665813917184)]

Plot the sweep points on a globeaxis, just to get an idea.

scatter(reduce(vcat, sweep_points_timeseries); axis = (; type = GlobeAxis))

Sample a raster, at the points at which the satellite's sensor would "sample".

using Rasters, RasterDataSources
import ArchGDAL, NCDatasets # to activate extensions on rasters for file IO
worldclim_file = RasterDataSources.getraster(WorldClim{Climate}, :tmin; month = 1)
ras = Raster(worldclim_file)
ras = replace_missing(ras, NaN)

Re-sample the raster to the sampling timeseries. This is currently a manual process, but we should just make everything a raster here.

xys = reduce(hcat, sweep_points_timeseries)
vals = map(xys) do xy
    ras[X(Contains(xy[1])), Y(Contains(xy[2]))]
end

Finally, we plot the results on a GeoMakie GlobeAxis.

f, a, p = surface(
    first.(xys), last.(xys), zeros(size(xys));
    color = vals,
    shading = NoShading,
    axis = (; type = GlobeAxis, show_axis = false)
)
# Create a "background" mesh such that the satellite traces
# on the other side of the globe (from the camera's point of view) are faded out.
bg = meshimage!(a, -180..180, -90..90, reshape([RGBAf(1,1,1,0.5)], 1, 1); uv_transform = :rotr90, zlevel = -100_000, reset_limits = false)
lines!(a, GeoMakie.coastlines(); color = (:black, 0.5), linewidth = 1)
f

---

This page was generated using Literate.jl.