prop/lib/microwaveprop_web/live/skewt_svg.ex
Graham McIntire 77bc894f56
feat(skewt): add wet-bulb / virtual-T traces and right-edge level rail
Three additions that bring the diagram closer to the SPC sounding
viewer's full-information chart:

  * Wet-bulb temperature trace, plotted as a thinner blue line under
    the dewpoint trace. Comes from a Stipanuk-style bracket between
    the level's dewpoint and dry-bulb temperatures (already used by
    `SkewtParams` for the wet-bulb-zero level — now public so the
    renderer can call it per level).

  * Environment virtual-temperature trace and parcel virtual-T
    trajectory, both as thin pinkish dashed lines. The parcel
    virtual-T is the actual buoyancy curve CAPE integrates against,
    so showing it next to the dry-bulb parcel ascent makes the
    moisture correction visible. New `SkewtParams.virtual_temp_c/2`
    and `/3` helpers, plus `SkewtParams.mixing_ratio/2`.

  * Critical-level rail on the right edge. Each of LCL / LFC / EL /
    0 °C / WBZ that resolves on the current profile gets a short
    horizontal tick at the right of the plot and a `<label> <p> mb`
    tag just outside it — matches the SPC chart's `143 mb (mix)` /
    `158 mb (mix)` annotation column. SkewtLive computes the
    pressures (and converts the height-keyed 0 °C and WBZ from
    SkewtParams back to pressure via `pressure_at_height/2`) before
    handing them to `SkewtSvg.render/2`.

Wind-derived overlays (hodograph, wind barbs) remain off — HRRR
profile data is wind-free above 10 m AGL.

mix test: 27 / 27 across the four /skewt-related test files,
2,902 / 2,902 + 221 properties on the full suite. mix credo --strict
clean (one cyclomatic-complexity refactor extracted along the way).
2026-04-25 14:24:28 -05:00

602 lines
20 KiB
Elixir
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

defmodule MicrowavepropWeb.SkewtSvg do
@moduledoc """
Render an HRRR vertical profile as a Skew-T-Log-P diagram in SVG.
HRRR's persisted profile stores `pres`, `hght`, `tmpc`, `dwpc` per
pressure level — same canonical shape `Microwaveprop.Weather.SoundingParams`
consumes. This renderer ignores anything else.
Wind barbs are deliberately omitted: HRRR persists wind only at 10 m
AGL, not on the pressure-level profile, so a per-level barb column
would be misleading.
"""
alias Microwaveprop.Weather.SkewtParams
# Plot canvas (SVG user units).
@width 720
@height 640
@left 60
@right 660
@top 30
@bottom 580
@p_bottom 1050.0
@p_top 100.0
@t_min -40.0
@t_max 40.0
# Skew factor — pixels of horizontal offset added per pixel of
# height. 1.0 puts isotherms at 45°.
@skew 1.0
@doc """
Render a Skew-T diagram for `profile` (a list of `%{"pres", "hght",
"tmpc", "dwpc"}` maps) and return an inline SVG string suitable for
`Phoenix.HTML.raw/1`.
When `:parcel_trace` is supplied (typically from
`Microwaveprop.Weather.SkewtParams.derive/1`) a dashed grey parcel-
ascent line is overlaid on the temperature trace, matching the
surface-parcel rendition the SPC sounding viewer publishes.
"""
@spec render([map()], keyword()) :: iodata()
def render(profile, opts \\ []) when is_list(profile) do
cleaned = clean_profile(profile)
parcel = Keyword.get(opts, :parcel_trace, [])
critical = Keyword.get(opts, :critical_levels, [])
plot_w = @right - @left
plot_h = @bottom - @top
[
~s|<svg viewBox="0 0 #{@width} #{@height}" xmlns="http://www.w3.org/2000/svg" class="w-full h-auto" role="img" aria-label="Skew-T-Log-P diagram">|,
style_block(),
# Define a rectangular clip the plot interior shares so any
# grid/trace polyline whose vertices stretch past the frame
# (e.g. skewed isotherms exiting the right edge near the top)
# is visually cropped. The frame and axis labels live *outside*
# this group, so the box outline still reads clean.
~s|<clipPath id="skewt-plot-clip"><rect x="#{@left}" y="#{@top}" width="#{plot_w}" height="#{plot_h}"/></clipPath>|,
~s|<g clip-path="url(#skewt-plot-clip)">|,
isobars(),
isotherms(),
dry_adiabats(),
moist_adiabats(),
mixing_ratio_lines(),
virtual_temp_trace(cleaned),
wet_bulb_trace(cleaned),
profile_traces(cleaned),
parcel_trace_line(parcel),
parcel_virtual_trace(parcel),
"</g>",
frame(),
axis_labels(),
critical_level_markers(critical),
cursor_layer(),
"</svg>"
]
end
# ── Wet-bulb / virtual-T traces (NOAA SPC supplementary lines) ─────
defp wet_bulb_trace([]), do: ""
defp wet_bulb_trace(profile) do
pts =
profile
|> Enum.sort_by(& &1.pres, :desc)
|> Enum.flat_map(fn lvl ->
case SkewtParams.wet_bulb_c(%{
tmpc: lvl.tmpc,
dwpc: lvl.dwpc,
pres: lvl.pres
}) do
nil -> []
t_w -> [{temperature_x(t_w, pressure_y(lvl.pres)), pressure_y(lvl.pres)}]
end
end)
|> clip_polyline()
polyline_or_empty(pts, "wet-bulb-line")
end
defp virtual_temp_trace([]), do: ""
defp virtual_temp_trace(profile) do
pts =
profile
|> Enum.sort_by(& &1.pres, :desc)
|> Enum.map(fn lvl ->
t_v = SkewtParams.virtual_temp_c(lvl.tmpc, lvl.dwpc, lvl.pres)
y = pressure_y(lvl.pres)
{temperature_x(t_v, y), y}
end)
|> clip_polyline()
polyline_or_empty(pts, "virtual-temp-line")
end
# Parcel ascent uses the saturated mixing ratio at its own
# temperature above the LCL — that's the virtual-T trajectory used
# in the actual CAPE integration. SPC draws this as a thin red
# dashed line slightly to the right of the dry-bulb parcel trace.
defp parcel_virtual_trace([]), do: ""
defp parcel_virtual_trace(trace) do
pts =
trace
|> Enum.sort_by(& &1.pres, :desc)
|> Enum.map(fn %{pres: p, tmpc: t} ->
t_v = SkewtParams.virtual_temp_c(t, p)
y = pressure_y(p)
{temperature_x(t_v, y), y}
end)
|> clip_polyline()
polyline_or_empty(pts, "parcel-virtual-line")
end
# Critical-level annotations on the right edge — each entry is a
# `%{label, pressure_mb}` map; the label appears outside the frame
# at the right so it doesn't obscure the trace, with a short
# in-plot tick. Mirrors SPC's `143 mb`, `158 mb (mix)`, etc. rail.
defp critical_level_markers([]), do: ""
defp critical_level_markers(levels) do
Enum.flat_map(levels, fn
%{label: label, pressure_mb: p} when is_number(p) ->
y = pressure_y(p)
if y >= @top and y <= @bottom do
[render_critical_marker(label, p, y)]
else
[]
end
_ ->
[]
end)
end
defp render_critical_marker(label, pressure_mb, y) do
[
~s|<line class="critical-tick" x1="#{@right - 8}" y1="#{r(y)}" x2="#{@right}" y2="#{r(y)}"/>|,
~s|<text class="critical-label" x="#{@right + 4}" y="#{r(y) + 3}">#{label} #{round(pressure_mb)} mb</text>|
]
end
# Surface-parcel ascent line (dashed grey), matching the SPC viewer's
# convention. The trace is a list of `%{pres, tmpc}` from surface up
# through every level above the LCL.
defp parcel_trace_line([]), do: ""
defp parcel_trace_line(trace) do
pts =
trace
|> Enum.sort_by(& &1.pres, :desc)
|> Enum.map(fn %{pres: p, tmpc: t} ->
y = pressure_y(p)
{temperature_x(t, y), y}
end)
|> clip_polyline()
case pts do
[] -> ""
_ -> ~s|<polyline class="parcel-line" points="#{points_str(pts)}" />|
end
end
defp points_str(points) do
Enum.map_join(points, " ", fn {x, y} -> "#{r(x)},#{r(y)}" end)
end
# Crosshair layer the JS hook (`assets/js/skewt_hook.js`) drives.
# Initially `display: none`; the hook flips display + updates element
# coordinates on `mousemove`.
defp cursor_layer do
"""
<g class="skewt-cursor" data-cursor style="display:none">
<line class="skewt-cursor-line" x1="#{@left}" x2="#{@right}" y1="0" y2="0"/>
<circle class="skewt-cursor-dot-t" data-cursor-dot-t cx="0" cy="0" r="3"/>
<circle class="skewt-cursor-dot-d" data-cursor-dot-d cx="0" cy="0" r="3"/>
<rect class="skewt-cursor-readout-bg" data-cursor-readout-bg x="0" y="0" width="148" height="62" rx="3"/>
<text class="skewt-cursor-readout" data-cursor-readout x="0" y="0">
<tspan data-cursor-line-1 x="0" dy="14"></tspan>
<tspan data-cursor-line-2 x="0" dy="14"></tspan>
<tspan data-cursor-line-3 x="0" dy="14"></tspan>
<tspan data-cursor-line-4 x="0" dy="14"></tspan>
</text>
</g>
"""
end
# NOAA SPC Skew-T-Log-P colour set, copied from the SPC sounding viewer
# (https://www.spc.noaa.gov/exper/soundings/) so the chart reads the
# same way as the public reference product:
#
# * Temperature trace .................. solid red (#ff0000)
# * Dewpoint trace ..................... solid green (#00aa00)
# * Isotherms .......................... dashed teal (#22d3a4-ish)
# * Dry adiabats ....................... solid orange (#f59e0b)
# * Moist (saturated) adiabats ......... dashed sky-blue (#38bdf8)
# * Mixing-ratio lines ................. dotted darker green (#16a34a),
# values 0.420 g/kg, drawn only
# in the lower atmosphere
# * Isobars ............................ light grey horizontal lines
defp style_block do
"""
<defs>
<style>
.frame { fill: none; stroke: currentColor; stroke-width: 1.2; }
.grid-major { fill: none; stroke: currentColor; stroke-opacity: 0.45; stroke-width: 0.6; }
.grid-minor { fill: none; stroke: currentColor; stroke-opacity: 0.18; stroke-width: 0.4; }
.grid-iso { fill: none; stroke: #14b8a6; stroke-opacity: 0.55; stroke-width: 0.5; stroke-dasharray: 4 3; }
.grid-dry { fill: none; stroke: #f59e0b; stroke-opacity: 0.65; stroke-width: 0.5; }
.grid-moist { fill: none; stroke: #38bdf8; stroke-opacity: 0.55; stroke-width: 0.5; stroke-dasharray: 4 3; }
.grid-mix { fill: none; stroke: #16a34a; stroke-opacity: 0.55; stroke-width: 0.5; stroke-dasharray: 1 3; }
.axis { font: 10px ui-sans-serif, system-ui, sans-serif; fill: currentColor; }
.iso-label { font: 9px ui-sans-serif, system-ui, sans-serif; fill: #14b8a6; fill-opacity: 0.85; }
.mix-label { font: 9px ui-sans-serif, system-ui, sans-serif; fill: #16a34a; fill-opacity: 0.85; }
.pres-label { font: 9px ui-sans-serif, system-ui, sans-serif; fill: currentColor; fill-opacity: 0.7; }
.temp-line { fill: none; stroke: #ff0000; stroke-width: 2.5; stroke-linejoin: round; stroke-linecap: round; }
.dew-line { fill: none; stroke: #00aa00; stroke-width: 2.5; stroke-linejoin: round; stroke-linecap: round; }
.parcel-line { fill: none; stroke: #475569; stroke-width: 1.5; stroke-dasharray: 5 4; stroke-linejoin: round; stroke-linecap: round; }
.parcel-virtual-line { fill: none; stroke: #f87171; stroke-width: 1; stroke-dasharray: 2 3; stroke-linejoin: round; stroke-linecap: round; }
.virtual-temp-line { fill: none; stroke: #fb7185; stroke-width: 1; stroke-dasharray: 1 3; stroke-linejoin: round; stroke-linecap: round; }
.wet-bulb-line { fill: none; stroke: #2563eb; stroke-width: 1.4; stroke-linejoin: round; stroke-linecap: round; opacity: 0.85; }
.critical-tick { stroke: currentColor; stroke-width: 1.2; }
.critical-label { font: 10px ui-sans-serif, system-ui, sans-serif; fill: currentColor; fill-opacity: 0.85; }
.skewt-cursor-line { stroke: currentColor; stroke-opacity: 0.5; stroke-width: 0.8; stroke-dasharray: 3 3; pointer-events: none; }
.skewt-cursor-dot-t { fill: #ff0000; pointer-events: none; }
.skewt-cursor-dot-d { fill: #00aa00; pointer-events: none; }
.skewt-cursor-readout-bg { fill: rgba(15,23,42,0.85); pointer-events: none; }
.skewt-cursor-readout { font: 11px ui-monospace, SFMono-Regular, monospace; fill: #f8fafc; pointer-events: none; }
</style>
</defs>
"""
end
# ── Geometry exposed to the hover hook ──────────────────────────────
@doc """
Plot geometry constants the JS hover hook needs to invert pixel-Y
back to pressure and to project temperature values onto the skewed
diagram. Returned as a plain map so it serializes to JSON via
`Jason.encode!/1` straight onto a `data-` attribute.
"""
@spec geometry() :: map()
def geometry do
%{
width: @width,
height: @height,
left: @left,
right: @right,
top: @top,
bottom: @bottom,
p_bottom: @p_bottom,
p_top: @p_top,
t_min: @t_min,
t_max: @t_max,
skew: @skew
}
end
# ── Coordinate transforms ───────────────────────────────────────────
@doc false
def pressure_y(p) when is_number(p) do
# High pressure → bottom of plot, low pressure → top.
@bottom -
(@bottom - @top) * (:math.log(@p_bottom) - :math.log(p)) /
(:math.log(@p_bottom) - :math.log(@p_top))
end
@doc false
def temperature_x(t, y) when is_number(t) and is_number(y) do
plot_w = @right - @left
base_x = @left + (t - @t_min) / (@t_max - @t_min) * plot_w
base_x + @skew * (@bottom - y)
end
defp clamp_x(x), do: x |> max(@left) |> min(@right)
defp clip_segment({x1, y1}, {x2, y2}) do
cond do
x1 < @left and x2 < @left -> nil
x1 > @right and x2 > @right -> nil
true -> {{clamp_x(x1), y1}, {clamp_x(x2), y2}}
end
end
# ── Frame & axes ────────────────────────────────────────────────────
defp frame do
~s|<rect class="frame" x="#{@left}" y="#{@top}" width="#{@right - @left}" height="#{@bottom - @top}" />|
end
defp axis_labels do
pressures = [1000, 850, 700, 500, 400, 300, 200, 100]
pres_labels =
for p <- pressures do
y = pressure_y(p)
~s|<text class="pres-label" x="#{@left - 6}" y="#{y + 3}" text-anchor="end">#{p}</text>|
end
temps = -40..40//10
temp_labels =
for t <- temps do
x = temperature_x(t, @bottom)
if x >= @left and x <= @right do
~s|<text class="axis" x="#{x}" y="#{@bottom + 14}" text-anchor="middle">#{t}</text>|
else
""
end
end
[
pres_labels,
temp_labels,
~s|<text class="axis" x="#{@left - 38}" y="#{(@top + @bottom) / 2}" text-anchor="middle" transform="rotate(-90 #{@left - 38} #{(@top + @bottom) / 2})">Pressure (mb)</text>|,
~s|<text class="axis" x="#{(@left + @right) / 2}" y="#{@bottom + 30}" text-anchor="middle">Temperature (°C)</text>|
]
end
# ── Background grid ─────────────────────────────────────────────────
defp isobars do
pressures = [1000, 925, 850, 700, 500, 400, 300, 250, 200, 150, 100]
for p <- pressures do
y = pressure_y(p)
cls = if p in [1000, 850, 700, 500, 300, 200, 100], do: "grid-major", else: "grid-minor"
~s|<line class="#{cls}" x1="#{@left}" y1="#{y}" x2="#{@right}" y2="#{y}" />|
end
end
defp isotherms do
# Skewed isotherms every 10°C across the plot. Each line goes from
# (T at bottom) to (T at top), with skew already baked into
# temperature_x. Clip to the plot box.
for t <- -100..40//10 do
render_isotherm(t)
end
end
defp render_isotherm(t) do
p1 = {temperature_x(t, @bottom), @bottom}
p2 = {temperature_x(t, @top), @top}
case clip_segment(p1, p2) do
nil -> ""
clipped -> isotherm_with_label(t, clipped)
end
end
defp isotherm_with_label(t, {{x1, y1}, {x2, y2}}) do
line = ~s|<line class="grid-iso" x1="#{r(x1)}" y1="#{r(y1)}" x2="#{r(x2)}" y2="#{r(y2)}" />|
[line, isotherm_label(t, x1)]
end
defp isotherm_label(t, x1) when t in [-40, -20, 0, 20, 40] do
if x1 > @left + 2 and x1 < @right - 12 do
~s|<text class="iso-label" x="#{r(x1) + 2}" y="#{@bottom - 2}">#{t}°</text>|
else
""
end
end
defp isotherm_label(_t, _x1), do: ""
defp dry_adiabats do
# Potential-temperature curves: T_K(p) = θ * (p / 1000)^(R/cp).
# Sample every theta = -30..150°C every 10° at 50-mb pressure steps.
for theta_c <- -30..150//10 do
theta_k = theta_c + 273.15
pts =
sample_pressures()
|> Enum.map(fn p ->
t_k = theta_k * :math.pow(p / 1000.0, 0.2854)
t_c = t_k - 273.15
y = pressure_y(p)
x = temperature_x(t_c, y)
{x, y}
end)
|> clip_polyline()
polyline_or_empty(pts, "grid-dry")
end
end
defp moist_adiabats do
# Saturated pseudo-adiabats — quick approximation by integrating
# dT/dp via the moist-adiabatic lapse rate. Good enough to put
# green dashes in the right place; not used for parcel ascent.
for theta_e_c <- 0..40//4 do
pts =
(theta_e_c + 273.15)
|> moist_adiabat_curve()
|> Enum.map(fn {p, t_c} ->
y = pressure_y(p)
x = temperature_x(t_c, y)
{x, y}
end)
|> clip_polyline()
polyline_or_empty(pts, "grid-moist")
end
end
defp moist_adiabat_curve(theta_e_k) do
# Integrate from 1000 mb upward using dT/dz = Γ_m and dp/dz from
# hydrostatic. Cheap approximation good enough for plotting.
starting_t = theta_e_k - 273.15
sample_pressures()
|> Enum.reduce({[], starting_t, 1000.0}, fn p, {acc, t, prev_p} ->
dp = p - prev_p
# Δz from hydrostatic with mean T_K ≈ t + 273
mean_t_k = max(t + 273.15, 200.0)
dz = -dp * 287.0 * mean_t_k / (9.81 * p)
gamma_m = moist_lapse_c_per_m(t, p)
new_t = t + dz * gamma_m
{[{p, new_t} | acc], new_t, p}
end)
|> elem(0)
|> Enum.reverse()
end
defp moist_lapse_c_per_m(t_c, p_mb) do
t_k = t_c + 273.15
# Saturation mixing ratio (g/g) via Buck's eq.
es = 6.1121 * :math.exp((18.678 - t_c / 234.5) * (t_c / (257.14 + t_c)))
qs = 0.622 * es / (p_mb - es)
lv = 2.5e6
cp = 1004.0
rd = 287.0
rv = 461.0
g = 9.81
num = g * (1.0 + lv * qs / (rd * t_k))
den = cp + lv * lv * qs / (rv * t_k * t_k)
-num / den
end
# NOAA SPC mixing-ratio set, in g/kg. Same values quoted along the
# bottom of the SPC sounding charts (0.4 / 0.7 / 1 / 2 / 3 / 5 / 8 /
# 12 / 16 / 20). Lines drawn from the surface up to 600 mb only —
# above that the saturation values get unhelpfully close together.
@mixing_ratios [0.4, 0.7, 1, 2, 3, 5, 8, 12, 16, 20]
@mix_top_pres 600.0
@mix_label_pres 700.0
defp mixing_ratio_lines do
Enum.map(@mixing_ratios, fn w ->
pts =
sample_pressures()
|> Enum.filter(&(&1 >= @mix_top_pres))
|> Enum.map(fn p ->
# Solve for T from w = 0.622 * e_s(T) / (p - e_s).
e = w / 1000.0 * p / (0.622 + w / 1000.0)
t_d = dewpoint_from_vapor(e)
y = pressure_y(p)
x = temperature_x(t_d, y)
{x, y}
end)
|> clip_polyline()
label = mix_label(w)
[polyline_or_empty(pts, "grid-mix"), label]
end)
end
# Anchor each mixing-ratio line with its g/kg value at ~700 mb,
# mirroring SPC's row of green numbers along the diagram's lower band.
defp mix_label(w) do
e = w / 1000.0 * @mix_label_pres / (0.622 + w / 1000.0)
t_d = dewpoint_from_vapor(e)
y = pressure_y(@mix_label_pres)
x = temperature_x(t_d, y)
if x >= @left and x <= @right do
label =
if w >= 1 do
w |> trunc() |> Integer.to_string()
else
:erlang.float_to_binary(w * 1.0, decimals: 1)
end
~s|<text class="mix-label" x="#{r(x)}" y="#{r(y - 2)}" text-anchor="middle">#{label}</text>|
else
""
end
end
# Inverse Buck — find T given vapor pressure e (hPa).
defp dewpoint_from_vapor(e) when e <= 0, do: -80.0
defp dewpoint_from_vapor(e) do
# Magnus form: T = 243.5 * ln(e/6.112) / (17.67 - ln(e/6.112))
ratio = :math.log(e / 6.112)
243.5 * ratio / (17.67 - ratio)
end
defp sample_pressures do
# 1050 down to 100 in 25-mb steps.
100..1050//25 |> Enum.to_list() |> Enum.reverse()
end
defp clip_polyline(points) do
Enum.filter(points, fn {x, y} -> x >= @left - 200 and x <= @right + 200 and y >= @top and y <= @bottom end)
end
defp polyline_or_empty([], _cls), do: ""
defp polyline_or_empty(points, cls) do
pts =
Enum.map_join(points, " ", fn {x, y} -> "#{r(x)},#{r(y)}" end)
~s|<polyline class="#{cls}" points="#{pts}" />|
end
# ── Profile traces ──────────────────────────────────────────────────
defp profile_traces([]), do: ""
defp profile_traces(profile) do
sorted = Enum.sort_by(profile, & &1.pres, :desc)
temp_pts =
sorted
|> Enum.map(fn p ->
y = pressure_y(p.pres)
x = temperature_x(p.tmpc, y)
{x, y}
end)
|> Enum.filter(fn {_x, y} -> y >= @top and y <= @bottom end)
dew_pts =
sorted
|> Enum.filter(& &1.dwpc)
|> Enum.map(fn p ->
y = pressure_y(p.pres)
x = temperature_x(p.dwpc, y)
{x, y}
end)
|> Enum.filter(fn {_x, y} -> y >= @top and y <= @bottom end)
[
polyline_or_empty(dew_pts, "dew-line"),
polyline_or_empty(temp_pts, "temp-line")
]
end
defp clean_profile(profile) do
profile
|> Enum.map(&normalize_level/1)
|> Enum.filter(&(&1.pres && &1.tmpc))
end
defp normalize_level(%{} = level) do
%{
pres: pick_num(level, ["pres", "pres_mb", :pres, :pres_mb]),
hght: pick_num(level, ["hght", "hght_m", :hght, :hght_m]),
tmpc: pick_num(level, ["tmpc", :tmpc]),
dwpc: pick_num(level, ["dwpc", :dwpc])
}
end
defp pick_num(map, keys) do
keys
|> Enum.find_value(fn k -> Map.get(map, k) end)
|> coerce_num()
end
defp coerce_num(n) when is_number(n), do: n * 1.0
defp coerce_num(_), do: nil
defp r(n) when is_float(n), do: Float.round(n, 1)
defp r(n), do: n
end