Graham McIntire
7f1a7fb369
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.
297 lines
9.3 KiB
Elixir
297 lines
9.3 KiB
Elixir
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
|