graham/prop
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Add skew-T log-P diagram to the contact detail page

Browse source 
Render the HRRR pressure-level profile as a full skew-T log-P
diagram on the contact detail page: skewed isotherms, isobars,
dry adiabats, saturation mixing ratio lines, plus the T and Td
traces. Math lives in MicrowavepropWeb.SkewT (Magnus formula,
dry adiabat potential temperature, log-P projection) and is
exercised by a dedicated test module. The atmospheric profile
section now expands by default so the chart is visible without
an extra click.
This commit is contained in:
Graham McIntire 2026-04-14 16:52:15 -05:00
parent 28b6ca0ac4
commit 7f1a7fb369
No known key found for this signature in database
GPG key ID: F4ABF488E6029E59
3 changed files with 601 additions and 2 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

178
lib/microwaveprop_web/live/contact_live/show.ex
View file

@ -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
</div>
</div>
<%= if profile.profile && profile.profile != [] do %>
<div class="overflow-x-auto">
<.skew_t_chart profile={profile.profile} />
<div class="overflow-x-auto mt-4">
<table class="table table-xs table-zebra">
<thead>
<tr>
@ -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"""
<div :if={@chart && @usable_points >= 2} class="mt-3">
<div class="text-xs font-semibold uppercase tracking-wider opacity-60 mb-1">
Skew-T log-P
</div>
<svg
viewBox={"0 0 #{@chart.width} #{@chart.height}"}
class="w-full max-w-[600px] bg-base-100 rounded border border-base-300"
role="img"
aria-label="Skew-T log-P diagram"
>
<%!-- Plot background --%>
<rect
x={@chart.padding_left}
y={@chart.padding_top}
width={@chart.plot_width}
height={@chart.plot_height}
fill="currentColor"
class="text-base-200/40"
/>
<%!-- Clip background curves to the plot area so skewed lines
don't bleed over the axes. --%>
<defs>
<clipPath id="skewt-clip">
<rect
x={@chart.padding_left}
y={@chart.padding_top}
width={@chart.plot_width}
height={@chart.plot_height}
/>
</clipPath>
</defs>
<g clip-path="url(#skewt-clip)">
<%!-- Dry adiabats (orange curves) --%>
<path
:for={adiabat <- @chart.dry_adiabats}
d={adiabat.path}
fill="none"
stroke="#f59e0b"
stroke-width="0.6"
stroke-opacity="0.45"
/>
<%!-- Mixing ratio lines (dashed teal) --%>
<path
:for={line <- @chart.mixing_ratio_lines}
d={line.path}
fill="none"
stroke="#0d9488"
stroke-width="0.6"
stroke-opacity="0.5"
stroke-dasharray="2 3"
/>
<%!-- Isotherms (gray) --%>
<line
:for={iso <- @chart.isotherms}
x1={iso.x1}
y1={iso.y1}
x2={iso.x2}
y2={iso.y2}
stroke="currentColor"
stroke-width="0.5"
class="text-base-content/25"
/>
<%!-- Isobars (gray, across plot) --%>
<line
:for={bar <- @chart.isobars}
x1={bar.x_start}
y1={bar.y}
x2={bar.x_end}
y2={bar.y}
stroke="currentColor"
stroke-width="0.5"
class="text-base-content/25"
/>
<%!-- Dewpoint trace --%>
<path
d={@chart.td_path}
fill="none"
stroke={@chart.dewpoint_color}
stroke-width="2"
stroke-linecap="round"
stroke-linejoin="round"
/>
<circle
:for={{x, y} <- @chart.td_points}
cx={x}
cy={y}
r="2"
fill={@chart.dewpoint_color}
/>
<%!-- Temperature trace --%>
<path
d={@chart.t_path}
fill="none"
stroke={@chart.temp_color}
stroke-width="2"
stroke-linecap="round"
stroke-linejoin="round"
/>
<circle
:for={{x, y} <- @chart.t_points}
cx={x}
cy={y}
r="2"
fill={@chart.temp_color}
/>
</g>
<%!-- Pressure axis labels (left side) --%>
<text
:for={bar <- @chart.isobars}
x={@chart.padding_left - 6}
y={bar.label_y}
font-size="10"
text-anchor="end"
fill="currentColor"
class="text-base-content/60"
>
{bar.label}
</text>
<%!-- Temperature axis labels (bottom, only those inside the visible band) --%>
<text
:for={iso <- @chart.isotherms}
:if={iso.label_x && iso.t_c >= -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
</text>
<%!-- Legend --%>
<g transform={"translate(#{@chart.padding_left + 8}, #{@chart.padding_top + 12})"}>
<line x1="0" y1="0" x2="18" y2="0" stroke={@chart.temp_color} stroke-width="2" />
<text x="22" y="3" font-size="10" fill="currentColor" class="text-base-content/70">
Temperature
</text>
<line x1="0" y1="14" x2="18" y2="14" stroke={@chart.dewpoint_color} stroke-width="2" />
<text x="22" y="17" font-size="10" fill="currentColor" class="text-base-content/70">
Dewpoint
</text>
</g>
</svg>
<div class="text-[11px] opacity-60 mt-1">
Orange curves: dry adiabats · dashed teal: saturation mixing ratio ·
thin gray: isotherms and isobars
</div>
</div>
"""
end
end

297
lib/microwaveprop_web/skew_t.ex Normal file
View file

@ -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

128
test/microwaveprop_web/skew_t_test.exs Normal file
View file

@ -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