diff --git a/lib/microwaveprop_web/live/contact_live/show.ex b/lib/microwaveprop_web/live/contact_live/show.ex index db0b8512..bd6618a5 100644 --- a/lib/microwaveprop_web/live/contact_live/show.ex +++ b/lib/microwaveprop_web/live/contact_live/show.ex @@ -15,6 +15,7 @@ defmodule MicrowavepropWeb.ContactLive.Show do alias Microwaveprop.Workers.HrrrFetchWorker alias Microwaveprop.Workers.SolarIndexWorker alias Microwaveprop.Workers.TerrainProfileWorker + alias MicrowavepropWeb.SkewT require Logger @@ -49,7 +50,7 @@ defmodule MicrowavepropWeb.ContactLive.Show do propagation_analysis: nil, data_sources: nil, terrain_expanded: false, - hrrr_profile_expanded: false, + hrrr_profile_expanded: true, obs_sort_by: "station_name", obs_sort_order: "asc", sounding_sort_by: "station_name", @@ -1168,7 +1169,8 @@ defmodule MicrowavepropWeb.ContactLive.Show do <%= if profile.profile && profile.profile != [] do %> -
+ <.skew_t_chart profile={profile.profile} /> +
@@ -2125,4 +2127,176 @@ defmodule MicrowavepropWeb.ContactLive.Show do defp maybe_put_band(params, band) when is_integer(band), do: Map.put(params, "band", Integer.to_string(band)) defp maybe_put_band(params, band) when is_binary(band), do: Map.put(params, "band", band) defp maybe_put_band(params, _), do: params + + # Render a skew-T log-P diagram from a pressure-level profile + # (list of `%{"pres", "tmpc", "dwpc"}` maps). Returns empty output + # when the profile has fewer than three usable levels. + attr :profile, :list, required: true + + defp skew_t_chart(assigns) do + chart = SkewT.build(assigns.profile) + usable_points = if chart, do: length(chart.t_points), else: 0 + assigns = assign(assigns, chart: chart, usable_points: usable_points) + + ~H""" +
= 2} class="mt-3"> +
+ Skew-T log-P +
+ + <%!-- Plot background --%> + + + <%!-- Clip background curves to the plot area so skewed lines + don't bleed over the axes. --%> + + + + + + + + <%!-- Dry adiabats (orange curves) --%> + + + <%!-- Mixing ratio lines (dashed teal) --%> + + + <%!-- Isotherms (gray) --%> + + + <%!-- Isobars (gray, across plot) --%> + + + <%!-- Dewpoint trace --%> + + + + <%!-- Temperature trace --%> + + + + + <%!-- Pressure axis labels (left side) --%> + + {bar.label} + + + <%!-- Temperature axis labels (bottom, only those inside the visible band) --%> + = -40 && iso.t_c <= 40 && rem(iso.t_c, 20) == 0} + x={iso.label_x} + y={iso.label_y} + font-size="10" + text-anchor="middle" + fill="currentColor" + class="text-base-content/60" + > + {iso.t_c}°C + + + <%!-- Legend --%> + + + + Temperature + + + + Dewpoint + + + +
+ Orange curves: dry adiabats · dashed teal: saturation mixing ratio · + thin gray: isotherms and isobars +
+
+ """ + end end diff --git a/lib/microwaveprop_web/skew_t.ex b/lib/microwaveprop_web/skew_t.ex new file mode 100644 index 00000000..48675792 --- /dev/null +++ b/lib/microwaveprop_web/skew_t.ex @@ -0,0 +1,297 @@ +defmodule MicrowavepropWeb.SkewT do + @moduledoc """ + Skew-T log-P diagram renderer. Builds the full set of SVG + primitives (isobars, isotherms, dry adiabats, saturation mixing + ratio lines, and temperature/dewpoint traces) for a single HRRR + or ERA5 pressure-level profile so the contact detail template can + draw it declaratively. + + ## Meteorology + + * Magnus formula for saturation vapor pressure: + `es(T) = 6.112 * exp(17.67 * T / (T + 243.5))` with T in °C. + * Saturation mixing ratio: `ws = 622 * es / (p - es)` in g/kg. + * Dry adiabat (constant potential temperature θ, K): + `T(p) = θ * (p / 1000)^(R/cp)` with `R/cp ≈ 0.2854`. + + ## Skew transform + + Pressure is on a log axis (top to bottom, larger p at the bottom). + Isotherms are skewed 45° by shifting each (T, p) point right as + pressure decreases, so a vertical line on the chart corresponds + to a temperature that rises with height — a hallmark of the + classic skew-T log-P diagram. + """ + + @dew_point_line_color "#16a34a" + @temp_line_color "#dc2626" + + @dry_adiabat_thetas_k 250..360//10 + @mixing_ratios_g_kg [0.4, 1.0, 2.0, 4.0, 8.0, 12.0, 16.0, 20.0] + + @type projected :: {x :: float(), y :: float()} + + @type chart :: %{ + width: number(), + height: number(), + padding_top: number(), + padding_bottom: number(), + padding_left: number(), + padding_right: number(), + plot_width: number(), + plot_height: number(), + t_min: number(), + t_max: number(), + p_bot: number(), + p_top: number(), + skew: number() + } + + @default_opts [ + width: 560, + height: 620, + padding_top: 24, + padding_bottom: 36, + padding_left: 44, + padding_right: 20, + t_min: -40.0, + t_max: 40.0, + p_bot: 1050.0, + p_top: 100.0, + skew: 0.55 + ] + + @doc """ + Build a chart config from overrideable keyword options. Returned + struct has the plot dimensions pre-computed. + """ + @spec chart_config(keyword()) :: chart() + def chart_config(opts \\ []) do + opts = Keyword.merge(@default_opts, opts) + + width = Keyword.fetch!(opts, :width) + height = Keyword.fetch!(opts, :height) + padding_top = Keyword.fetch!(opts, :padding_top) + padding_bottom = Keyword.fetch!(opts, :padding_bottom) + padding_left = Keyword.fetch!(opts, :padding_left) + padding_right = Keyword.fetch!(opts, :padding_right) + + plot_width = width - padding_left - padding_right + plot_height = height - padding_top - padding_bottom + + %{ + width: width, + height: height, + padding_top: padding_top, + padding_bottom: padding_bottom, + padding_left: padding_left, + padding_right: padding_right, + plot_width: plot_width, + plot_height: plot_height, + t_min: Keyword.fetch!(opts, :t_min), + t_max: Keyword.fetch!(opts, :t_max), + p_bot: Keyword.fetch!(opts, :p_bot), + p_top: Keyword.fetch!(opts, :p_top), + skew: Keyword.fetch!(opts, :skew) + } + end + + @doc """ + Build the skew-T primitives for `profile` — a list of maps with + string keys `"pres"`, `"tmpc"`, `"dwpc"`, and optionally `"hght"`. + Returns `nil` when the profile is missing or empty. + """ + @spec build([map()] | nil, keyword()) :: map() | nil + def build(profile, opts \\ []) + def build(nil, _), do: nil + def build([], _), do: nil + + def build(profile, opts) when is_list(profile) do + chart = chart_config(opts) + + t_points = + profile + |> Enum.filter(fn l -> is_number(l["tmpc"]) and is_number(l["pres"]) end) + |> Enum.sort_by(&(-&1["pres"])) + |> Enum.map(fn l -> project(l["tmpc"], l["pres"], chart) end) + + td_points = + profile + |> Enum.filter(fn l -> is_number(l["dwpc"]) and is_number(l["pres"]) end) + |> Enum.sort_by(&(-&1["pres"])) + |> Enum.map(fn l -> project(l["dwpc"], l["pres"], chart) end) + + %{ + width: chart.width, + height: chart.height, + padding_top: chart.padding_top, + padding_bottom: chart.padding_bottom, + padding_left: chart.padding_left, + padding_right: chart.padding_right, + plot_width: chart.plot_width, + plot_height: chart.plot_height, + isobars: isobars(chart), + isotherms: isotherms(chart), + dry_adiabats: dry_adiabats(chart), + mixing_ratio_lines: mixing_ratio_lines(chart), + t_points: t_points, + td_points: td_points, + t_path: points_to_path(t_points), + td_path: points_to_path(td_points), + temp_color: @temp_line_color, + dewpoint_color: @dew_point_line_color + } + end + + # ── Meteorology ────────────────────────────────────────────── + + @doc """ + Saturation vapor pressure (hPa / mb) for temperature `t_c` in °C. + Magnus formula with the Tetens coefficients used by IFS/GFS. + """ + @spec saturation_vapor_pressure(number()) :: float() + def saturation_vapor_pressure(t_c) when is_number(t_c) do + 6.112 * :math.exp(17.67 * t_c / (t_c + 243.5)) + end + + @doc """ + Invert the saturation-mixing-ratio relation: given a mixing + ratio `ws_g_kg` and pressure `p_mb`, return the temperature + (°C) whose saturation vapor pressure yields that mixing ratio. + """ + @spec temperature_from_mixing_ratio(number(), number()) :: float() + def temperature_from_mixing_ratio(ws_g_kg, p_mb) when is_number(ws_g_kg) and is_number(p_mb) and ws_g_kg > 0 do + es = ws_g_kg * p_mb / (622.0 + ws_g_kg) + ln = :math.log(es / 6.112) + 243.5 * ln / (17.67 - ln) + end + + @doc """ + Temperature (°C) along the dry adiabat with potential temperature + `theta_k` (K) at pressure `p_mb` (mb). + """ + @spec dry_adiabat_temperature(number(), number()) :: float() + def dry_adiabat_temperature(theta_k, p_mb) when is_number(theta_k) and is_number(p_mb) do + theta_k * :math.pow(p_mb / 1000.0, 0.2854) - 273.15 + end + + # ── Projection ─────────────────────────────────────────────── + + @doc """ + Project a `(t_c, p_mb)` point onto chart pixel coordinates. + """ + @spec project(number(), number(), chart()) :: projected() + def project(t_c, p_mb, chart) when is_number(t_c) and is_number(p_mb) do + p_ratio = pressure_ratio(p_mb, chart) + + y = chart.padding_top + (1.0 - p_ratio) * chart.plot_height + + x_base = + chart.padding_left + + (t_c - chart.t_min) / (chart.t_max - chart.t_min) * chart.plot_width + + x = x_base + chart.skew * p_ratio * chart.plot_width + {x, y} + end + + defp pressure_ratio(p_mb, chart) do + # 0 at p_bot (chart bottom), 1 at p_top (chart top). Log-P scale. + :math.log(chart.p_bot / p_mb) / :math.log(chart.p_bot / chart.p_top) + end + + # ── Background curves ──────────────────────────────────────── + + @isobar_levels_mb [1000, 850, 700, 500, 400, 300, 250, 200, 150, 100] + + defp isobars(chart) do + for p <- @isobar_levels_mb, p <= chart.p_bot and p >= chart.p_top do + {_x, y} = project(0.0, p * 1.0, chart) + + %{ + y: y, + x_start: chart.padding_left, + x_end: chart.padding_left + chart.plot_width, + label: "#{p}", + label_y: y + 3 + } + end + end + + defp isotherms(chart) do + # Draw isotherms every 10°C across the data area. Each one is a + # straight line between its (t, p_bot) and (t, p_top) projections, + # clipped to the plot rectangle by the svg viewport. + for t_c <- -120..40//10 do + {x_bot, y_bot} = project(t_c * 1.0, chart.p_bot, chart) + {x_top, y_top} = project(t_c * 1.0, chart.p_top, chart) + + label_x = + if t_c >= chart.t_min and t_c <= chart.t_max do + {x, _} = project(t_c * 1.0, chart.p_bot, chart) + x + end + + %{ + x1: x_bot, + y1: y_bot, + x2: x_top, + y2: y_top, + t_c: t_c, + label_x: label_x, + label_y: chart.padding_top + chart.plot_height + 14 + } + end + end + + defp dry_adiabats(chart) do + pressures = log_pressure_samples(chart.p_bot, chart.p_top, 24) + + for theta <- @dry_adiabat_thetas_k do + points = + Enum.map(pressures, fn p -> + t = dry_adiabat_temperature(theta * 1.0, p) + project(t, p, chart) + end) + + %{theta_k: theta, path: points_to_path(points)} + end + end + + defp mixing_ratio_lines(chart) do + # Mixing ratio isohumes only make sense in the lower troposphere. + pressures = log_pressure_samples(chart.p_bot, 400.0, 18) + + for ws <- @mixing_ratios_g_kg do + points = + Enum.map(pressures, fn p -> + t = temperature_from_mixing_ratio(ws, p) + project(t, p, chart) + end) + + %{ws_g_kg: ws, path: points_to_path(points)} + end + end + + defp log_pressure_samples(p_from, p_to, count) do + log_from = :math.log(p_from) + log_to = :math.log(p_to) + step = (log_to - log_from) / (count - 1) + + for i <- 0..(count - 1) do + :math.exp(log_from + i * step) + end + end + + defp points_to_path([]), do: "" + + defp points_to_path([{x0, y0} | rest]) do + initial = "M#{format_coord(x0)},#{format_coord(y0)}" + + body = + Enum.map_join(rest, " ", fn {x, y} -> "L#{format_coord(x)},#{format_coord(y)}" end) + + if body == "", do: initial, else: initial <> " " <> body + end + + defp format_coord(n) when is_number(n), do: :erlang.float_to_binary(n * 1.0, decimals: 2) +end diff --git a/test/microwaveprop_web/skew_t_test.exs b/test/microwaveprop_web/skew_t_test.exs new file mode 100644 index 00000000..0417cc5f --- /dev/null +++ b/test/microwaveprop_web/skew_t_test.exs @@ -0,0 +1,128 @@ +defmodule MicrowavepropWeb.SkewTTest do + use ExUnit.Case, async: true + + alias MicrowavepropWeb.SkewT + + describe "saturation_vapor_pressure/1 (Magnus)" do + test "reference values are within 1% of textbook" do + # At 0°C the Magnus formula gives ≈ 6.112 hPa. + assert_in_delta SkewT.saturation_vapor_pressure(0.0), 6.112, 0.05 + + # At 20°C es ≈ 23.37 hPa. + assert_in_delta SkewT.saturation_vapor_pressure(20.0), 23.37, 0.2 + + # At -20°C es ≈ 1.254 hPa. + assert_in_delta SkewT.saturation_vapor_pressure(-20.0), 1.254, 0.05 + end + end + + describe "temperature_from_mixing_ratio/2" do + test "round-trips mixing ratios" do + # Pick a pressure and temperature, compute ws, invert, expect the original. + p = 850.0 + + for t <- [-10.0, 0.0, 15.0, 25.0] do + es = SkewT.saturation_vapor_pressure(t) + ws_g_kg = 622.0 * es / (p - es) + recovered = SkewT.temperature_from_mixing_ratio(ws_g_kg, p) + assert_in_delta recovered, t, 0.05 + end + end + end + + describe "dry_adiabat_temperature/2" do + test "θ == T at 1000 mb" do + # At the reference pressure the potential temperature equals the + # actual temperature — so θ = 300 K ⇒ T(1000) = 300 - 273.15 ≈ 26.85°C. + assert_in_delta SkewT.dry_adiabat_temperature(300.0, 1000.0), 26.85, 0.02 + end + + test "a parcel lifted dry-adiabatically cools" do + # Lifting from 1000 mb (θ = 290 K, T ≈ 16.85°C) to 700 mb should + # drop the temperature by ~26°C (roughly 9.8°C/km, ~3 km). + t_surface = SkewT.dry_adiabat_temperature(290.0, 1000.0) + t_700 = SkewT.dry_adiabat_temperature(290.0, 700.0) + assert t_surface - t_700 > 20.0 + end + end + + describe "project/3" do + setup do + chart = SkewT.chart_config(width: 400, height: 500) + %{chart: chart} + end + + test "p_bot lands at the chart bottom and p_top at the top", %{chart: chart} do + {_x_bot, y_bot} = SkewT.project(0.0, chart.p_bot, chart) + {_x_top, y_top} = SkewT.project(0.0, chart.p_top, chart) + + # SVG y grows downward, so p_bot (bottom) has a larger y than p_top. + assert y_bot > y_top + assert_in_delta y_bot, chart.padding_top + chart.plot_height, 0.5 + assert_in_delta y_top, chart.padding_top, 0.5 + end + + test "isotherms skew right as pressure decreases", %{chart: chart} do + # Same temperature at two pressures should have a larger x at + # the lower pressure (top of chart) because of the skew. + {x_bot, _} = SkewT.project(0.0, chart.p_bot, chart) + {x_top, _} = SkewT.project(0.0, chart.p_top, chart) + assert x_top > x_bot + end + + test "temperature range spans the bottom of the chart", %{chart: chart} do + # At p_bot, t_min lands at the left edge of the data area. + {x_min, _} = SkewT.project(chart.t_min, chart.p_bot, chart) + {x_max, _} = SkewT.project(chart.t_max, chart.p_bot, chart) + assert_in_delta x_min, chart.padding_left, 0.5 + assert_in_delta x_max, chart.padding_left + chart.plot_width, 0.5 + end + end + + describe "build/2" do + test "with no profile returns nil" do + assert SkewT.build(nil) == nil + assert SkewT.build([]) == nil + end + + test "returns a map of rendering primitives for a valid profile" do + profile = [ + %{"pres" => 1000.0, "hght" => 100.0, "tmpc" => 25.0, "dwpc" => 18.0}, + %{"pres" => 850.0, "hght" => 1500.0, "tmpc" => 15.0, "dwpc" => 10.0}, + %{"pres" => 700.0, "hght" => 3100.0, "tmpc" => 5.0, "dwpc" => -2.0}, + %{"pres" => 500.0, "hght" => 5600.0, "tmpc" => -15.0, "dwpc" => -25.0} + ] + + chart = SkewT.build(profile) + + assert is_map(chart) + assert chart.width > 0 + assert chart.height > 0 + # Background grid + assert length(chart.isobars) > 0 + assert length(chart.isotherms) > 0 + assert length(chart.dry_adiabats) > 0 + assert length(chart.mixing_ratio_lines) > 0 + # Traces + assert chart.t_points != [] + assert chart.td_points != [] + # The first plotted t_point should correspond to the highest pressure + # (surface-most) profile entry. + assert length(chart.t_points) == 4 + end + + test "filters out profile levels with missing temperature or dewpoint independently" do + profile = [ + %{"pres" => 1000.0, "tmpc" => 25.0, "dwpc" => 18.0}, + %{"pres" => 900.0, "tmpc" => nil, "dwpc" => 15.0}, + %{"pres" => 850.0, "tmpc" => 15.0, "dwpc" => nil}, + %{"pres" => 700.0, "tmpc" => 5.0, "dwpc" => -2.0} + ] + + chart = SkewT.build(profile) + # 3 levels have tmpc, 3 have dwpc — the traces are independent. + assert length(chart.t_points) == 3 + assert length(chart.td_points) == 3 + end + end +end