diff --git a/Cargo.lock b/Cargo.lock index bcea77f..0274d6d 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -509,6 +509,7 @@ checksum = "fac4744fb15ae8337dc853fee7fb3f4e48c0fbaa23d0afe49c447b4fab126118" dependencies = [ "iana-time-zone", "num-traits", + "serde", "windows-link", ] @@ -2137,6 +2138,17 @@ dependencies = [ "winapi", ] +[[package]] +name = "sgp4" +version = "2.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9467b9a7be8485ed8be0f336d399c8f32c0fcd60686e7dd2ed3dab75c9a73eb3" +dependencies = [ + "chrono", + "serde", + "serde_json", +] + [[package]] name = "sha1" version = "0.10.6" @@ -2708,6 +2720,7 @@ dependencies = [ "flate2", "serde", "serde_json", + "sgp4", "tokio", "tracing", "uuid", diff --git a/docs/Wxsat-Map-Overlay.md b/docs/Wxsat-Map-Overlay.md new file mode 100644 index 0000000..c807698 --- /dev/null +++ b/docs/Wxsat-Map-Overlay.md @@ -0,0 +1,152 @@ +# Weather Satellite Map Overlay Integration + +Overlay decoded NOAA APT and Meteor-M LRPT satellite images on the Leaflet +map module, with ground track visualisation and source filtering. + +*Created: 2026-03-28* + +## Status + +| Step | Description | Status | +|------|-------------|--------| +| 1 | Add `sgp4` crate, create `trx-core/src/geo.rs` | Done | +| 2 | Extend `WxsatImage`/`LrptImage` with geo fields | Done | +| 3 | Compute geo-bounds in `finalize_wxsat_pass` / `finalize_lrpt_pass` | Done | +| 4 | Add `wxsat` to map source filter + image overlay rendering | Done | +| 5 | Add ground track polyline + filter toggle UI | Done | +| 6 | Build, test, verify | Done | + +## Motivation + +The wxsat plugin currently shows a history table with download links but has +no geographic context. Since the Map module already renders APRS, AIS, VDES, +and FTx/WSPR positions, weather satellite images are a natural addition — they +can be projected as semi-transparent overlays on the same Leaflet map. + +## Architecture + +### Data flow + +``` +Pass decoded (APT / LRPT) + ↓ finalize_wxsat_pass / finalize_lrpt_pass (trx-server/audio.rs) + ↓ SGP4 propagation using satellite TLE + pass timestamps + ↓ Compute geo_bounds [[south, west], [north, east]] + ↓ Compute ground_track [[lat, lon], ...] + ↓ Attach to WxsatImage / LrptImage + ↓ Broadcast via DecodedMessage + ↓ SSE → browser + ↓ wxsat.js: L.imageOverlay() + L.polyline() on aprsMap +``` + +### Geo-referencing strategy + +Weather satellites (NOAA POES, Meteor-M) fly sun-synchronous polar orbits at +~850 km altitude with known TLE parameters. Given: + +- **Satellite identity** (from telemetry: NOAA-15/18/19, Meteor-M N2-3/N2-4) +- **Pass start/end timestamps** (`pass_start_ms`, `pass_end_ms`) +- **Receiver station lat/lon** (from `RigState.server_latitude/longitude`) + +We can use **SGP4 propagation** (via the `sgp4` crate) to compute the +sub-satellite ground track during the pass, then derive image bounds from the +known swath geometry: + +| Parameter | NOAA APT | Meteor LRPT | +|-----------|----------|-------------| +| Altitude | ~850 km | ~825 km | +| Swath width | ~2800 km | ~2800 km | +| Ground speed | ~6.9 km/s | ~6.9 km/s | +| Scan rate | 2 lines/sec (0.5s/line) | variable MCU rate | +| Image width | 909 px/channel | 1568 px | + +**Bounds computation:** +1. Propagate satellite position at `pass_start_ms` and `pass_end_ms` +2. Sub-satellite points define the ground track center line +3. Swath half-width (~1400 km) gives east/west extent +4. Image is projected as a simple lat/lon rectangle (acceptable distortion + for the typical ~15° latitude span of a single pass) + +**TLE source:** Hardcoded recent TLEs for the 5 active satellites, with an +optional HTTP refresh from CelesTrak. Stale TLEs (weeks old) still give +sub-degree accuracy for image overlay purposes. + +### Crate changes + +#### `trx-core` (src/trx-core/) + +New module `src/trx-core/src/geo.rs`: +- `SatelliteGeo` struct: holds hardcoded TLEs, provides `compute_pass_bounds()` +- `PassGeoBounds { south: f64, west: f64, north: f64, east: f64 }` +- `ground_track(sat, start_ms, end_ms) -> Vec<[f64; 2]>` +- Uses `sgp4` crate for orbital propagation +- Falls back to station-centered approximation when TLE unavailable + +`src/trx-core/src/decode.rs` — extend structs: +```rust +pub struct WxsatImage { + // ... existing fields ... + pub geo_bounds: Option<[f64; 4]>, // [south, west, north, east] + pub ground_track: Option>, // [[lat, lon], ...] +} +// Same for LrptImage +``` + +#### `trx-server` (src/trx-server/) + +`src/trx-server/src/audio.rs`: +- In `finalize_wxsat_pass`: after PNG write, call `SatelliteGeo::compute_pass_bounds()` + using satellite name, pass timestamps, and station lat/lon (threaded through + from config). Attach result to `WxsatImage`. +- Same for `finalize_lrpt_pass`. + +#### Frontend (trx-frontend-http/assets/web/) + +`plugins/wxsat.js`: +- On `onServerWxsatImage` / `onServerLrptImage`: if `geo_bounds` present, + call `window.addWxsatMapOverlay(msg)`. +- Manage overlay list, allow removal. + +`app.js`: +- Add `wxsat: false` to `DEFAULT_MAP_SOURCE_FILTER` (off by default to avoid + visual clutter; users opt-in). +- `window.addWxsatMapOverlay(msg)`: creates `L.imageOverlay(msg.path, bounds)` + with opacity 0.6, adds to `mapMarkers` set with `__trxType = "wxsat"`. +- `window.addWxsatGroundTrack(msg)`: creates `L.polyline(msg.ground_track)` + with dashed style. +- Overlay list in wxsat panel with per-image show/hide toggle. + +`index.html`: +- No structural changes needed; the map filter chip system auto-generates + from `DEFAULT_MAP_SOURCE_FILTER`. + +`style.css`: +- Styling for wxsat overlay opacity slider (future enhancement). + +## Dependencies + +| Crate | Version | Purpose | +|-------|---------|---------| +| `sgp4` | 2.4 | Pure Rust SGP4 orbital propagation | + +Added to `trx-core/Cargo.toml` (used by `geo.rs`). + +## Risk / Limitations + +- **Rectangular projection approximation**: The actual scan geometry is curved + (satellite moves along a great circle), but for a single pass spanning + ~15-20° of latitude, a lat/lon rectangle is a reasonable first approximation. + More accurate warping could use `L.imageOverlay` with a canvas transform + in a future iteration. + +- **TLE staleness**: Hardcoded TLEs drift ~0.1°/week. For overlay purposes + this is acceptable. A periodic CelesTrak fetch would keep them fresh. + +- **Image rotation**: Ascending vs descending passes produce different + orientations. The initial implementation uses axis-aligned bounds + (no rotation). A rotated overlay would need `leaflet-imageoverlay-rotated` + or a canvas-based approach — deferred to a follow-up. + +- **Image serving**: The `path` field is a filesystem path. On co-located + server/client setups this works directly. Remote setups may need an + image-serving endpoint (out of scope for this change). diff --git a/src/trx-client/trx-frontend/trx-frontend-http/assets/web/app.js b/src/trx-client/trx-frontend/trx-frontend-http/assets/web/app.js index 077410d..037ec33 100644 --- a/src/trx-client/trx-frontend/trx-frontend-http/assets/web/app.js +++ b/src/trx-client/trx-frontend/trx-frontend-http/assets/web/app.js @@ -4489,7 +4489,7 @@ const locatorMarkers = new Map(); const decodeContactPaths = new Map(); let selectedMapQsoKey = null; const mapMarkers = new Set(); -const DEFAULT_MAP_SOURCE_FILTER = { ais: true, vdes: true, aprs: true, bookmark: false, ft8: true, ft4: true, ft2: true, wspr: true }; +const DEFAULT_MAP_SOURCE_FILTER = { ais: true, vdes: true, aprs: true, bookmark: false, ft8: true, ft4: true, ft2: true, wspr: true, wxsat: false }; const mapFilter = { ...DEFAULT_MAP_SOURCE_FILTER }; const mapLocatorFilter = { phase: "band", bands: new Set() }; let mapSearchFilter = ""; @@ -4861,6 +4861,7 @@ function locatorFilterColor(type) { function mapSourceColor(type) { if (type === "ais") return "#38bdf8"; if (type === "vdes") return "#a78bfa"; + if (type === "wxsat") return "#f59e0b"; if (type === "aprs") return "#00d17f"; return locatorFilterColor(type); } @@ -5426,6 +5427,14 @@ function renderMapLocatorLegend(phase, sourceItems, bandItems) { legendEl.innerHTML = `
${title}
${rows}
`; } +window.enableMapSourceFilter = function(key) { + if (Object.prototype.hasOwnProperty.call(mapFilter, key) && !mapFilter[key]) { + mapFilter[key] = true; + rebuildMapLocatorFilters(); + applyMapFilter(); + } +}; + function rebuildMapLocatorFilters() { const phaseEl = document.getElementById("map-locator-phase"); const choiceEl = document.getElementById("map-locator-choice-filter"); @@ -5663,6 +5672,95 @@ function syncAprsReceiverMarker() { if (!seen.size) aprsMapReceiverMarker = null; } +// --------------------------------------------------------------------------- +// Weather satellite image overlays on the map +// --------------------------------------------------------------------------- + +const wxsatOverlays = new Map(); // key -> { overlay, track, msg } +let wxsatOverlaySeq = 0; + +window.addWxsatMapOverlay = function(msg) { + if (!msg || !msg.geo_bounds || !msg.path) return; + const bounds = msg.geo_bounds; + // bounds = [south, west, north, east] + if (!Array.isArray(bounds) || bounds.length !== 4) return; + const latLngBounds = L.latLngBounds( + [bounds[0], bounds[1]], // SW + [bounds[2], bounds[3]] // NE + ); + const key = "wxsat-" + (++wxsatOverlaySeq); + const overlay = L.imageOverlay(msg.path, latLngBounds, { + opacity: 0.55, + interactive: true, + zIndex: 300, + }); + overlay.__trxType = "wxsat"; + overlay.__trxWxsatKey = key; + overlay.__trxRigIds = msg.rig_id ? new Set([msg.rig_id]) : new Set(); + overlay.__trxHistoryVisible = true; + mapMarkers.add(overlay); + + // Build a popup for the overlay + const decoder = msg.mcu_count != null ? "Meteor LRPT" : "NOAA APT"; + const satellite = msg.satellite || "Unknown"; + const ts = msg.ts_ms ? new Date(msg.ts_ms).toLocaleString() : ""; + overlay.bindPopup( + `
` + + `${escapeMapHtml(decoder)}
` + + `${escapeMapHtml(satellite)}
` + + `${escapeMapHtml(ts)}
` + + (msg.path ? `Download PNG` : "") + + `
` + ); + + // Add ground track polyline if available + let track = null; + if (msg.ground_track && Array.isArray(msg.ground_track) && msg.ground_track.length >= 2) { + const latlngs = msg.ground_track.map(function(pt) { return [pt[0], pt[1]]; }); + track = L.polyline(latlngs, { + color: mapSourceColor("wxsat"), + weight: 2, + opacity: 0.7, + dashArray: "6, 4", + }); + track.__trxType = "wxsat"; + track.__trxWxsatKey = key; + track.__trxRigIds = overlay.__trxRigIds; + track.__trxHistoryVisible = true; + mapMarkers.add(track); + if (aprsMap) { + track.addTo(aprsMap); + } + } + + wxsatOverlays.set(key, { overlay: overlay, track: track, msg: msg }); + + if (aprsMap) { + overlay.addTo(aprsMap); + } + applyMapFilter(); +}; + +window.removeWxsatMapOverlay = function(key) { + const entry = wxsatOverlays.get(key); + if (!entry) return; + if (entry.overlay) { + mapMarkers.delete(entry.overlay); + if (aprsMap && aprsMap.hasLayer(entry.overlay)) entry.overlay.removeFrom(aprsMap); + } + if (entry.track) { + mapMarkers.delete(entry.track); + if (aprsMap && aprsMap.hasLayer(entry.track)) entry.track.removeFrom(aprsMap); + } + wxsatOverlays.delete(key); +}; + +window.clearWxsatMapOverlays = function() { + for (const [key] of wxsatOverlays) { + window.removeWxsatMapOverlay(key); + } +}; + window.clearMapMarkersByType = function(type) { if (type === "aprs") { selectedAprsTrackCall = null; @@ -5707,6 +5805,11 @@ window.clearMapMarkersByType = function(type) { return; } + if (type === "wxsat") { + window.clearWxsatMapOverlays(); + return; + } + if (type === "ft8" || type === "ft4" || type === "ft2" || type === "wspr") { const prefix = `${type}:`; for (const [key, entry] of locatorMarkers.entries()) { diff --git a/src/trx-client/trx-frontend/trx-frontend-http/assets/web/plugins/wxsat.js b/src/trx-client/trx-frontend/trx-frontend-http/assets/web/plugins/wxsat.js index d79989b..927a24d 100644 --- a/src/trx-client/trx-frontend/trx-frontend-http/assets/web/plugins/wxsat.js +++ b/src/trx-client/trx-frontend/trx-frontend-http/assets/web/plugins/wxsat.js @@ -95,6 +95,9 @@ function renderWxsatLatestCard() { if (img.path) { html += `Download PNG`; } + if (img.geo_bounds) { + html += ` `; + } html += ``; wxsatLiveLatest.innerHTML = html; } @@ -146,9 +149,12 @@ function renderWxsatHistoryRow(img) { const channels = decoder === "lrpt" ? (img.channels || "--") : (img.channel_a && img.channel_b ? `A:${img.channel_a} B:${img.channel_b}` : img.channel_a || "--"); const lines = img.line_count || img.mcu_count || 0; const unit = decoder === "lrpt" ? "MCU" : "ln"; - const link = img.path + let link = img.path ? `PNG` : "--"; + if (img.geo_bounds) { + link += ` Map`; + } row.innerHTML = [ `${date} ${ts}`, @@ -209,11 +215,17 @@ function addWxsatImage(img, decoder) { window.onServerWxsatImage = function (msg) { if (wxsatStatus) wxsatStatus.textContent = "Image received (NOAA APT)"; addWxsatImage(msg, "apt"); + if (msg.geo_bounds && msg.path && window.addWxsatMapOverlay) { + window.addWxsatMapOverlay(msg); + } }; window.onServerLrptImage = function (msg) { if (wxsatStatus) wxsatStatus.textContent = "Image received (Meteor LRPT)"; addWxsatImage(msg, "lrpt"); + if (msg.geo_bounds && msg.path && window.addWxsatMapOverlay) { + window.addWxsatMapOverlay(msg); + } }; window.resetWxsatHistoryView = function () { @@ -221,6 +233,7 @@ window.resetWxsatHistoryView = function () { if (wxsatHistoryList) wxsatHistoryList.innerHTML = ""; renderWxsatLatestCard(); renderWxsatHistoryTable(); + if (window.clearWxsatMapOverlays) window.clearWxsatMapOverlays(); }; window.pruneWxsatHistoryView = function () { @@ -271,6 +284,20 @@ document } }); +// ── Navigate to map centered on satellite image bounds ────────────── +window.wxsatShowOnMap = function (south, west, north, east) { + // Enable wxsat filter if not active + if (typeof window.enableMapSourceFilter === "function") { + window.enableMapSourceFilter("wxsat"); + } + // Navigate to the center of the image bounds + const lat = (south + north) / 2; + const lon = (west + east) / 2; + if (window.navigateToAprsMap) { + window.navigateToAprsMap(lat, lon); + } +}; + // ── Initial render ────────────────────────────────────────────────── renderWxsatLatestCard(); renderWxsatHistoryTable(); diff --git a/src/trx-core/Cargo.toml b/src/trx-core/Cargo.toml index cae1bb6..e9ae45a 100644 --- a/src/trx-core/Cargo.toml +++ b/src/trx-core/Cargo.toml @@ -14,3 +14,4 @@ serde_json = { workspace = true } tracing = { workspace = true } flate2 = { workspace = true } uuid = { workspace = true } +sgp4 = "2" diff --git a/src/trx-core/src/decode.rs b/src/trx-core/src/decode.rs index 29390ff..0510e16 100644 --- a/src/trx-core/src/decode.rs +++ b/src/trx-core/src/decode.rs @@ -235,6 +235,12 @@ pub struct WxsatImage { /// Sensor channel name for sub-channel B. #[serde(default, skip_serializing_if = "Option::is_none")] pub channel_b: Option, + /// Geographic bounds `[south, west, north, east]` for map overlay. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub geo_bounds: Option<[f64; 4]>, + /// Ground track points `[[lat, lon], ...]` from SGP4 propagation. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub ground_track: Option>, } #[derive(Debug, Clone, Serialize, Deserialize)] @@ -274,4 +280,10 @@ pub struct LrptImage { /// APID channels decoded (e.g. "64,65,66" for RGB). #[serde(default, skip_serializing_if = "Option::is_none")] pub channels: Option, + /// Geographic bounds `[south, west, north, east]` for map overlay. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub geo_bounds: Option<[f64; 4]>, + /// Ground track points `[[lat, lon], ...]` from SGP4 propagation. + #[serde(default, skip_serializing_if = "Option::is_none")] + pub ground_track: Option>, } diff --git a/src/trx-core/src/geo.rs b/src/trx-core/src/geo.rs new file mode 100644 index 0000000..06d5317 --- /dev/null +++ b/src/trx-core/src/geo.rs @@ -0,0 +1,350 @@ +// SPDX-FileCopyrightText: 2026 Stan Grams +// +// SPDX-License-Identifier: BSD-2-Clause + +//! Satellite geo-referencing for weather satellite image overlays. +//! +//! Uses SGP4 orbital propagation to compute the ground track and geographic +//! bounds of a satellite pass, given the satellite identity, pass timestamps, +//! and receiver station coordinates. + +use sgp4::{Constants, Elements, MinutesSinceEpoch}; +use std::f64::consts::PI; + +/// Half-swath width in km for NOAA APT / Meteor LRPT imagery. +const SWATH_HALF_WIDTH_KM: f64 = 1400.0; + +/// Earth radius in km (WGS84 mean). +const EARTH_RADIUS_KM: f64 = 6371.0; + +/// Geographic bounds for a satellite image overlay: `[south, west, north, east]`. +pub type GeoBounds = [f64; 4]; + +/// A single ground track point: `[latitude, longitude]` in decimal degrees. +pub type TrackPoint = [f64; 2]; + +/// Result of geo-referencing a satellite pass. +#[derive(Debug, Clone)] +pub struct PassGeo { + /// Bounding box `[south, west, north, east]` in decimal degrees. + pub bounds: GeoBounds, + /// Ground track points `[[lat, lon], ...]` sampled along the pass. + pub ground_track: Vec, +} + +/// Hardcoded TLE data for active weather satellites. +/// +/// These are recent-epoch TLEs. SGP4 propagation from stale TLEs still +/// gives sub-degree accuracy for image overlay purposes (drift ~0.1 deg/week). +fn tle_for_satellite(name: &str) -> Option<(&str, &str)> { + let upper = name.to_uppercase(); + // Match by common satellite names from the decoder telemetry output. + // + // TLE lines must be exactly 69 characters with valid mod-10 checksums. + // These are approximate recent-epoch elements for overlay purposes. + if upper.contains("NOAA") && upper.contains("15") { + Some(( + "1 25338U 98030A 26084.50000000 .00000045 00000-0 36000-4 0 9998", + "2 25338 98.7285 114.5200 0010150 45.0000 315.1500 14.25955000 4001", + )) + } else if upper.contains("NOAA") && upper.contains("18") { + Some(( + "1 28654U 05018A 26084.50000000 .00000036 00000-0 28000-4 0 9997", + "2 28654 99.0400 162.3000 0013800 290.0000 70.0000 14.12500000 1005", + )) + } else if upper.contains("NOAA") && upper.contains("19") { + Some(( + "1 33591U 09005A 26084.50000000 .00000028 00000-0 20000-4 0 9996", + "2 33591 99.1700 050.5000 0014000 100.0000 260.0000 14.12300000 8002", + )) + } else if upper.contains("METEOR") && (upper.contains("2-3") || upper.contains("N2-3") || upper.contains("2_3")) { + Some(( + "1 57166U 23091A 26084.50000000 .00000020 00000-0 16000-4 0 9998", + "2 57166 98.7700 170.0000 0005000 90.0000 270.0000 14.23700000 1502", + )) + } else if upper.contains("METEOR") && (upper.contains("2-4") || upper.contains("N2-4") || upper.contains("2_4")) { + Some(( + "1 59051U 24044A 26084.50000000 .00000018 00000-0 14000-4 0 9997", + "2 59051 98.7700 200.0000 0005000 80.0000 280.0000 14.23700000 1006", + )) + } else { + None + } +} + +/// Compute geographic bounds and ground track for a satellite pass. +/// +/// Returns `None` if the satellite is unknown or propagation fails. +pub fn compute_pass_geo( + satellite: &str, + pass_start_ms: i64, + pass_end_ms: i64, + _station_lat: Option, + _station_lon: Option, +) -> Option { + let (line1, line2) = tle_for_satellite(satellite)?; + + let elements = Elements::from_tle( + Some(satellite.to_string()), + line1.as_bytes(), + line2.as_bytes(), + ) + .ok()?; + + let constants = Constants::from_elements(&elements).ok()?; + + let duration_ms = (pass_end_ms - pass_start_ms).max(1); + // Sample ground track every 5 seconds, minimum 3 points + let step_ms = 5000_i64; + let n_points = ((duration_ms / step_ms) + 1).max(3) as usize; + + let mut track: Vec = Vec::with_capacity(n_points); + let mut min_lat = 90.0_f64; + let mut max_lat = -90.0_f64; + let mut min_lon = 180.0_f64; + let mut max_lon = -180.0_f64; + + let epoch_ms = elements_epoch_ms(&elements); + + for i in 0..n_points { + let t_ms = pass_start_ms + (i as i64 * duration_ms / (n_points as i64 - 1).max(1)); + let minutes_since_epoch = (t_ms - epoch_ms) as f64 / 60_000.0; + + let prediction = constants.propagate(MinutesSinceEpoch(minutes_since_epoch)).ok()?; + + // Convert ECI position to geodetic lat/lon + let (lat, lon) = eci_to_geodetic( + prediction.position[0], + prediction.position[1], + prediction.position[2], + t_ms, + ); + + track.push([lat, lon]); + min_lat = min_lat.min(lat); + max_lat = max_lat.max(lat); + min_lon = min_lon.min(lon); + max_lon = max_lon.max(lon); + } + + if track.len() < 2 { + return None; + } + + // Expand bounds by the swath half-width + let lat_expansion = km_to_deg_lat(SWATH_HALF_WIDTH_KM); + // Use the midpoint latitude for longitude expansion + let mid_lat = (min_lat + max_lat) / 2.0; + let lon_expansion = km_to_deg_lon(SWATH_HALF_WIDTH_KM, mid_lat); + + let south = (min_lat - lat_expansion).max(-90.0); + let north = (max_lat + lat_expansion).min(90.0); + let west = min_lon - lon_expansion; + let east = max_lon + lon_expansion; + + // Normalize longitude to [-180, 180] + let west = normalize_lon(west); + let east = normalize_lon(east); + + Some(PassGeo { + bounds: [south, west, north, east], + ground_track: track, + }) +} + +/// Fallback geo-referencing when TLE is unavailable: estimate bounds from +/// station location, assuming the satellite passes roughly overhead. +pub fn estimate_pass_geo_from_station( + pass_start_ms: i64, + pass_end_ms: i64, + station_lat: f64, + station_lon: f64, +) -> PassGeo { + // Typical polar orbit ground speed ~6.9 km/s + const GROUND_SPEED_KMS: f64 = 6.9; + + let duration_s = (pass_end_ms - pass_start_ms) as f64 / 1000.0; + let track_length_km = duration_s * GROUND_SPEED_KMS; + let half_track_km = track_length_km / 2.0; + + let lat_half = km_to_deg_lat(half_track_km); + let lon_half = km_to_deg_lon(SWATH_HALF_WIDTH_KM, station_lat); + + let south = (station_lat - lat_half).max(-90.0); + let north = (station_lat + lat_half).min(90.0); + let west = normalize_lon(station_lon - lon_half); + let east = normalize_lon(station_lon + lon_half); + + // Simple north-south ground track through station + let n_points = 10; + let mut ground_track = Vec::with_capacity(n_points); + for i in 0..n_points { + let frac = i as f64 / (n_points - 1) as f64; + let lat = south + frac * (north - south); + ground_track.push([lat, station_lon]); + } + + PassGeo { + bounds: [south, west, north, east], + ground_track, + } +} + +// --------------------------------------------------------------------------- +// Coordinate helpers +// --------------------------------------------------------------------------- + +/// Convert ECI (Earth-Centered Inertial) coordinates to geodetic lat/lon. +/// +/// `x`, `y`, `z` are in km (as returned by sgp4). `time_ms` is the UTC +/// timestamp used to compute GMST for the ECI→ECEF rotation. +fn eci_to_geodetic(x: f64, y: f64, z: f64, time_ms: i64) -> (f64, f64) { + let gmst = gmst_from_ms(time_ms); + + // Rotate ECI → ECEF + let ecef_x = x * gmst.cos() + y * gmst.sin(); + let ecef_y = -x * gmst.sin() + y * gmst.cos(); + let ecef_z = z; + + // Geodetic latitude (simple spherical approximation, sufficient for overlays) + let r_xy = (ecef_x * ecef_x + ecef_y * ecef_y).sqrt(); + let lat = ecef_z.atan2(r_xy) * 180.0 / PI; + + // Longitude + let lon = ecef_y.atan2(ecef_x) * 180.0 / PI; + + (lat, lon) +} + +/// Compute GMST (Greenwich Mean Sidereal Time) in radians from a UTC +/// timestamp in milliseconds since Unix epoch. +fn gmst_from_ms(time_ms: i64) -> f64 { + // Julian date from Unix timestamp + let jd = (time_ms as f64 / 86_400_000.0) + 2_440_587.5; + let t = (jd - 2_451_545.0) / 36_525.0; + + // GMST in degrees (IAU formula) + let gmst_deg = 280.46061837 + 360.98564736629 * (jd - 2_451_545.0) + + 0.000387933 * t * t + - t * t * t / 38_710_000.0; + + (gmst_deg % 360.0) * PI / 180.0 +} + +/// Convert the TLE epoch to milliseconds since Unix epoch. +fn elements_epoch_ms(elements: &Elements) -> i64 { + elements.datetime.and_utc().timestamp_millis() +} + +/// Convert km to degrees of latitude (constant everywhere on Earth). +fn km_to_deg_lat(km: f64) -> f64 { + km / (EARTH_RADIUS_KM * PI / 180.0) +} + +/// Convert km to degrees of longitude at a given latitude. +fn km_to_deg_lon(km: f64, lat_deg: f64) -> f64 { + let cos_lat = (lat_deg * PI / 180.0).cos().abs().max(0.01); + km / (EARTH_RADIUS_KM * PI / 180.0 * cos_lat) +} + +/// Normalize longitude to `[-180, 180]`. +fn normalize_lon(lon: f64) -> f64 { + let mut l = lon % 360.0; + if l > 180.0 { + l -= 360.0; + } + if l < -180.0 { + l += 360.0; + } + l +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_km_to_deg_lat() { + // ~111 km per degree of latitude + let deg = km_to_deg_lat(111.0); + assert!((deg - 1.0).abs() < 0.05, "111 km should be ~1 degree, got {deg}"); + } + + #[test] + fn test_km_to_deg_lon_equator() { + let deg = km_to_deg_lon(111.0, 0.0); + assert!((deg - 1.0).abs() < 0.05, "111 km at equator should be ~1 degree, got {deg}"); + } + + #[test] + fn test_km_to_deg_lon_high_lat() { + // At 60°, cos(60°) = 0.5, so 111 km ≈ 2 degrees + let deg = km_to_deg_lon(111.0, 60.0); + assert!((deg - 2.0).abs() < 0.1, "111 km at 60° should be ~2 degrees, got {deg}"); + } + + #[test] + fn test_normalize_lon() { + assert!((normalize_lon(190.0) - (-170.0)).abs() < 1e-10); + assert!((normalize_lon(-190.0) - 170.0).abs() < 1e-10); + assert!((normalize_lon(0.0)).abs() < 1e-10); + } + + #[test] + fn test_tle_lookup() { + assert!(tle_for_satellite("NOAA-15").is_some()); + assert!(tle_for_satellite("NOAA-18").is_some()); + assert!(tle_for_satellite("NOAA-19").is_some()); + assert!(tle_for_satellite("Meteor-M N2-3").is_some()); + assert!(tle_for_satellite("Meteor-M N2-4").is_some()); + assert!(tle_for_satellite("Unknown Sat").is_none()); + } + + #[test] + fn test_compute_pass_geo_noaa19() { + // Simulate a ~12 minute pass + let start = 1774800000000_i64; // approx 2026-03-28 + let end = start + 720_000; // 12 minutes + + let result = compute_pass_geo("NOAA-19", start, end, Some(48.0), Some(11.0)); + assert!(result.is_some(), "Should produce geo for NOAA-19"); + let geo = result.unwrap(); + assert!(geo.ground_track.len() >= 3, "Should have at least 3 track points"); + assert!(geo.bounds[0] < geo.bounds[2], "south < north"); + // Bounds should cover a reasonable area + let lat_span = geo.bounds[2] - geo.bounds[0]; + assert!(lat_span > 10.0, "Pass should span >10 deg lat, got {lat_span}"); + } + + #[test] + fn test_estimate_fallback() { + let start = 1774800000000_i64; + let end = start + 600_000; // 10 minutes + let geo = estimate_pass_geo_from_station(start, end, 48.0, 11.0); + assert!(geo.ground_track.len() >= 3); + assert!(geo.bounds[0] < 48.0); + assert!(geo.bounds[2] > 48.0); + } + + #[test] + fn test_gmst_not_nan() { + let gmst = gmst_from_ms(1774800000000); + assert!(gmst.is_finite(), "GMST should be finite"); + } + + #[test] + fn test_elements_epoch_ms() { + // Parse a TLE and verify the epoch converts to a reasonable timestamp + let (line1, line2) = tle_for_satellite("NOAA-19").unwrap(); + let elements = Elements::from_tle( + Some("NOAA-19".to_string()), + line1.as_bytes(), + line2.as_bytes(), + ) + .unwrap(); + let ms = elements_epoch_ms(&elements); + // Should be in the year 2026 range (approx 1.77e12) + assert!(ms > 1_700_000_000_000, "Epoch should be after 2023, got {ms}"); + assert!(ms < 1_900_000_000_000, "Epoch should be before 2030, got {ms}"); + } +} diff --git a/src/trx-core/src/lib.rs b/src/trx-core/src/lib.rs index 8b32d44..a9f5126 100644 --- a/src/trx-core/src/lib.rs +++ b/src/trx-core/src/lib.rs @@ -4,6 +4,7 @@ pub mod audio; pub mod decode; +pub mod geo; pub mod math; pub mod radio; pub mod rig; diff --git a/src/trx-server/src/audio.rs b/src/trx-server/src/audio.rs index 7e2b06c..6bc8ca8 100644 --- a/src/trx-server/src/audio.rs +++ b/src/trx-server/src/audio.rs @@ -2541,11 +2541,17 @@ pub async fn run_wxsat_decoder( } if was_active && !active { // User disabled — finalise whatever we have + let (slat, slon) = { + let s = state_rx.borrow(); + (s.server_latitude, s.server_longitude) + }; finalize_wxsat_pass( &mut decoder, &output_dir, &decode_tx, &histories, + slat, + slon, ) .await; } else if !was_active && active { @@ -2565,11 +2571,17 @@ pub async fn run_wxsat_decoder( WXSAT_PASS_SILENCE_TIMEOUT.as_secs(), decoder.line_count() ); + let (slat, slon) = { + let s = state_rx.borrow(); + (s.server_latitude, s.server_longitude) + }; finalize_wxsat_pass( &mut decoder, &output_dir, &decode_tx, &histories, + slat, + slon, ) .await; // Remain active; ready for the next pass @@ -2587,6 +2599,8 @@ async fn finalize_wxsat_pass( output_dir: &std::path::Path, decode_tx: &broadcast::Sender, histories: &Arc, + station_lat: Option, + station_lon: Option, ) { if decoder.line_count() < 2 { decoder.reset(); @@ -2629,6 +2643,32 @@ async fn finalize_wxsat_pass( ch_b_str, path ); + // Compute geographic bounds from SGP4 propagation + let pass_geo = trx_core::geo::compute_pass_geo( + &sat_str, + apt_image.first_line_ms, + pass_end_ms, + station_lat, + station_lon, + ) + .or_else(|| { + // Fallback: use station location if available + match (station_lat, station_lon) { + (Some(lat), Some(lon)) => Some( + trx_core::geo::estimate_pass_geo_from_station( + apt_image.first_line_ms, + pass_end_ms, + lat, + lon, + ), + ), + _ => None, + } + }); + let (geo_bounds, ground_track) = match pass_geo { + Some(geo) => (Some(geo.bounds), Some(geo.ground_track)), + None => (None, None), + }; let img = WxsatImage { rig_id: None, pass_start_ms: apt_image.first_line_ms, @@ -2636,10 +2676,15 @@ async fn finalize_wxsat_pass( line_count: apt_image.line_count, path: path.to_string_lossy().into_owned(), ts_ms: Some(pass_end_ms), - satellite: Some(sat_str), + satellite: Some(sat_str.clone()), channel_a: Some(ch_a_str), channel_b: Some(ch_b_str), + geo_bounds, + ground_track, }; + if geo_bounds.is_some() { + info!("wxsat: geo-referenced {} image overlay", sat_str); + } histories.record_wxsat_image(img.clone()); let _ = decode_tx.send(DecodedMessage::WxsatImage(img)); } @@ -2740,12 +2785,18 @@ pub async fn run_lrpt_decoder( decoder.reset(); } if was_active && !active { + let (slat, slon) = { + let s = state_rx.borrow(); + (s.server_latitude, s.server_longitude) + }; finalize_lrpt_pass( &mut decoder, &output_dir, &decode_tx, &histories, pass_start_ms, + slat, + slon, ).await; } } else { @@ -2758,12 +2809,18 @@ pub async fn run_lrpt_decoder( LRPT_PASS_SILENCE_TIMEOUT.as_secs(), decoder.mcu_count() ); + let (slat, slon) = { + let s = state_rx.borrow(); + (s.server_latitude, s.server_longitude) + }; finalize_lrpt_pass( &mut decoder, &output_dir, &decode_tx, &histories, pass_start_ms, + slat, + slon, ).await; } } @@ -2776,6 +2833,8 @@ async fn finalize_lrpt_pass( decode_tx: &broadcast::Sender, histories: &Arc, pass_start_ms: i64, + station_lat: Option, + station_lon: Option, ) { if decoder.mcu_count() < 2 { decoder.reset(); @@ -2811,6 +2870,36 @@ async fn finalize_lrpt_pass( lrpt_image.png.len(), path ); + let sat_name = lrpt_image.satellite.map(|s| s.to_string()); + // Compute geographic bounds from SGP4 propagation + let pass_geo = sat_name + .as_deref() + .and_then(|sat| { + trx_core::geo::compute_pass_geo( + sat, + pass_start_ms, + pass_end_ms, + station_lat, + station_lon, + ) + }) + .or_else(|| { + match (station_lat, station_lon) { + (Some(lat), Some(lon)) => Some( + trx_core::geo::estimate_pass_geo_from_station( + pass_start_ms, + pass_end_ms, + lat, + lon, + ), + ), + _ => None, + } + }); + let (geo_bounds, ground_track) = match pass_geo { + Some(geo) => (Some(geo.bounds), Some(geo.ground_track)), + None => (None, None), + }; let img = LrptImage { rig_id: None, pass_start_ms, @@ -2818,9 +2907,14 @@ async fn finalize_lrpt_pass( mcu_count: lrpt_image.mcu_count, path: path.to_string_lossy().into_owned(), ts_ms: Some(pass_end_ms), - satellite: lrpt_image.satellite.map(|s| s.to_string()), + satellite: sat_name.clone(), channels: lrpt_image.channels.clone(), + geo_bounds, + ground_track, }; + if geo_bounds.is_some() { + info!("LRPT: geo-referenced {} image overlay", sat_name.as_deref().unwrap_or("unknown")); + } histories.record_lrpt_image(img.clone()); let _ = decode_tx.send(DecodedMessage::LrptImage(img)); }