From f5ec0972cce95122fe85634c61c353f971dc0405 Mon Sep 17 00:00:00 2001 From: Gabriel Wu <13583761+lucifer1004@users.noreply.github.com> Date: Fri, 5 Jun 2026 12:44:53 +0800 Subject: [PATCH 1/2] docs(examples): add rain radar notebook --- CHANGELOG.md | 10 +- docs/make.jl | 31 ++- docs/src/examples/rain_radar.md | 309 +++++++++++++++++++++++++++++ examples/README.md | 16 +- examples/rain_radar.jl | 333 ++++++++++++++++++++++++++++++++ 5 files changed, 688 insertions(+), 11 deletions(-) create mode 100644 docs/src/examples/rain_radar.md create mode 100644 examples/rain_radar.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index fd4ce7c..2fcec2f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,9 +25,13 @@ next release is expected to be `v0.5.0`. `Iterative(EBCM; stable = true)`. - **Performance benchmark suite** and a PR benchmark CI workflow (AirspeedVelocity). -- **`examples/` folder** with five runnable Pluto notebooks — shapes gallery, - angular scattering & polarization, orientation averaging, solver landscape, and - spectral sensitivity — rendered into the documentation site. +- **`examples/` folder** with six runnable Pluto notebooks — shapes gallery, + solver landscape, angular scattering & polarization, orientation averaging, + and spectral sensitivity, plus a rain-radar observables example adapted from + [`@xiongyuup`](https://github.com/xiongyuup)'s + [PR #3](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/pull/3) + with water/ice refractive-index fits and a multi-frequency + brightness-temperature lookup — rendered into the documentation site. ### Changed diff --git a/docs/make.jl b/docs/make.jl index 7b5986b..5a3100d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,9 +11,32 @@ using Documenter # the examples/ environment in its own Pluto worker. const NB_DIR = normpath(joinpath(@__DIR__, "..", "examples")) const NB_OUT = joinpath(@__DIR__, "src", "examples") -const NOTEBOOKS = ["shapes_gallery.jl", "angular_scattering.jl", - "orientation_averaging.jl", "solver_landscape.jl", "spectral_sensitivity.jl"] -const NB_PAGES = ["examples/" * replace(nb, ".jl" => ".md") for nb in NOTEBOOKS] +const NOTEBOOKS = ["shapes_gallery.jl", "solver_landscape.jl", + "angular_scattering.jl", "orientation_averaging.jl", "spectral_sensitivity.jl", + "rain_radar.jl"] + +function notebook_title(nb) + in_markdown = false + for line in eachline(joinpath(NB_DIR, nb)) + s = strip(line) + if !in_markdown + startswith(s, "md\"\"\"") || continue + in_markdown = true + s = strip(s[6:end]) + elseif startswith(s, "\"\"\"") + in_markdown = false + continue + end + startswith(s, "# ") && return String(strip(s[3:end])) + end + String(replace(splitext(nb)[1], "_" => " ")) +end + +const NB_PAGE_PATHS = ["examples/" * replace(nb, ".jl" => ".md") for nb in NOTEBOOKS] +const NB_PAGES = [ + notebook_title(nb) => path + for (nb, path) in zip(NOTEBOOKS, NB_PAGE_PATHS) +] function build_examples() mkpath(NB_OUT) @@ -47,7 +70,7 @@ makedocs(; prettyurls = get(ENV, "CI", "false") == "true", canonical = "https://JuliaRemoteSensing.github.io/TransitionMatrices.jl", edit_link = "main", - size_threshold_ignore = NB_PAGES, # rendered notebooks embed large plots + size_threshold_ignore = NB_PAGE_PATHS, # rendered notebooks embed large plots assets = String[]), pages = [ "Home" => "index.md", diff --git a/docs/src/examples/rain_radar.md b/docs/src/examples/rain_radar.md new file mode 100644 index 0000000..865cea8 --- /dev/null +++ b/docs/src/examples/rain_radar.md @@ -0,0 +1,309 @@ +```@raw html + + + + + +

Rain radar observables

This notebook turns single-particle T-matrices into simple rain radar and brightness-temperature observables. It is adapted from the radiative-transfer prototype contributed by @xiongyuup in PR #3.

The original contribution demonstrates a larger workflow with WRF input, hydrometeor microphysics, radar moments, and microwave radiative transfer. This example keeps only the self-contained core: water refractive index, rain-drop axis ratio, spheroidal T-matrices, a Marshall-Palmer-type drop-size distribution, and integrated dual-polarization moments.

+ +
begin
+    import Pkg
+    Pkg.activate(@__DIR__)
+    Pkg.develop(; path = dirname(@__DIR__))
+    Pkg.instantiate()
+    using TransitionMatrices, Plots
+end
+ + + +

Water, ice, and drop models

The dielectric fits below are compact microwave models for liquid water and ice. Rain drops are represented as oblate spheroids with diameter \(D\) and axis ratio \(c/a\).

+ +
begin
+    c0_m_per_s() = 2.99792458e8
+    c0_mm_per_s() = 1000 * c0_m_per_s()
+
+    function water_refractive_index(freq_hz, T_celsius)
+        λ_m = c0_m_per_s() / freq_hz
+        ε0 = 8.854e-12
+        εs = 78.54 * (1 - 4.579e-3 * (T_celsius - 25) +
+                      1.19e-5 * (T_celsius - 25)^2 -
+                      2.8e-8 * (T_celsius - 25)^3)
+        ε∞ = 5.27137 + 2.16474e-2 * T_celsius - 1.31198e-3 * T_celsius^2
+        α = -16.8129 / (T_celsius + 273) + 6.09265e-2
+        λs = 3.3836e-6 * exp(2513.98 / (T_celsius + 273))
+        σ = 1.1117e-4
+
+        x = (λs / λ_m)^(1 - α)
+        denominator = 1 + 2x * sinpi(α / 2) + (λs / λ_m)^(2 - 2α)
+        εreal = ε∞ + (εs - ε∞) * (1 + x * sinpi(α / 2)) / denominator
+        εimag = (εs - ε∞) * x * cospi(α / 2) / denominator +
+                σ * λ_m / (2π * 2.99792458e8 * ε0)
+
+        sqrt(complex(εreal, εimag))
+    end
+
+    function ice_refractive_index(freq_hz, T_celsius)
+        λ_m = c0_m_per_s() / freq_hz
+        ε0 = 8.854e-12
+        εs = 203.168 + 2.5 * T_celsius + 0.15T_celsius^2
+        ε∞ = 3.168
+        α = 0.288 + 5.2e-3 * T_celsius + 2.3e-4 * T_celsius^2
+        λs = 9.990288e-6 * exp(6643.5 / (T_celsius + 273))
+        σ = 1.1146e-13 * exp(-6291.2 / (T_celsius + 273))
+
+        x = (λs / λ_m)^(1 - α)
+        denominator = 1 + 2x * sinpi(α / 2) + (λs / λ_m)^(2 - 2α)
+        εreal = ε∞ + (εs - ε∞) * (1 + x * sinpi(α / 2)) / denominator
+        εimag = (εs - ε∞) * x * cospi(α / 2) / denominator +
+                σ * λ_m / (2π * c0_m_per_s() * ε0)
+
+        sqrt(complex(εreal, εimag))
+    end
+
+    function rain_axis_ratio(D_mm)
+        D_mm < 0.7 && return 1.0
+        D_mm < 1.5 && return 1.173 - 0.5265D_mm + 0.4698D_mm^2 -
+                            0.1317D_mm^3 - 8.5e-3D_mm^4
+        1.065 - 6.25e-2D_mm - 3.99e-3D_mm^2 +
+            7.66e-4D_mm^3 - 4.095e-5D_mm^4
+    end
+end
+
rain_axis_ratio (generic function with 1 method)
+ +
begin
+    round_complex(z; digits = 4) =
+        complex(round(real(z); digits), round(imag(z); digits))
+
+    material_summary = [
+        (; material = "water", freq_GHz = 9.375, T = 0.0,
+            m = round_complex(water_refractive_index(9.375e9, 0.0))),
+        (; material = "ice", freq_GHz = 9.375, T = -10.0,
+            m = round_complex(ice_refractive_index(9.375e9, -10.0)))
+    ]
+end
+
2-element Vector{@NamedTuple{material::String, freq_GHz::Float64, T::Float64, m::ComplexF64}}:
+ (material = "water", freq_GHz = 9.375, T = 0.0, m = 7.252 + 2.8558im)
+ (material = "ice", freq_GHz = 9.375, T = -10.0, m = 1.7799 + 0.0001im)
+ + +

Single-particle scattering table

Each diameter is solved as a spheroid with the IITM backend. The backward amplitudes feed reflectivity and differential reflectivity; the forward amplitudes feed differential phase; extinction feeds the brightness-temperature toy calculation.

+ +
begin
+    function scattering_table(Ds, freq_hz; temperature_celsius = 0.0,
+            solver = IITM(6, 10, 16))
+        λ = c0_mm_per_s() / freq_hz
+        m = water_refractive_index(freq_hz, temperature_celsius)
+        map(Ds) do D
+            a = D / 2
+            axis_ratio = rain_axis_ratio(D)
+            c = a * axis_ratio
+            shape = Spheroid{Float64, ComplexF64}(a, c, ComplexF64(m))
+            T = transition_matrix(shape, λ, solver)
+            Sback = amplitude_matrix(T, π / 2, 0.0, π / 2, π; λ)
+            Sfwd = amplitude_matrix(T, π / 2, 0.0, π / 2, 0.0; λ)
+            (; D_mm = D,
+                axis_ratio,
+                Shh_back = Sback[2, 2],
+                Svv_back = Sback[1, 1],
+                Shh_forward = Sfwd[2, 2],
+                Svv_forward = Sfwd[1, 1],
+                Cext_mm2 = calc_Cext(T, λ))
+        end
+    end
+
+    spacing(table) = length(table) > 1 ? table[2].D_mm - table[1].D_mm : 1.0
+end
+
spacing (generic function with 1 method)
+ +
begin
+    freq_hz = 9.375e9
+    temperature_celsius = 0.0
+    Ds = collect(0.75:0.75:3.75)
+    solver = IITM(6, 10, 16)
+
+    table = scattering_table(Ds, freq_hz; temperature_celsius, solver)
+    table_summary = [
+        (; D = row.D_mm,
+            axis_ratio = round(row.axis_ratio; digits = 3),
+            Cext = round(row.Cext_mm2; digits = 4))
+        for row in table
+    ]
+end
+
5-element Vector{@NamedTuple{D::Float64, axis_ratio::Float64, Cext::Float64}}:
+ (D = 0.75, axis_ratio = 0.984, Cext = 0.0052)
+ (D = 1.5, axis_ratio = 0.965, Cext = 0.0705)
+ (D = 2.25, axis_ratio = 0.912, Cext = 0.4337)
+ (D = 3.0, axis_ratio = 0.859, Cext = 1.8037)
+ (D = 3.75, axis_ratio = 0.807, Cext = 5.3644)
+ +
table_summary
+
5-element Vector{@NamedTuple{D::Float64, axis_ratio::Float64, Cext::Float64}}:
+ (D = 0.75, axis_ratio = 0.984, Cext = 0.0052)
+ (D = 1.5, axis_ratio = 0.965, Cext = 0.0705)
+ (D = 2.25, axis_ratio = 0.912, Cext = 0.4337)
+ (D = 3.0, axis_ratio = 0.859, Cext = 1.8037)
+ (D = 3.75, axis_ratio = 0.807, Cext = 5.3644)
+ + +

Bulk rain observables

For compactness this uses \(N(D, R) = N_0 \exp[-\Lambda(R)D]\) with \(D\) in millimeters. The constants are chosen to produce a plausible monotonic example, not a complete retrieval algorithm.

+ +
begin
+    drop_distribution(D_mm, rain_rate_mm_h) =
+        rain_rate_mm_h == 0 ? 0.0 : 8.0e3 * exp(-4.1 * rain_rate_mm_h^(-0.21) * D_mm)
+
+    function dualpol_moments(table, freq_hz, rain_rate_mm_h;
+            temperature_celsius = 0.0)
+        λ = c0_mm_per_s() / freq_hz
+        ΔD = spacing(table)
+        m = water_refractive_index(freq_hz, temperature_celsius)
+        Kw2 = abs2((m^2 - 1) / (m^2 + 2))
+        weights = [
+            drop_distribution(row.D_mm, rain_rate_mm_h) * ΔD
+            for row in table
+        ]
+
+        hh = sum(abs2(row.Shh_back) * weights[i] for (i, row) in pairs(table))
+        vv = sum(abs2(row.Svv_back) * weights[i] for (i, row) in pairs(table))
+        Zhh = hh * λ^4 / (π^5 * Kw2)
+        Zvv = vv * λ^4 / (π^5 * Kw2)
+        cov = sum(conj(row.Shh_back) * row.Svv_back * weights[i]
+            for (i, row) in pairs(table))
+        ρhv = abs(cov) / sqrt(hh * vv)
+        kdp = 180 / π * 1e-3 * λ *
+              sum(real(row.Shh_forward - row.Svv_forward) * weights[i]
+                  for (i, row) in pairs(table))
+
+        (; R = rain_rate_mm_h,
+            ZH = round(10log10(Zhh); digits = 2),
+            ZDR = round(10log10(Zhh / Zvv); digits = 2),
+            ρhv = round(ρhv; digits = 4),
+            KDP = round(kdp; digits = 4))
+    end
+
+    function brightness_temperature(table, rain_rate_mm_h;
+            path_km = 5.0, Tatm_K = 270.0)
+        ΔD = spacing(table)
+        α = sum(drop_distribution(row.D_mm, rain_rate_mm_h) *
+                row.Cext_mm2 * ΔD for row in table) / 1e6
+        τ = α * path_km
+        Tatm_K * (1 - exp(-τ))
+    end
+end
+
brightness_temperature (generic function with 1 method)
+ +
begin
+    rain_rates = [1.0, 10.0, 50.0]
+    radar = [
+        dualpol_moments(table, freq_hz, R; temperature_celsius)
+        for R in rain_rates
+    ]
+    tb = [
+        (; R, TB = round(brightness_temperature(table, R); digits = 4))
+        for R in rain_rates
+    ]
+    (; radar, tb)
+end
+
(radar = [(R = 1.0, ZH = 13.23, ZDR = 0.56, ρhv = 0.9987, KDP = 0.0412), (R = 10.0, ZH = 27.33, ZDR = 1.23, ρhv = 0.9963, KDP = 0.5779), (R = 50.0, ZH = 36.03, ZDR = 1.66, ρhv = 0.9962, KDP = 3.1766)], tb = [(R = 1.0, TB = 0.0036), (R = 10.0, TB = 0.0418), (R = 50.0, TB = 0.2254)])
+ + +

Multi-frequency brightness-temperature lookup

The larger PR #3 prototype writes lookup tables to MAT files. Here the same idea stays in memory: build one scattering table per frequency, integrate over the rain distribution, and return a compact table that Pluto renders directly.

+ +
begin
+    function brightness_temperature_lut(Ds, frequencies_hz, rain_rates;
+            temperature_celsius = 0.0, solver = IITM(6, 10, 16),
+            path_km = 5.0, Tatm_K = 270.0)
+        rows = NamedTuple[]
+        for f in frequencies_hz
+            local_table = scattering_table(Ds, f; temperature_celsius, solver)
+            for R in rain_rates
+                push!(rows, (; freq_GHz = round(f / 1e9; digits = 3), R,
+                    TB = round(brightness_temperature(local_table, R;
+                        path_km, Tatm_K); digits = 4)))
+            end
+        end
+        rows
+    end
+
+    lookup_frequencies = [10e9, 18e9]
+    tb_lut = brightness_temperature_lut(Ds, lookup_frequencies, rain_rates;
+        temperature_celsius, solver)
+end
+
6-element Vector{NamedTuple}:
+ (freq_GHz = 10.0, R = 1.0, TB = 0.0042)
+ (freq_GHz = 10.0, R = 10.0, TB = 0.0496)
+ (freq_GHz = 10.0, R = 50.0, TB = 0.2675)
+ (freq_GHz = 18.0, R = 1.0, TB = 0.0167)
+ (freq_GHz = 18.0, R = 10.0, TB = 0.1989)
+ (freq_GHz = 18.0, R = 50.0, TB = 0.9967)
+ +
tb_lut
+
6-element Vector{NamedTuple}:
+ (freq_GHz = 10.0, R = 1.0, TB = 0.0042)
+ (freq_GHz = 10.0, R = 10.0, TB = 0.0496)
+ (freq_GHz = 10.0, R = 50.0, TB = 0.2675)
+ (freq_GHz = 18.0, R = 1.0, TB = 0.0167)
+ (freq_GHz = 18.0, R = 10.0, TB = 0.1989)
+ (freq_GHz = 18.0, R = 50.0, TB = 0.9967)
+ + +

Plots

+ +
let
+    p1 = plot(getfield.(table_summary, :D), getfield.(table_summary, :Cext);
+        marker = :circle, lw = 2, label = "Cext",
+        xlabel = "drop diameter D (mm)", ylabel = "extinction cross section (mm²)")
+    p2 = plot(getfield.(radar, :R), getfield.(radar, :ZH);
+        marker = :circle, lw = 2, xscale = :log10, label = "ZH",
+        xlabel = "rain rate (mm h⁻¹)", ylabel = "reflectivity (dBZ)")
+    p3 = plot(getfield.(radar, :R), getfield.(radar, :ZDR);
+        marker = :circle, lw = 2, xscale = :log10, label = "ZDR",
+        xlabel = "rain rate (mm h⁻¹)", ylabel = "differential reflectivity (dB)")
+    p4 = plot(getfield.(tb, :R), getfield.(tb, :TB);
+        marker = :circle, lw = 2, xscale = :log10, label = "TB",
+        xlabel = "rain rate (mm h⁻¹)", ylabel = "brightness temperature (K)")
+    plot(p1, p2, p3, p4; layout = (2, 2), size = (820, 620))
+end
+ + +
let
+    p = plot(; xscale = :log10, xlabel = "rain rate (mm h⁻¹)",
+        ylabel = "brightness temperature (K)", legend = :topleft)
+    for f in getfield.(tb_lut, :freq_GHz) |> unique
+        subset = filter(row -> row.freq_GHz == f, tb_lut)
+        plot!(p, getfield.(subset, :R), getfield.(subset, :TB);
+            marker = :circle, lw = 2, label = "$(f) GHz")
+    end
+    p
+end
+ + + +

References

+ + +``` diff --git a/examples/README.md b/examples/README.md index 27bebf8..d1664f0 100644 --- a/examples/README.md +++ b/examples/README.md @@ -20,6 +20,10 @@ so it works from a fresh checkout (the first run installs dependencies). - [`shapes_gallery.jl`](shapes_gallery.jl) — **quick-start gallery**: define a spheroid, cylinder, Chebyshev particle, and N-fold prism; compute each T-matrix and its `Qsca`/`Qext`/`g`/`ω`; notes on which solver to use. +- [`solver_landscape.jl`](solver_landscape.jl) — **solver landscape & + convergence**: the two-layer API (`EBCM`, `IITM`, `Iterative`, `stable`) on a + high-aspect spheroid, showing classic EBCM diverging with `nₘₐₓ` while the + stabilized path holds. - [`angular_scattering.jl`](angular_scattering.jl) — **angular scattering & polarization**: the orientation-averaged scattering (Mueller) matrix vs angle (phase function `F₁₁`, linear polarization `-F₁₂/F₁₁`, `F₂₂/F₁₁`, `F₄₄/F₁₁`), @@ -27,10 +31,6 @@ so it works from a fresh checkout (the first run installs dependencies). - [`orientation_averaging.jl`](orientation_averaging.jl) — **orientation averaging**: Mishchenko's analytic `RandomOrientationTransitionMatrix` vs the numerical `orientation_average`, and the numerical average converging to it. -- [`solver_landscape.jl`](solver_landscape.jl) — **solver landscape & - convergence**: the two-layer API (`EBCM`, `IITM`, `Iterative`, `stable`) on a - high-aspect spheroid, showing classic EBCM diverging with `nₘₐₓ` while the - stabilized path holds. - [`spectral_sensitivity.jl`](spectral_sensitivity.jl) — **fast wavelength / refractive-index sensitivities** via the Sh-matrix moment-separation backend. `prepare_sh` runs the geometry quadrature once; differentiating the cheap @@ -38,3 +38,11 @@ so it works from a fresh checkout (the first run installs dependencies). across a whole spectrum at a large speedup over differentiating the from-scratch assembly. Includes plots of the spectra and their sensitivities, and a timing comparison. +- [`rain_radar.jl`](rain_radar.jl) — **rain radar observables**: adapt a + radiative-transfer prototype contributed by + [`@xiongyuup`](https://github.com/xiongyuup) in + [PR #3](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/pull/3) + into a self-contained example. It uses spheroidal IITM T-matrices, a + rain-drop axis-ratio model, water/ice refractive-index fits, and a drop-size + distribution to compute dual-pol radar moments and single- / multi-frequency + brightness-temperature integrals without external WRF/NetCDF/MAT data. diff --git a/examples/rain_radar.jl b/examples/rain_radar.jl new file mode 100644 index 0000000..9509043 --- /dev/null +++ b/examples/rain_radar.jl @@ -0,0 +1,333 @@ +### A Pluto.jl notebook ### +# v1.0.1 + +using Markdown +using InteractiveUtils + +# ╔═╡ 14fa0b92-232b-4a39-b509-88df2453872a +begin + import Pkg + Pkg.activate(@__DIR__) + Pkg.develop(; path = dirname(@__DIR__)) + Pkg.instantiate() + using TransitionMatrices, Plots +end + +# ╔═╡ 65279bf2-9dd8-4a9f-837d-2b15670cc301 +md""" +# Rain radar observables + +This notebook turns single-particle T-matrices into simple rain radar and +brightness-temperature observables. It is adapted from the radiative-transfer +prototype contributed by +[`@xiongyuup`](https://github.com/xiongyuup) in +[`PR #3`](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/pull/3). + +The original contribution demonstrates a larger workflow with WRF input, +hydrometeor microphysics, radar moments, and microwave radiative transfer. This +example keeps only the self-contained core: water refractive index, rain-drop +axis ratio, spheroidal T-matrices, a Marshall-Palmer-type drop-size +distribution, and integrated dual-polarization moments. +""" + +# ╔═╡ 96b07bb0-d7f7-4fdb-9c58-7e43a6a0bf31 +md""" +## Water, ice, and drop models + +The dielectric fits below are compact microwave models for liquid water and +ice. Rain drops are represented as oblate spheroids with diameter ``D`` and +axis ratio ``c/a``. +""" + +# ╔═╡ 109c0c9d-967c-489c-846e-9194a0b6e8a2 +begin + c0_m_per_s() = 2.99792458e8 + c0_mm_per_s() = 1000 * c0_m_per_s() + + function water_refractive_index(freq_hz, T_celsius) + λ_m = c0_m_per_s() / freq_hz + ε0 = 8.854e-12 + εs = 78.54 * (1 - 4.579e-3 * (T_celsius - 25) + + 1.19e-5 * (T_celsius - 25)^2 - + 2.8e-8 * (T_celsius - 25)^3) + ε∞ = 5.27137 + 2.16474e-2 * T_celsius - 1.31198e-3 * T_celsius^2 + α = -16.8129 / (T_celsius + 273) + 6.09265e-2 + λs = 3.3836e-6 * exp(2513.98 / (T_celsius + 273)) + σ = 1.1117e-4 + + x = (λs / λ_m)^(1 - α) + denominator = 1 + 2x * sinpi(α / 2) + (λs / λ_m)^(2 - 2α) + εreal = ε∞ + (εs - ε∞) * (1 + x * sinpi(α / 2)) / denominator + εimag = (εs - ε∞) * x * cospi(α / 2) / denominator + + σ * λ_m / (2π * 2.99792458e8 * ε0) + + sqrt(complex(εreal, εimag)) + end + + function ice_refractive_index(freq_hz, T_celsius) + λ_m = c0_m_per_s() / freq_hz + ε0 = 8.854e-12 + εs = 203.168 + 2.5 * T_celsius + 0.15T_celsius^2 + ε∞ = 3.168 + α = 0.288 + 5.2e-3 * T_celsius + 2.3e-4 * T_celsius^2 + λs = 9.990288e-6 * exp(6643.5 / (T_celsius + 273)) + σ = 1.1146e-13 * exp(-6291.2 / (T_celsius + 273)) + + x = (λs / λ_m)^(1 - α) + denominator = 1 + 2x * sinpi(α / 2) + (λs / λ_m)^(2 - 2α) + εreal = ε∞ + (εs - ε∞) * (1 + x * sinpi(α / 2)) / denominator + εimag = (εs - ε∞) * x * cospi(α / 2) / denominator + + σ * λ_m / (2π * c0_m_per_s() * ε0) + + sqrt(complex(εreal, εimag)) + end + + function rain_axis_ratio(D_mm) + D_mm < 0.7 && return 1.0 + D_mm < 1.5 && return 1.173 - 0.5265D_mm + 0.4698D_mm^2 - + 0.1317D_mm^3 - 8.5e-3D_mm^4 + 1.065 - 6.25e-2D_mm - 3.99e-3D_mm^2 + + 7.66e-4D_mm^3 - 4.095e-5D_mm^4 + end +end + +# ╔═╡ 07e1ca0a-3bf8-45a9-ae4d-8c51971c6f57 +begin + round_complex(z; digits = 4) = + complex(round(real(z); digits), round(imag(z); digits)) + + material_summary = [ + (; material = "water", freq_GHz = 9.375, T = 0.0, + m = round_complex(water_refractive_index(9.375e9, 0.0))), + (; material = "ice", freq_GHz = 9.375, T = -10.0, + m = round_complex(ice_refractive_index(9.375e9, -10.0))) + ] +end + +# ╔═╡ 5d8211a5-3b7a-49a3-8249-474d7a5dc833 +md""" +## Single-particle scattering table + +Each diameter is solved as a spheroid with the IITM backend. The backward +amplitudes feed reflectivity and differential reflectivity; the forward +amplitudes feed differential phase; extinction feeds the brightness-temperature +toy calculation. +""" + +# ╔═╡ 0c68e61c-e0d7-43ea-9c06-f144d042aa9b +begin + function scattering_table(Ds, freq_hz; temperature_celsius = 0.0, + solver = IITM(6, 10, 16)) + λ = c0_mm_per_s() / freq_hz + m = water_refractive_index(freq_hz, temperature_celsius) + map(Ds) do D + a = D / 2 + axis_ratio = rain_axis_ratio(D) + c = a * axis_ratio + shape = Spheroid{Float64, ComplexF64}(a, c, ComplexF64(m)) + T = transition_matrix(shape, λ, solver) + Sback = amplitude_matrix(T, π / 2, 0.0, π / 2, π; λ) + Sfwd = amplitude_matrix(T, π / 2, 0.0, π / 2, 0.0; λ) + (; D_mm = D, + axis_ratio, + Shh_back = Sback[2, 2], + Svv_back = Sback[1, 1], + Shh_forward = Sfwd[2, 2], + Svv_forward = Sfwd[1, 1], + Cext_mm2 = calc_Cext(T, λ)) + end + end + + spacing(table) = length(table) > 1 ? table[2].D_mm - table[1].D_mm : 1.0 +end + +# ╔═╡ 8fce0f52-2fac-4798-84f8-75d084e27d85 +begin + freq_hz = 9.375e9 + temperature_celsius = 0.0 + Ds = collect(0.75:0.75:3.75) + solver = IITM(6, 10, 16) + + table = scattering_table(Ds, freq_hz; temperature_celsius, solver) + table_summary = [ + (; D = row.D_mm, + axis_ratio = round(row.axis_ratio; digits = 3), + Cext = round(row.Cext_mm2; digits = 4)) + for row in table + ] +end + +# ╔═╡ aaef0a83-7657-4c40-a875-bc7f18a7b42e +table_summary + +# ╔═╡ d56cbf67-bad5-45aa-a2b4-2aa4e184fa8d +md""" +## Bulk rain observables + +For compactness this uses +``N(D, R) = N_0 \exp[-\Lambda(R)D]`` with ``D`` in millimeters. The constants are +chosen to produce a plausible monotonic example, not a complete retrieval +algorithm. +""" + +# ╔═╡ d78d1d83-508b-4521-ae44-9068d15f4d14 +begin + drop_distribution(D_mm, rain_rate_mm_h) = + rain_rate_mm_h == 0 ? 0.0 : 8.0e3 * exp(-4.1 * rain_rate_mm_h^(-0.21) * D_mm) + + function dualpol_moments(table, freq_hz, rain_rate_mm_h; + temperature_celsius = 0.0) + λ = c0_mm_per_s() / freq_hz + ΔD = spacing(table) + m = water_refractive_index(freq_hz, temperature_celsius) + Kw2 = abs2((m^2 - 1) / (m^2 + 2)) + weights = [ + drop_distribution(row.D_mm, rain_rate_mm_h) * ΔD + for row in table + ] + + hh = sum(abs2(row.Shh_back) * weights[i] for (i, row) in pairs(table)) + vv = sum(abs2(row.Svv_back) * weights[i] for (i, row) in pairs(table)) + Zhh = hh * λ^4 / (π^5 * Kw2) + Zvv = vv * λ^4 / (π^5 * Kw2) + cov = sum(conj(row.Shh_back) * row.Svv_back * weights[i] + for (i, row) in pairs(table)) + ρhv = abs(cov) / sqrt(hh * vv) + kdp = 180 / π * 1e-3 * λ * + sum(real(row.Shh_forward - row.Svv_forward) * weights[i] + for (i, row) in pairs(table)) + + (; R = rain_rate_mm_h, + ZH = round(10log10(Zhh); digits = 2), + ZDR = round(10log10(Zhh / Zvv); digits = 2), + ρhv = round(ρhv; digits = 4), + KDP = round(kdp; digits = 4)) + end + + function brightness_temperature(table, rain_rate_mm_h; + path_km = 5.0, Tatm_K = 270.0) + ΔD = spacing(table) + α = sum(drop_distribution(row.D_mm, rain_rate_mm_h) * + row.Cext_mm2 * ΔD for row in table) / 1e6 + τ = α * path_km + Tatm_K * (1 - exp(-τ)) + end +end + +# ╔═╡ 383e1078-e216-4810-b83d-cd7a6590749b +begin + rain_rates = [1.0, 10.0, 50.0] + radar = [ + dualpol_moments(table, freq_hz, R; temperature_celsius) + for R in rain_rates + ] + tb = [ + (; R, TB = round(brightness_temperature(table, R); digits = 4)) + for R in rain_rates + ] + (; radar, tb) +end + +# ╔═╡ 38b8df2b-2f47-4f81-b027-f46f5e6e4b75 +md""" +## Multi-frequency brightness-temperature lookup + +The larger PR #3 prototype writes lookup tables to MAT files. Here the same +idea stays in memory: build one scattering table per frequency, integrate over +the rain distribution, and return a compact table that Pluto renders directly. +""" + +# ╔═╡ d6d595e6-4e3f-4c4d-9412-032aef0d940e +begin + function brightness_temperature_lut(Ds, frequencies_hz, rain_rates; + temperature_celsius = 0.0, solver = IITM(6, 10, 16), + path_km = 5.0, Tatm_K = 270.0) + rows = NamedTuple[] + for f in frequencies_hz + local_table = scattering_table(Ds, f; temperature_celsius, solver) + for R in rain_rates + push!(rows, (; freq_GHz = round(f / 1e9; digits = 3), R, + TB = round(brightness_temperature(local_table, R; + path_km, Tatm_K); digits = 4))) + end + end + rows + end + + lookup_frequencies = [10e9, 18e9] + tb_lut = brightness_temperature_lut(Ds, lookup_frequencies, rain_rates; + temperature_celsius, solver) +end + +# ╔═╡ 20ae5289-72e6-4f64-b4e2-b61b6eb62473 +tb_lut + +# ╔═╡ 88f653ee-6d22-4568-8285-57acf2d6c644 +md""" +## Plots +""" + +# ╔═╡ a16f870a-5a48-4638-8491-524afc5f265f +let + p1 = plot(getfield.(table_summary, :D), getfield.(table_summary, :Cext); + marker = :circle, lw = 2, label = "Cext", + xlabel = "drop diameter D (mm)", ylabel = "extinction cross section (mm²)") + p2 = plot(getfield.(radar, :R), getfield.(radar, :ZH); + marker = :circle, lw = 2, xscale = :log10, label = "ZH", + xlabel = "rain rate (mm h⁻¹)", ylabel = "reflectivity (dBZ)") + p3 = plot(getfield.(radar, :R), getfield.(radar, :ZDR); + marker = :circle, lw = 2, xscale = :log10, label = "ZDR", + xlabel = "rain rate (mm h⁻¹)", ylabel = "differential reflectivity (dB)") + p4 = plot(getfield.(tb, :R), getfield.(tb, :TB); + marker = :circle, lw = 2, xscale = :log10, label = "TB", + xlabel = "rain rate (mm h⁻¹)", ylabel = "brightness temperature (K)") + plot(p1, p2, p3, p4; layout = (2, 2), size = (820, 620)) +end + +# ╔═╡ dcf844c1-6d42-47e3-9760-a4d094646796 +let + p = plot(; xscale = :log10, xlabel = "rain rate (mm h⁻¹)", + ylabel = "brightness temperature (K)", legend = :topleft) + for f in getfield.(tb_lut, :freq_GHz) |> unique + subset = filter(row -> row.freq_GHz == f, tb_lut) + plot!(p, getfield.(subset, :R), getfield.(subset, :TB); + marker = :circle, lw = 2, label = "$(f) GHz") + end + p +end + +# ╔═╡ 52c0c938-0a88-41a9-bab8-ed9abf42adee +md""" +## References + +- H. Li, Y. Xiong, and Y. Chen, "Simulation of Complex Meteorological Target + Echoes for Airborne Dual-Polarization Weather Radar Based on Invariant + Imbedding T-Matrix," *IEEE Transactions on Geoscience and Remote Sensing*, + 62, 5105817, 2024. +- G. Zhang, *Weather Radar Polarimetry*, China Meteorological Press, 2018, + pp. 39-40. +- B. R. Brown, M. M. Bell, and A. J. Frambach, "Validation of Simulated + Hurricane Drop Size Distributions Using Polarimetric Radar," *Geophysical + Research Letters*, 42, 2016. +""" + +# ╔═╡ Cell order: +# ╟─65279bf2-9dd8-4a9f-837d-2b15670cc301 +# ╠═14fa0b92-232b-4a39-b509-88df2453872a +# ╟─96b07bb0-d7f7-4fdb-9c58-7e43a6a0bf31 +# ╠═109c0c9d-967c-489c-846e-9194a0b6e8a2 +# ╠═07e1ca0a-3bf8-45a9-ae4d-8c51971c6f57 +# ╟─5d8211a5-3b7a-49a3-8249-474d7a5dc833 +# ╠═0c68e61c-e0d7-43ea-9c06-f144d042aa9b +# ╠═8fce0f52-2fac-4798-84f8-75d084e27d85 +# ╠═aaef0a83-7657-4c40-a875-bc7f18a7b42e +# ╟─d56cbf67-bad5-45aa-a2b4-2aa4e184fa8d +# ╠═d78d1d83-508b-4521-ae44-9068d15f4d14 +# ╠═383e1078-e216-4810-b83d-cd7a6590749b +# ╟─38b8df2b-2f47-4f81-b027-f46f5e6e4b75 +# ╠═d6d595e6-4e3f-4c4d-9412-032aef0d940e +# ╠═20ae5289-72e6-4f64-b4e2-b61b6eb62473 +# ╟─88f653ee-6d22-4568-8285-57acf2d6c644 +# ╠═a16f870a-5a48-4638-8491-524afc5f265f +# ╠═dcf844c1-6d42-47e3-9760-a4d094646796 +# ╟─52c0c938-0a88-41a9-bab8-ed9abf42adee From 95d408eaf9f7c0744107e99e3066ebcf6655f5fb Mon Sep 17 00:00:00 2001 From: Gabriel Wu <13583761+lucifer1004@users.noreply.github.com> Date: Fri, 5 Jun 2026 13:23:40 +0800 Subject: [PATCH 2/2] feat(examples): add radiative transfer workflow Add an optional WRF-like dual-pol radar workflow with synthetic NetCDF data generation, and bump the package plus dependent compat bounds to 0.5. --- .gitignore | 5 + CHANGELOG.md | 11 +- Project.toml | 2 +- README.md | 2 +- docs/src/examples/rain_radar.md | 2 +- docs/src/index.md | 2 +- examples/README.md | 8 + examples/radiative_transfer/Project.toml | 10 + examples/radiative_transfer/README.md | 89 +++++ .../radiative_transfer/generate_fake_wrf.jl | 156 ++++++++ examples/radiative_transfer/shared.jl | 25 ++ .../radiative_transfer/wrf_dualpol_radar.jl | 370 ++++++++++++++++++ examples/rain_radar.jl | 2 +- .../EBCMPrecisionLossEstimators/Project.toml | 2 +- 14 files changed, 677 insertions(+), 9 deletions(-) create mode 100644 examples/radiative_transfer/Project.toml create mode 100644 examples/radiative_transfer/README.md create mode 100644 examples/radiative_transfer/generate_fake_wrf.jl create mode 100644 examples/radiative_transfer/shared.jl create mode 100644 examples/radiative_transfer/wrf_dualpol_radar.jl diff --git a/.gitignore b/.gitignore index 80545b7..49cfb7c 100644 --- a/.gitignore +++ b/.gitignore @@ -21,6 +21,11 @@ Manifest.toml /examples/*.md !/examples/README.md +# Generated data/cache for optional advanced examples. +/examples/radiative_transfer/*.nc +/examples/radiative_transfer/*.mat +/examples/radiative_transfer/*.jls + # Editor and OS metadata .DS_Store .idea/ diff --git a/CHANGELOG.md b/CHANGELOG.md index 2fcec2f..62fcdeb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,8 +6,9 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [Unreleased] -Work on top of `v0.4.0`. Includes a breaking change to solver selection, so the -next release is expected to be `v0.5.0`. +No unreleased changes yet. + +## [0.5.0] - 2026-06-05 ### Added @@ -32,6 +33,9 @@ next release is expected to be `v0.5.0`. [PR #3](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/pull/3) with water/ice refractive-index fits and a multi-frequency brightness-temperature lookup — rendered into the documentation site. +- **Optional radiative-transfer workflow example** with synthetic WRF-like + NetCDF generation and a WRF-style dual-pol radar pipeline under + `examples/radiative_transfer/`. ### Changed @@ -80,7 +84,8 @@ next release is expected to be `v0.5.0`. - Initial release. -[Unreleased]: https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/compare/v0.4.0...HEAD +[Unreleased]: https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/compare/v0.5.0...HEAD +[0.5.0]: https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/compare/v0.4.0...v0.5.0 [0.4.0]: https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/compare/v0.3.1...v0.4.0 [0.3.1]: https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/compare/v0.3.0...v0.3.1 [0.3.0]: https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/compare/v0.2.0...v0.3.0 diff --git a/Project.toml b/Project.toml index 46853cb..bdb4d81 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "TransitionMatrices" uuid = "057c4241-e127-4181-840e-6b4b92e6eef5" -version = "0.4.0" +version = "0.5.0" authors = ["Gabriel Wu and contributors"] [deps] diff --git a/README.md b/README.md index 0cf1ed6..037cf94 100644 --- a/README.md +++ b/README.md @@ -40,4 +40,4 @@ The precision types `Double64`, `Float128`, `ComplexF128`, `Arb`, and `Acb` are re-exported by `TransitionMatrices.jl` and can be directly used after `using TransitionMatrices`. -The `0.4` compatibility line uses `Quadmath.jl` 1.x and `Wigxjpf.jl` 0.3.x. +The `0.5` compatibility line uses `Quadmath.jl` 1.x and `Wigxjpf.jl` 0.3.x. diff --git a/docs/src/examples/rain_radar.md b/docs/src/examples/rain_radar.md index 865cea8..a3508d8 100644 --- a/docs/src/examples/rain_radar.md +++ b/docs/src/examples/rain_radar.md @@ -62,7 +62,7 @@ end denominator = 1 + 2x * sinpi(α / 2) + (λs / λ_m)^(2 - 2α) εreal = ε∞ + (εs - ε∞) * (1 + x * sinpi(α / 2)) / denominator εimag = (εs - ε∞) * x * cospi(α / 2) / denominator + - σ * λ_m / (2π * 2.99792458e8 * ε0) + σ * λ_m / (2π * c0_m_per_s() * ε0) sqrt(complex(εreal, εimag)) end diff --git a/docs/src/index.md b/docs/src/index.md index bafb7a3..b66fbe7 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -41,4 +41,4 @@ The precision types `Double64`, `Float128`, `ComplexF128`, `Arb`, and `Acb` are re-exported by `TransitionMatrices.jl` and can be directly used after `using TransitionMatrices`. -The `0.4` compatibility line uses `Quadmath.jl` 1.x and `Wigxjpf.jl` 0.3.x. +The `0.5` compatibility line uses `Quadmath.jl` 1.x and `Wigxjpf.jl` 0.3.x. diff --git a/examples/README.md b/examples/README.md index d1664f0..03fbbcc 100644 --- a/examples/README.md +++ b/examples/README.md @@ -46,3 +46,11 @@ so it works from a fresh checkout (the first run installs dependencies). rain-drop axis-ratio model, water/ice refractive-index fits, and a drop-size distribution to compute dual-pol radar moments and single- / multi-frequency brightness-temperature integrals without external WRF/NetCDF/MAT data. + +## Advanced Workflows + +- [`radiative_transfer/`](radiative_transfer/) — **WRF-like dual-pol radar + workflow**: an optional standalone environment with synthetic NetCDF data + generation and a WRF-style rain radar pipeline. It keeps heavy NetCDF + workflow dependencies out of the main package and is intended as a smoke-test + and adaptation point for real WRF output. diff --git a/examples/radiative_transfer/Project.toml b/examples/radiative_transfer/Project.toml new file mode 100644 index 0000000..0a7b518 --- /dev/null +++ b/examples/radiative_transfer/Project.toml @@ -0,0 +1,10 @@ +[deps] +NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" +TransitionMatrices = "057c4241-e127-4181-840e-6b4b92e6eef5" + +[compat] +NCDatasets = "0.14" +SpecialFunctions = "2" +TransitionMatrices = "0.5" +julia = "1.10" diff --git a/examples/radiative_transfer/README.md b/examples/radiative_transfer/README.md new file mode 100644 index 0000000..c16ce48 --- /dev/null +++ b/examples/radiative_transfer/README.md @@ -0,0 +1,89 @@ +# Radiative Transfer Workflow + +This optional example adapts the WRF-driven dual-polarization radar prototype +contributed by [`@xiongyuup`](https://github.com/xiongyuup) in +[`PR #3`](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/pull/3). +It is a workflow example, not part of the `TransitionMatrices.jl` library API. + +The workflow reads WRF-like NetCDF fields, precomputes a rain-drop scattering +table with `TransitionMatrices.jl`, integrates a gamma drop-size distribution at +each grid cell, and writes dual-pol radar fields to NetCDF. + +## Setup + +Run from the repository root: + +```sh +julia --project=examples/radiative_transfer -e 'import Pkg; Pkg.develop(path="."); Pkg.instantiate()' +``` + +The dependencies live only in this example environment. They are not added to +the main package. + +## Synthetic Smoke Test + +The generator creates a tiny deterministic WRF-like file. It is only for +checking that I/O and the compute pipeline work; it is not validation data. + +```sh +julia --project=examples/radiative_transfer \ + examples/radiative_transfer/generate_fake_wrf.jl \ + --output examples/radiative_transfer/fake_wrfout.nc +``` + +Then run the radar workflow: + +```sh +JULIA_NUM_THREADS=4 julia --project=examples/radiative_transfer \ + examples/radiative_transfer/wrf_dualpol_radar.jl \ + --input examples/radiative_transfer/fake_wrfout.nc \ + --output examples/radiative_transfer/fake_dualpol.nc \ + --dsd double +``` + +The default diameter grid and IITM settings are intentionally small so the smoke +test finishes quickly. Use the CLI options to increase resolution for real +experiments. + +## Inputs + +The workflow expects these variables: + +- `XLAT`, `XLONG`, `ZNU`, `XTIME` +- `QRAIN`, `QVAPOR`, `PB`, `P`, `T` +- `QNRAIN` for `--dsd double` + +All 4-D microphysics and thermodynamic variables use +`west_east × south_north × bottom_top × Time` ordering. + +## Outputs + +The output NetCDF contains: + +- `ZH`, `ZV`, `ZDR`, `RHOHV`, `KDP` +- copied `QRAIN` and derived/copied `QNRAIN` +- coordinate fields `XLAT`, `XLONG`, `ZNU`, `XTIME` + +Dry or invalid cells are filled with `-999`. + +## Notes + +- `--dsd single` derives the number concentration from a fixed intercept + parameter. +- `--dsd double` uses `QNRAIN` from the input file. +- A serialized scattering table cache is written by default to + `scattering_table.jls`; delete it or change `--cache` when changing + discretization settings beyond the checked metadata. +- Full-bin WRF microphysics is not implemented here. That remains a larger + follow-up because it needs a different input contract and more careful unit + handling. + +## References + +- H. Li, Y. Xiong, and Y. Chen, "Simulation of Complex Meteorological Target + Echoes for Airborne Dual-Polarization Weather Radar Based on Invariant + Imbedding T-Matrix," *IEEE Transactions on Geoscience and Remote Sensing*, + 62, 5105817, 2024. +- B. R. Brown, M. M. Bell, and A. J. Frambach, "Validation of Simulated + Hurricane Drop Size Distributions Using Polarimetric Radar," *Geophysical + Research Letters*, 42, 2016. diff --git a/examples/radiative_transfer/generate_fake_wrf.jl b/examples/radiative_transfer/generate_fake_wrf.jl new file mode 100644 index 0000000..357cb65 --- /dev/null +++ b/examples/radiative_transfer/generate_fake_wrf.jl @@ -0,0 +1,156 @@ +#!/usr/bin/env julia + +using NCDatasets + +include(joinpath(@__DIR__, "shared.jl")) +using .RadiativeTransferShared: MASS_COEFF, RAIN_N0, R_DRY_AIR, WATER_DENSITY, + air_density, write_var + +function usage() + println(""" + Generate a tiny deterministic WRF-like NetCDF file for smoke testing. + + Usage: + julia --project=examples/radiative_transfer \\ + examples/radiative_transfer/generate_fake_wrf.jl [options] + + Options: + --output PATH Output NetCDF path [fake_wrfout.nc] + --nx N west_east grid size [8] + --ny N south_north grid size [6] + --nz N bottom_top grid size [4] + --nt N time size [1] + --help Show this message + """) +end + +function parse_args(args) + opts = Dict( + "output" => "fake_wrfout.nc", + "nx" => "8", + "ny" => "6", + "nz" => "4", + "nt" => "1") + i = 1 + while i <= length(args) + arg = args[i] + if arg == "--help" || arg == "-h" + usage() + exit(0) + elseif haskey(opts, replace(arg, "--" => "")) + i == length(args) && error("missing value for $arg") + opts[replace(arg, "--" => "")] = args[i + 1] + i += 2 + else + error("unknown argument: $arg") + end + end + nx = parse(Int, opts["nx"]) + ny = parse(Int, opts["ny"]) + nz = parse(Int, opts["nz"]) + nt = parse(Int, opts["nt"]) + for (name, value) in (("nx", nx), ("ny", ny), ("nz", nz), ("nt", nt)) + value > 0 || error("invalid --$name: expected a positive integer, got $value") + end + (; output = opts["output"], nx, ny, nz, nt) +end + +function synthetic_fields(nx, ny, nz, nt) + xlat = Array{Float32}(undef, nx, ny, nt) + xlon = similar(xlat) + znu = collect(Float32, range(1.0, 0.15; length = nz)) + xtime = collect(Float32, 0:10:(10 * (nt - 1))) + + shape = (nx, ny, nz, nt) + qrain = zeros(Float32, shape) + qnrain = zeros(Float32, shape) + qvapor = Array{Float32}(undef, shape) + pb = Array{Float32}(undef, shape) + p = Array{Float32}(undef, shape) + theta = Array{Float32}(undef, shape) + refl = fill(Float32(-999), shape) + + cx = (nx + 1) / 2 + cy = (ny + 1) / 2 + sx = max(nx / 4, 1) + sy = max(ny / 4, 1) + + for t in 1:nt + for j in 1:ny, i in 1:nx + xlat[i, j, t] = 34.0f0 + 0.04f0 * Float32(j - 1) + xlon[i, j, t] = -98.0f0 + 0.04f0 * Float32(i - 1) + end + + for k in 1:nz, j in 1:ny, i in 1:nx + height_factor = exp(-0.45 * (k - 1)) + blob = exp(-((i - cx)^2 / (2sx^2) + (j - cy)^2 / (2sy^2) + + (k - 1)^2 / 2.5)) + rain = blob > 0.035 ? 8.0e-4 * blob * height_factor : 0.0 + + qvapor[i, j, k, t] = Float32(0.011 * exp(-0.12 * (k - 1))) + pb[i, j, k, t] = Float32(93000 * exp(-0.10 * (k - 1))) + p[i, j, k, t] = Float32(250 * sinpi((i - 1) / max(nx - 1, 1)) * + cospi((j - 1) / max(ny - 1, 1))) + theta[i, j, k, t] = Float32(2.0 - 1.5 * (k - 1)) + qrain[i, j, k, t] = Float32(rain) + refl[i, j, k, t] = Float32(rain > 0 ? 10 + 45 * sqrt(blob) : -999) + + if rain > 0 + rho = air_density(qvapor[i, j, k, t], pb[i, j, k, t], + p[i, j, k, t], theta[i, j, k, t]) + qrv = rain * rho + lambda = (MASS_COEFF * RAIN_N0 * 6 / qrv)^(1 / 4) + number_density = RAIN_N0 / lambda + qnrain[i, j, k, t] = Float32(number_density / rho) + end + end + end + + (; xlat, xlon, znu, xtime, qrain, qnrain, qvapor, pb, p, theta, refl) +end + +function write_fake_wrf(path, fields) + rm(path; force = true) + ds = NCDataset(path, "c") + try + nx, ny, nz, nt = size(fields.qrain) + defDim(ds, "west_east", nx) + defDim(ds, "south_north", ny) + defDim(ds, "bottom_top", nz) + defDim(ds, "Time", nt) + ds.attrib["title"] = "Synthetic WRF-like smoke-test data" + ds.attrib["source"] = "TransitionMatrices.jl examples/radiative_transfer" + ds.attrib["note"] = "Not validation data; generated only to exercise the pipeline." + + write_var(ds, "XLAT", fields.xlat, ("west_east", "south_north", "Time"); + units = "degree_north") + write_var(ds, "XLONG", fields.xlon, ("west_east", "south_north", "Time"); + units = "degree_east") + write_var(ds, "ZNU", fields.znu, ("bottom_top",); + description = "synthetic eta coordinate") + write_var(ds, "XTIME", fields.xtime, ("Time",); units = "minutes") + + dims4 = ("west_east", "south_north", "bottom_top", "Time") + write_var(ds, "QRAIN", fields.qrain, dims4; units = "kg kg-1") + write_var(ds, "QNRAIN", fields.qnrain, dims4; units = "kg-1") + write_var(ds, "QVAPOR", fields.qvapor, dims4; units = "kg kg-1") + write_var(ds, "PB", fields.pb, dims4; units = "Pa") + write_var(ds, "P", fields.p, dims4; units = "Pa") + write_var(ds, "T", fields.theta, dims4; units = "K", + description = "WRF-style perturbation potential temperature") + write_var(ds, "REFL_10CM", fields.refl, dims4; units = "dBZ", + description = "synthetic placeholder reflectivity") + finally + close(ds) + end +end + +function main(args = ARGS) + opts = parse_args(args) + fields = synthetic_fields(opts.nx, opts.ny, opts.nz, opts.nt) + write_fake_wrf(opts.output, fields) + println("wrote synthetic WRF-like file: ", opts.output) + println("grid: $(opts.nx) x $(opts.ny) x $(opts.nz) x $(opts.nt)") +end + +abspath(PROGRAM_FILE) == abspath(@__FILE__) && main() diff --git a/examples/radiative_transfer/shared.jl b/examples/radiative_transfer/shared.jl new file mode 100644 index 0000000..7f8596d --- /dev/null +++ b/examples/radiative_transfer/shared.jl @@ -0,0 +1,25 @@ +module RadiativeTransferShared + +using NCDatasets: defVar + +export R_DRY_AIR, RAIN_N0, WATER_DENSITY, MASS_COEFF, air_density, write_var + +const R_DRY_AIR = 287.15 +const RAIN_N0 = 8.0e3 +const WATER_DENSITY = 997.0 +const MASS_COEFF = pi * WATER_DENSITY / 6 * 1e-9 + +function air_density(qvapor, pb, p, theta_perturbation) + pressure = pb + p + tempk = (theta_perturbation + 300.0) * (pressure / 100000.0)^(2.0 / 7.0) + 0.622 * pressure / (R_DRY_AIR * tempk * (qvapor + 0.622)) +end + +function write_var(ds, name, data, dims; units = "", description = "") + v = defVar(ds, name, data, dims) + isempty(units) || (v.attrib["units"] = units) + isempty(description) || (v.attrib["description"] = description) + return v +end + +end diff --git a/examples/radiative_transfer/wrf_dualpol_radar.jl b/examples/radiative_transfer/wrf_dualpol_radar.jl new file mode 100644 index 0000000..83729ba --- /dev/null +++ b/examples/radiative_transfer/wrf_dualpol_radar.jl @@ -0,0 +1,370 @@ +#!/usr/bin/env julia + +using NCDatasets +using Serialization +using SpecialFunctions +using TransitionMatrices + +include(joinpath(@__DIR__, "shared.jl")) +using .RadiativeTransferShared: MASS_COEFF, RAIN_N0, R_DRY_AIR, WATER_DENSITY, + air_density, write_var + +const C0_MM_PER_S = 2.99792458e11 +const C0_M_PER_S = 2.99792458e8 +const FILL = Float32(-999) + +function usage() + println(""" + Compute dual-pol rain radar fields from a WRF-like NetCDF file. + + Usage: + julia --project=examples/radiative_transfer \\ + examples/radiative_transfer/wrf_dualpol_radar.jl --input fake_wrfout.nc [options] + + Options: + --input PATH Input WRF-like NetCDF path + --output PATH Output NetCDF path [dualpol_radar.nc] + --dsd MODE DSD mode: single or double [double] + --freq HZ Radar frequency in Hz [9.375e9] + --diameters SPEC Diameter grid start:step:stop in mm [0.75:0.75:3.75] + --cache PATH Serialized scattering-table cache [scattering_table.jls] + --nmax N IITM order [6] + --nr N IITM radial grid [10] + --ntheta N IITM polar grid [16] + --help Show this message + """) +end + +function parse_args(args) + opts = Dict( + "input" => "", + "output" => "dualpol_radar.nc", + "dsd" => "double", + "freq" => "9.375e9", + "diameters" => "0.75:0.75:3.75", + "cache" => "scattering_table.jls", + "nmax" => "6", + "nr" => "10", + "ntheta" => "16") + i = 1 + while i <= length(args) + arg = args[i] + if arg == "--help" || arg == "-h" + usage() + exit(0) + elseif startswith(arg, "--") && haskey(opts, arg[3:end]) + i == length(args) && error("missing value for $arg") + opts[arg[3:end]] = args[i + 1] + i += 2 + else + error("unknown argument: $arg") + end + end + input = opts["input"] + output = opts["output"] + isempty(input) && error("--input is required") + abspath(input) == abspath(output) && + error("--output must differ from --input; both resolve to $(abspath(input))") + dsd = opts["dsd"] + dsd in ("single", "double") || error("--dsd must be single or double") + (; input, + output, + dsd, + freq_hz = parse(Float64, opts["freq"]), + diameters = parse_range(opts["diameters"]), + cache = opts["cache"], + solver = IITM(parse(Int, opts["nmax"]), parse(Int, opts["nr"]), + parse(Int, opts["ntheta"]))) +end + +function parse_range(spec) + parts = parse.(Float64, split(spec, ":")) + bins = if length(parts) == 3 + parts[2] != 0 || error("range step must be nonzero") + collect(parts[1]:parts[2]:parts[3]) + elseif length(parts) == 2 + collect(parts[1]:parts[2]) + else + error("range must be start:step:stop or start:stop") + end + isempty(bins) && error("range must produce at least one bin") + bins +end + +function water_refractive_index(freq_hz, T_celsius) + λ_m = C0_M_PER_S / freq_hz + ε0 = 8.854e-12 + εs = 78.54 * (1 - 4.579e-3 * (T_celsius - 25) + + 1.19e-5 * (T_celsius - 25)^2 - + 2.8e-8 * (T_celsius - 25)^3) + ε∞ = 5.27137 + 2.16474e-2 * T_celsius - 1.31198e-3 * T_celsius^2 + α = -16.8129 / (T_celsius + 273) + 6.09265e-2 + λs = 3.3836e-6 * exp(2513.98 / (T_celsius + 273)) + σ = 1.1117e-4 + + x = (λs / λ_m)^(1 - α) + denominator = 1 + 2x * sinpi(α / 2) + (λs / λ_m)^(2 - 2α) + εreal = ε∞ + (εs - ε∞) * (1 + x * sinpi(α / 2)) / denominator + εimag = (εs - ε∞) * x * cospi(α / 2) / denominator + + σ * λ_m / (2π * C0_M_PER_S * ε0) + + sqrt(complex(εreal, εimag)) +end + +function rain_axis_ratio(D_mm) + D_mm < 0.7 && return 1.0 + D_mm < 1.5 && return 1.173 - 0.5265D_mm + 0.4698D_mm^2 - + 0.1317D_mm^3 - 8.5e-3D_mm^4 + 1.065 - 6.25e-2D_mm - 3.99e-3D_mm^2 + + 7.66e-4D_mm^3 - 4.095e-5D_mm^4 +end + +function scattering_table(Ds, freq_hz; temperature_celsius = 0.0, + solver = IITM(6, 10, 16)) + λ = C0_MM_PER_S / freq_hz + m = water_refractive_index(freq_hz, temperature_celsius) + map(Ds) do D + a = D / 2 + c = a * rain_axis_ratio(D) + shape = Spheroid{Float64, ComplexF64}(a, c, ComplexF64(m)) + T = transition_matrix(shape, λ, solver) + Sback = amplitude_matrix(T, π / 2, 0.0, π / 2, π; λ) + Sfwd = amplitude_matrix(T, π / 2, 0.0, π / 2, 0.0; λ) + (; D_mm = D, + Shh_back = Sback[2, 2], + Svv_back = Sback[1, 1], + Shh_forward = Sfwd[2, 2], + Svv_forward = Sfwd[1, 1]) + end +end + +scattering_solver_params(solver) = (; solver = repr(solver)) + +function scattering_solver_params(solver::IITM) + (; nmax = solver.nₘₐₓ, + nr = solver.Nr, + ntheta = solver.Nϑ, + nphi = solver.Nφ, + rmin = solver.rₘᵢₙ) +end + +function scattering_cache_matches(payload, Ds, freq_hz, solver_params) + names = cache_propertynames(payload) + (:Ds in names && :freq_hz in names && :solver_params in names) || return false + payload.Ds == Ds && payload.freq_hz == freq_hz && + payload.solver_params == solver_params +end + +function cache_propertynames(payload) + try + propertynames(payload) + catch + () + end +end + +function cached_scattering_table(payload, Ds) + :table in cache_propertynames(payload) || return nothing + table = payload.table + isnothing(table) && return nothing + try + length(table) == length(Ds) ? table : nothing + catch + nothing + end +end + +function load_or_compute_scattering_table(cache, Ds, freq_hz; solver) + solver_params = scattering_solver_params(solver) + if !isempty(cache) && isfile(cache) + payload = try + deserialize(cache) + catch err + @warn "ignoring unreadable scattering table cache" cache error = sprint(showerror, err) + nothing + end + if !isnothing(payload) + if scattering_cache_matches(payload, Ds, freq_hz, solver_params) + table = cached_scattering_table(payload, Ds) + if !isnothing(table) + println("loaded scattering table cache: ", cache) + return table + end + @warn "ignoring scattering table cache without a valid table" cache + else + println("ignoring stale scattering table cache: ", cache) + end + end + end + + table = scattering_table(Ds, freq_hz; solver) + if !isempty(cache) + serialize(cache, (; Ds, freq_hz, solver_params, table)) + println("wrote scattering table cache: ", cache) + end + table +end + +struct WRFFields + xlat::Array{Float64, 3} + xlong::Array{Float64, 3} + znu::Vector{Float64} + xtime::Vector{Float64} + qrain::Array{Float64, 4} + qnrain::Union{Nothing, Array{Float64, 4}} + qvapor::Array{Float64, 4} + pb::Array{Float64, 4} + p::Array{Float64, 4} + theta::Array{Float64, 4} +end + +read_array(ds, name) = Float64.(ds[name][ntuple(_ -> Colon(), ndims(ds[name]))...]) + +function read_wrf_fields(path; dsd) + ds = NCDataset(path, "r") + try + qnrain = haskey(ds, "QNRAIN") ? read_array(ds, "QNRAIN") : nothing + dsd == "double" && isnothing(qnrain) && + error("double DSD mode requires QNRAIN in the input file") + WRFFields( + read_array(ds, "XLAT"), + read_array(ds, "XLONG"), + vec(read_array(ds, "ZNU")), + vec(read_array(ds, "XTIME")), + read_array(ds, "QRAIN"), + qnrain, + read_array(ds, "QVAPOR"), + read_array(ds, "PB"), + read_array(ds, "P"), + read_array(ds, "T")) + finally + close(ds) + end +end + +function dsd_parameters(qrv, qn_density, mode) + if mode == "single" + Λ = (MASS_COEFF * RAIN_N0 * gamma(4) / qrv)^(1 / 4) + return (; N0 = RAIN_N0, Λ) + end + + Λ = (MASS_COEFF * gamma(4) * qn_density / qrv)^(1 / 3) + N0 = qn_density * Λ + (; N0, Λ) +end + +function radar_moments(table, freq_hz, N0, Λ) + λ = C0_MM_PER_S / freq_hz + ΔD = length(table) > 1 ? table[2].D_mm - table[1].D_mm : 1.0 + m = water_refractive_index(freq_hz, 0.0) + Kw2 = abs2((m^2 - 1) / (m^2 + 2)) + weights = [N0 * exp(-Λ * row.D_mm) * ΔD for row in table] + + hh = sum(abs2(row.Shh_back) * weights[i] for (i, row) in pairs(table)) + vv = sum(abs2(row.Svv_back) * weights[i] for (i, row) in pairs(table)) + if !(hh > 0 && vv > 0) + return (; ZH = NaN, ZV = NaN, ZDR = NaN, RHOHV = NaN, KDP = NaN) + end + + Zhh = hh * λ^4 / (π^5 * Kw2) + Zvv = vv * λ^4 / (π^5 * Kw2) + cov = sum(conj(row.Shh_back) * row.Svv_back * weights[i] + for (i, row) in pairs(table)) + ρhv = abs(cov) / sqrt(hh * vv) + kdp = 180 / π * 1e-3 * λ * + sum(real(row.Shh_forward - row.Svv_forward) * weights[i] + for (i, row) in pairs(table)) + + (; ZH = 10log10(Zhh), + ZV = 10log10(Zvv), + ZDR = 10log10(Zhh / Zvv), + RHOHV = ρhv, + KDP = kdp) +end + +function valid_radar_moments(moments) + all(isfinite, (moments.ZH, moments.ZV, moments.ZDR, moments.RHOHV, moments.KDP)) +end + +function process_wrf_grid(fields::WRFFields, table, freq_hz; dsd) + shape = size(fields.qrain) + ZH = fill(FILL, shape) + ZV = fill(FILL, shape) + ZDR = fill(FILL, shape) + RHOHV = fill(FILL, shape) + KDP = fill(FILL, shape) + QNRAIN_OUT = isnothing(fields.qnrain) ? fill(FILL, shape) : Float32.(fields.qnrain) + + Threads.@threads for idx in CartesianIndices(fields.qrain) + qr = fields.qrain[idx] + qr > 1e-10 || continue + + rho = air_density(fields.qvapor[idx], fields.pb[idx], fields.p[idx], + fields.theta[idx]) + qrv = qr * rho + qn_density = dsd == "double" ? max(fields.qnrain[idx] * rho, 1e-9) : NaN + params = dsd_parameters(qrv, qn_density, dsd) + dsd == "single" && (QNRAIN_OUT[idx] = Float32(params.N0 / params.Λ / rho)) + + moments = radar_moments(table, freq_hz, params.N0, params.Λ) + if valid_radar_moments(moments) + ZH[idx] = Float32(moments.ZH) + ZV[idx] = Float32(moments.ZV) + ZDR[idx] = Float32(moments.ZDR) + RHOHV[idx] = Float32(moments.RHOHV) + KDP[idx] = Float32(moments.KDP) + else + ZH[idx] = FILL + ZV[idx] = FILL + ZDR[idx] = FILL + RHOHV[idx] = FILL + KDP[idx] = FILL + end + end + + (; ZH, ZV, ZDR, RHOHV, KDP, QNRAIN = QNRAIN_OUT) +end + +function write_output(path, fields::WRFFields, output) + rm(path; force = true) + ds = NCDataset(path, "c") + try + nx, ny, nz, nt = size(fields.qrain) + defDim(ds, "west_east", nx) + defDim(ds, "south_north", ny) + defDim(ds, "bottom_top", nz) + defDim(ds, "Time", nt) + ds.attrib["title"] = "Synthetic dual-pol radar output" + ds.attrib["source"] = "TransitionMatrices.jl radiative_transfer example" + + write_var(ds, "XLAT", Float32.(fields.xlat), ("west_east", "south_north", "Time")) + write_var(ds, "XLONG", Float32.(fields.xlong), ("west_east", "south_north", "Time")) + write_var(ds, "ZNU", Float32.(fields.znu), ("bottom_top",)) + write_var(ds, "XTIME", Float32.(fields.xtime), ("Time",)) + + dims4 = ("west_east", "south_north", "bottom_top", "Time") + write_var(ds, "ZH", output.ZH, dims4; units = "dBZ") + write_var(ds, "ZV", output.ZV, dims4; units = "dBZ") + write_var(ds, "ZDR", output.ZDR, dims4; units = "dB") + write_var(ds, "RHOHV", output.RHOHV, dims4) + write_var(ds, "KDP", output.KDP, dims4; units = "degree km-1") + write_var(ds, "QRAIN", Float32.(fields.qrain), dims4; units = "kg kg-1") + write_var(ds, "QNRAIN", output.QNRAIN, dims4; units = "kg-1") + finally + close(ds) + end +end + +function main(args = ARGS) + opts = parse_args(args) + fields = read_wrf_fields(opts.input; dsd = opts.dsd) + table = load_or_compute_scattering_table(opts.cache, opts.diameters, opts.freq_hz; + solver = opts.solver) + output = process_wrf_grid(fields, table, opts.freq_hz; dsd = opts.dsd) + write_output(opts.output, fields, output) + valid = count(!=(FILL), output.ZH) + println("wrote dual-pol radar output: ", opts.output) + println("valid rain grid cells: ", valid, " / ", length(output.ZH)) +end + +abspath(PROGRAM_FILE) == abspath(@__FILE__) && main() diff --git a/examples/rain_radar.jl b/examples/rain_radar.jl index 9509043..e145990 100644 --- a/examples/rain_radar.jl +++ b/examples/rain_radar.jl @@ -59,7 +59,7 @@ begin denominator = 1 + 2x * sinpi(α / 2) + (λs / λ_m)^(2 - 2α) εreal = ε∞ + (εs - ε∞) * (1 + x * sinpi(α / 2)) / denominator εimag = (εs - ε∞) * x * cospi(α / 2) / denominator + - σ * λ_m / (2π * 2.99792458e8 * ε0) + σ * λ_m / (2π * c0_m_per_s() * ε0) sqrt(complex(εreal, εimag)) end diff --git a/packages/EBCMPrecisionLossEstimators/Project.toml b/packages/EBCMPrecisionLossEstimators/Project.toml index d9e29bf..be74b49 100644 --- a/packages/EBCMPrecisionLossEstimators/Project.toml +++ b/packages/EBCMPrecisionLossEstimators/Project.toml @@ -11,5 +11,5 @@ TransitionMatrices = "057c4241-e127-4181-840e-6b4b92e6eef5" [compat] MLJ = "0.23" MLJXGBoostInterface = "0.3" -TransitionMatrices = "0.4" +TransitionMatrices = "0.5" julia = "1.8"