diff --git a/assets/js/propagation_map_hook.js b/assets/js/propagation_map_hook.js index bcef7b44..3c51e960 100644 --- a/assets/js/propagation_map_hook.js +++ b/assets/js/propagation_map_hook.js @@ -358,11 +358,43 @@ function buildPopupHTML(detail, viewshedLoading) { ${explanations} ${detail.factors.duct_info ? buildDuctInfoHTML(detail.factors.duct_info) : ""} + ${detail.rain_scatter && detail.rain_scatter.cells.length > 0 ? buildScatterHTML(detail.rain_scatter) : ""} ${detail.forecast ? buildForecastSvg(detail.forecast) : ""} ${viewshedStatus} ` } +function buildScatterHTML(scatter) { + const cls = scatter.classification + const cells = scatter.cells + const clsColors = { excellent: "#059669", good: "#0d9488", marginal: "#ca8a04", none: "#666" } + const clsColor = clsColors[cls] || "#666" + const clsLabel = cls.charAt(0).toUpperCase() + cls.slice(1) + + const top3 = cells.slice(0, 3) + const rows = top3.map(c => { + const dir = bearingLabel(c.bearing) + return `
+ ${c.dbz} dBZ at ${Math.round(c.distance_km)} km ${dir} (${c.scatter_db} dB) +
` + }).join("") + + return ` +
+
+ Rain Scatter + ${clsLabel} +
+ ${rows} + ${cells.length > 3 ? `
+${cells.length - 3} more cells
` : ""} +
` +} + +function bearingLabel(deg) { + const dirs = ["N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW"] + return dirs[Math.round(deg / 22.5) % 16] +} + function buildDuctInfoHTML(info) { const layers = info.ducts || [] let layerRows = "" @@ -530,6 +562,29 @@ export const PropagationMap = { this.detailPanel.innerHTML = buildPopupHTML(merged, this.viewshedLoading) this.showDetailPanel() + // Draw rain scatter cells on the map + if (this.scatterMarkers) { + this.scatterMarkers.clearLayers() + } else { + this.scatterMarkers = L.layerGroup().addTo(this.map) + } + if (detail.rain_scatter && detail.rain_scatter.cells.length > 0) { + for (const c of detail.rain_scatter.cells) { + const opacity = Math.min(0.9, Math.max(0.3, (c.scatter_db + 30) / 30)) + const color = c.dbz >= 45 ? "#dc2626" : c.dbz >= 35 ? "#ea580c" : "#ca8a04" + const marker = L.circleMarker([c.lat, c.lon], { + radius: Math.min(12, Math.max(4, c.dbz / 5)), + color: color, + fillColor: color, + fillOpacity: opacity * 0.5, + weight: 1.5, + opacity: opacity, + interactive: false + }) + this.scatterMarkers.addLayer(marker) + } + } + // Draw propagation reach polygon based on contiguous good cells if (this.gridLookup && this.clickedLatLng) { // Use MARGINAL threshold (50) as minimum for propagation reach @@ -677,6 +732,7 @@ export const PropagationMap = { this.detailPanel.innerHTML = "" } this.rangeCircles.clearLayers() + if (this.scatterMarkers) this.scatterMarkers.clearLayers() this.clickedLatLng = null this.lastDetail = null }, diff --git a/lib/microwaveprop/propagation/rain_scatter.ex b/lib/microwaveprop/propagation/rain_scatter.ex new file mode 100644 index 00000000..ffcda330 --- /dev/null +++ b/lib/microwaveprop/propagation/rain_scatter.ex @@ -0,0 +1,143 @@ +defmodule Microwaveprop.Propagation.RainScatter do + @moduledoc """ + Estimates rain scatter potential from NEXRAD composite reflectivity. + + Rain scatter enables microwave contacts at 100-300+ km by scattering + signals off precipitation cells. Signal strength depends on reflectivity + (rain intensity), frequency, and geometry (distance from each station + to the rain cell). + + Uses a simplified bistatic radar equation: + scatter_db ≈ 10*log10(Z) + 10*log10(V) - 20*log10(R1) - 20*log10(R2) + + frequency_gain - path_losses + + where Z is reflectivity factor (from dBZ), V is effective scattering + volume, and R1/R2 are distances from each endpoint to the cell. + """ + + @earth_radius_km 6371.0 + + # Minimum reflectivity to consider for scatter (dBZ) + @min_dbz 25.0 + + # Maximum scatter range from either endpoint (km) + @max_range_km 300.0 + + @doc """ + Find rain cells with scatter potential for a given point and band. + + Takes a list of `{lat, lon, dbz}` rain cells (from NEXRAD extraction), + the observer's position, and frequency in GHz. + + Returns a list of scatter cell maps sorted by potential (strongest first), + each with: lat, lon, dbz, distance_km, scatter_db (relative signal estimate), + and bearing from the observer. + """ + def find_scatter_cells(rain_cells, obs_lat, obs_lon, freq_ghz) do + rain_cells + |> Enum.filter(fn {_lat, _lon, dbz} -> dbz >= @min_dbz end) + |> Enum.map(fn {lat, lon, dbz} -> + dist_km = haversine_km(obs_lat, obs_lon, lat, lon) + bearing = bearing_deg(obs_lat, obs_lon, lat, lon) + + # Scatter signal estimate (relative dB) + # Higher reflectivity = more scattering targets + # Closer cells = stronger signal (inverse square from both endpoints) + # Higher frequency = stronger Rayleigh scattering (up to ~10 GHz, then Mie) + scatter_db = estimate_scatter_db(dbz, dist_km, freq_ghz) + + %{ + lat: Float.round(lat, 3), + lon: Float.round(lon, 3), + dbz: Float.round(dbz, 1), + distance_km: Float.round(dist_km, 1), + bearing: Float.round(bearing, 0), + scatter_db: Float.round(scatter_db, 1) + } + end) + |> Enum.filter(&(&1.distance_km <= @max_range_km and &1.distance_km >= 10)) + |> Enum.sort_by(&(-&1.scatter_db)) + |> Enum.take(20) + end + + @doc """ + Classify overall scatter potential for a point. + + Returns :excellent, :good, :marginal, or :none based on + the best available scatter cell. + """ + def classify(scatter_cells) do + case scatter_cells do + [] -> :none + [best | _] -> + cond do + best.scatter_db >= -10 -> :excellent + best.scatter_db >= -20 -> :good + best.scatter_db >= -30 -> :marginal + true -> :none + end + end + end + + # Simplified scatter signal estimate in relative dB. + # + # Based on bistatic radar equation for volume scattering: + # P_rx ∝ Z * σ_scatter * V / (R^4) + # + # For a single cell at distance R from the observer (assuming the + # other station is also near R for a rough estimate): + # scatter_db ≈ dBZ + freq_factor - 40*log10(R_km) + volume_term + defp estimate_scatter_db(dbz, dist_km, freq_ghz) do + # Reflectivity contribution (dBZ is already in log scale) + z_term = dbz + + # Frequency factor: Rayleigh scattering ∝ f^4 below ~10 GHz, + # transitions to Mie at higher frequencies (weaker dependence). + # Normalize to 10 GHz as reference. + freq_factor = + cond do + freq_ghz <= 10 -> 40 * :math.log10(max(freq_ghz, 0.5) / 10.0) + freq_ghz <= 50 -> 20 * :math.log10(freq_ghz / 10.0) + true -> 20 * :math.log10(50.0 / 10.0) + end + + # Path loss: R^4 for bistatic (R^2 each way) + # Normalize to 100 km reference distance + range_loss = 40 * :math.log10(max(dist_km, 1) / 100.0) + + # Effective scattering volume (~1 km^3 rain cell) + volume_term = 10.0 + + # Baseline offset to center the scale around useful values + baseline = -50.0 + + baseline + z_term + freq_factor - range_loss + volume_term + end + + defp haversine_km(lat1, lon1, lat2, lon2) do + dlat = :math.pi() * (lat2 - lat1) / 180.0 + dlon = :math.pi() * (lon2 - lon1) / 180.0 + rlat1 = :math.pi() * lat1 / 180.0 + rlat2 = :math.pi() * lat2 / 180.0 + + a = + :math.sin(dlat / 2) * :math.sin(dlat / 2) + + :math.cos(rlat1) * :math.cos(rlat2) * + :math.sin(dlon / 2) * :math.sin(dlon / 2) + + c = 2.0 * :math.atan2(:math.sqrt(a), :math.sqrt(1.0 - a)) + @earth_radius_km * c + end + + defp bearing_deg(lat1, lon1, lat2, lon2) do + rlat1 = :math.pi() * lat1 / 180.0 + rlat2 = :math.pi() * lat2 / 180.0 + dlon = :math.pi() * (lon2 - lon1) / 180.0 + + x = :math.sin(dlon) * :math.cos(rlat2) + y = :math.cos(rlat1) * :math.sin(rlat2) - :math.sin(rlat1) * :math.cos(rlat2) * :math.cos(dlon) + + bearing = :math.atan2(x, y) * 180.0 / :math.pi() + Float.round(:math.fmod(bearing + 360.0, 360.0), 1) + end +end diff --git a/lib/microwaveprop/weather/nexrad_client.ex b/lib/microwaveprop/weather/nexrad_client.ex index 199df5aa..3e8b5b7c 100644 --- a/lib/microwaveprop/weather/nexrad_client.ex +++ b/lib/microwaveprop/weather/nexrad_client.ex @@ -58,6 +58,78 @@ defmodule Microwaveprop.Weather.NexradClient do "#{@base_url}/#{date_path}/GIS/uscomp/n0q_#{file_ts}.png" end + @doc """ + Fetch the latest n0q frame and extract all rain cells with reflectivity + above `min_dbz` within `radius_km` of the given point. Returns + `{:ok, [{lat, lon, dbz}]}` or `{:error, reason}`. + + Samples every ~5 km (10 pixels) within the bounding box for efficiency. + """ + @spec fetch_rain_cells(float(), float(), float(), float()) :: {:ok, [{float(), float(), float()}]} | {:error, term()} + def fetch_rain_cells(lat, lon, radius_km \\ 300.0, min_dbz \\ 25.0) do + now = DateTime.utc_now() + rounded = round_to_5min(now) + url = frame_url(rounded) + + req_opts = Application.get_env(:microwaveprop, :nexrad_req_options, []) + + case Req.get(url, [receive_timeout: 60_000, retry: false] ++ req_opts) do + {:ok, %{status: 200, body: body}} when is_binary(body) -> + case decode_png_to_pixels(body) do + {:ok, pixels, width} -> + cells = extract_rain_cells(pixels, width, lat, lon, radius_km, min_dbz) + {:ok, cells} + + {:error, reason} -> + {:error, reason} + end + + {:ok, %{status: status}} -> + {:error, "NEXRAD n0q HTTP #{status}"} + + {:error, reason} -> + {:error, reason} + end + end + + defp extract_rain_cells(pixels, width, center_lat, center_lon, radius_km, min_dbz) do + height = div(byte_size(pixels), width) + + # Convert radius to approximate pixel count (~0.005 deg/px, ~0.5 km/px at mid-lat) + dlat = radius_km / 111.0 + dlon = radius_km / (111.0 * :math.cos(:math.pi() * center_lat / 180.0)) + + lat_min = center_lat - dlat + lat_max = center_lat + dlat + lon_min = center_lon - dlon + lon_max = center_lon + dlon + + {x_min, y_max_px} = latlon_to_pixel(lat_min, lon_min) + {x_max, y_min_px} = latlon_to_pixel(lat_max, lon_max) + + x_min = max(x_min, 0) + x_max = min(x_max, width - 1) + y_min_px = max(y_min_px, 0) + y_max_px = min(y_max_px, height - 1) + + # Sample every 10 pixels (~5 km) for efficiency + step = 10 + + for y <- y_min_px..y_max_px//step, + x <- x_min..x_max//step, + offset = y * width + x, + offset >= 0, + offset < byte_size(pixels), + <<_::binary-size(offset), pixel_val::8, _::binary>> = pixels, + pixel_val > 0, + dbz = pixel_to_dbz(pixel_val), + dbz >= min_dbz do + cell_lat = @lat_max - y * @deg_per_pixel + cell_lon = @lon_min + x * @deg_per_pixel + {Float.round(cell_lat, 3), Float.round(cell_lon, 3), Float.round(dbz, 1)} + end + end + @doc """ Convert (lat, lon) to pixel coordinates in the n0q image. diff --git a/lib/microwaveprop_web/live/map_live.ex b/lib/microwaveprop_web/live/map_live.ex index 8ec6dce8..59a10c8d 100644 --- a/lib/microwaveprop_web/live/map_live.ex +++ b/lib/microwaveprop_web/live/map_live.ex @@ -134,7 +134,9 @@ defmodule MicrowavepropWeb.MapLive do %{time: DateTime.to_iso8601(f.valid_time), score: f.score} end) - payload = if detail, do: Map.put(detail, :forecast, forecast_data), else: %{} + scatter = fetch_rain_scatter(lat, lon, band) + + payload = if detail, do: detail |> Map.put(:forecast, forecast_data) |> Map.put(:rain_scatter, scatter), else: %{} {:noreply, push_event(socket, "point_detail", payload)} end @@ -253,6 +255,23 @@ defmodule MicrowavepropWeb.MapLive do end end + defp fetch_rain_scatter(lat, lon, band_mhz) do + alias Microwaveprop.Propagation.RainScatter + alias Microwaveprop.Weather.NexradClient + + freq_ghz = band_mhz / 1000.0 + + case NexradClient.fetch_rain_cells(lat, lon) do + {:ok, rain_cells} -> + cells = RainScatter.find_scatter_cells(rain_cells, lat, lon, freq_ghz) + classification = RainScatter.classify(cells) + %{cells: cells, classification: to_string(classification)} + + {:error, _} -> + %{cells: [], classification: "none"} + end + end + defp band_info(band_mhz) do config = BandConfig.get(band_mhz)