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.
julia
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
julia
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.95316953931108, -0.21383296251610504), (72.44341434106293, -0.10672778602982819), (72.93373095997808, 0.0004017897759837355), (73.42404759167403, 0.1075313359553473), (73.91429243175713, 0.21463642358191728)]
[(71.79837594238133, 0.5028374122027154), (72.28863153942606, 0.6100252546414895), (72.77898176649367, 0.7171846620027664), (73.26935481103858, 0.8242911827619259), (73.75967884046462, 0.931320400452978)]
[(71.6433977098909, 1.219513361706255), (72.13373868013473, 1.3267894528939954), (72.62419714561938, 1.4339842635475444), (73.11470125603013, 1.5410733112212138), (73.6051791209236, 1.6480321835895793)]
[(71.48820192958809, 1.9361896469123565), (71.97870290955211, 2.0435596067709447), (72.46934432351597, 2.150795415124998), (72.96005425373382, 2.257872550235359), (73.45076072217834, 2.364766595584615)]
[(71.3327561694158, 2.652861002865055), (71.82349190840712, 2.7603304980312546), (72.31439111503418, 2.8676129307324056), (72.80538177435713, 2.9746837315202774), (73.29639179089084, 3.0815184713269894)]
[(71.17702708857452, 3.3695221375317663), (71.66807250302799, 3.4770968910121742), (72.15930453517106, 3.584431616561461), (72.65065104296089, 3.6915016885501646), (73.14203978339744, 3.7982826569622774)]
[(71.02098157928754, 4.0861677305483335), (71.5124118073706, 4.1938535314550185), (72.0040519414144, 4.301246269888317), (72.49582968217197, 4.4083212554952915), (72.98767260884216, 4.515054008861582)]
[(70.86458621368743, 4.802792431896424), (71.35647667026146, 4.910595145279111), (71.84860048139575, 5.018051677927566), (72.3408851603, 5.12513726613932), (72.83325807780953, 5.231827392585305)]
[(70.7078067058724, 5.519390860500472), (71.20023313790419, 5.627316437291022), (71.69291655590268, 5.734842616637162), (72.18578425490976, 5.841944552763036), (72.67876336651078, 5.9485976818286)]
[(70.5506088682089, 6.235957602743508), (71.04364741072816, 6.344012089829161), (71.53696677632965, 6.4516138494749775), (72.03049401093037, 6.5587379449957925), (72.5241559756064, 6.665359757350346)]
⋮
[(-76.82197859245451, -5.897124971978316), (-76.32916534201979, -5.790479822068719), (-75.83646182686881, -5.683389917626609), (-75.34394084406787, -5.575879804703733), (-74.85167502888359, -5.467974308458096)]
[(-76.97644969742406, -5.180675874815813), (-76.48422144108147, -5.073992139037255), (-75.99207977361321, -4.96691691082197), (-75.50009726212129, -4.85947469215276), (-75.00834633292871, -4.751690228646851)]
[(-77.13084943483918, -4.464180422994778), (-76.63913013013008, -4.357452476184785), (-76.1474743841167, -4.250386194336345), (-75.65595456470444, -4.1430060439353245), (-75.16464291978717, -4.035336699745295)]
[(-77.28521008111518, -3.7476436900383363), (-76.79392406139483, -3.6408659454247165), (-76.30267866463637, -3.533802933329968), (-75.81154609000964, -3.4264790932396565), (-75.3205984372257, -3.3189190376576723)]
[(-77.4395642993574, -3.0310707667802923), (-76.94863621769504, -2.924237666374511), (-76.45772589617594, -2.8171722909520263), (-75.96690539543613, -2.7098990614330734), (-75.47624669702888, -2.6024425365765467)]
[(-77.59394399667727, -2.3144667624356456), (-77.10329877032653, -2.2075727674808814), (-76.61264849288906, -2.100499429537505), (-76.12206511645637, -1.9932711593472217), (-75.63162053427352, -1.8859124703729877)]
[(-77.74838171241814, -1.5978368056566623), (-77.25794446776102, -1.4908763871151185), (-76.7674793912, -1.3837895117585668), (-76.27705835610469, -1.276600588502126), (-75.78675319713173, -1.1693340939097034)]
[(-77.90290948043597, -0.8811860455800905), (-77.41260549830511, -0.7741536746443439), (-76.92225091278901, -0.667047701736426), (-76.43191754823722, -0.5598925422751898), (-75.94167721035298, -0.4527126442835095)]
[(-78.05755937064266, -0.16451965288243084), (-77.56731403192795, -0.057409791495912604), (-77.0769953066533, 0.0497208338686246), (-76.58667499977423, 0.15684779296675874), (-76.09642491763586, 0.26394665798430633)]Plot the sweep points on a globeaxis, just to get an idea.
julia
scatter(reduce(vcat, sweep_points_timeseries); axis = (; type = GlobeAxis))
Sample a raster, at the points at which the satellite's sensor would "sample".
julia
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)┌ 2160×1080 Raster{Float64, 2} ┐
├──────────────────────────────┴───────────────────────────────────────── dims ┐
↓ X Projected{Float64} -180.0:0.16666666666666666:179.83333333333331 ForwardOrdered Regular Intervals{Start},
→ Y Projected{Float64} 89.83333333333333:-0.16666666666666666:-90.0 ReverseOrdered Regular Intervals{Start}
├──────────────────────────────────────────────────────────────────── metadata ┤
Metadata{Rasters.GDALsource} of Dict{String, Any} with 1 entry:
"filepath" => "/home/runner/.julia/artifacts/WorldClim/Climate/tmin/wc2.1_10m…
├────────────────────────────────────────────────────────────────────── raster ┤
missingval: NaN
extent: Extent(X = (-180.0, 179.99999999999997), Y = (-90.0, 90.0))
crs: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.25722...
└──────────────────────────────────────────────────────────────────────────────┘
↓ → 89.8333 89.6667 89.5 … -89.6667 -89.8333 -90.0
-180.0 NaN NaN NaN -22.1597 -19.9805 -20.338
-179.833 NaN NaN NaN -22.905 -21.058 -21.01
⋮ ⋱ ⋮
179.667 NaN NaN NaN -24.4398 -22.4585 -22.41
179.833 NaN NaN NaN … -22.9375 -20.7255 -21.0588Re-sample the raster to the sampling timeseries. This is currently a manual process, but we should just make everything a raster here.
julia
xys = reduce(hcat, sweep_points_timeseries)
vals = map(xys) do xy
ras[X(Contains(xy[1])), Y(Contains(xy[2]))]
end5×3001 Matrix{Float64}:
NaN NaN NaN NaN NaN NaN NaN … 19.1945 18.6513 14.9392 2.512
NaN NaN NaN NaN NaN NaN NaN 20.4805 19.8312 19.1225 13.2595
NaN NaN NaN NaN NaN NaN NaN 20.7752 20.477 20.1785 20.035
NaN NaN NaN 26.55 NaN NaN NaN 21.0618 20.752 20.6007 20.2668
NaN NaN NaN 26.5 NaN NaN NaN 21.2065 21.002 20.918 20.6645Finally, we plot the results on a GeoMakie GlobeAxis.
julia
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.