Skip to content

fix: exact geos inverse + red-band relaxation for limb-correct true color - #1

Merged
LeonJian merged 4 commits into
masterfrom
fix-modifiers-geos-red-parity
Aug 15, 2026
Merged

LeonJian merged 4 commits into
masterfrom
fix-modifiers-geos-red-parity

Conversation

@LeonJian

@LeonJian LeonJian commented Aug 1, 2026

Copy link
Copy Markdown
Owner

Summary

Fixes the insufficient edge contrast in the true_color_reproduction output image (bottom limb mean ~50/255, flat ~11-22 std) by replacing the flat tangent-plane geos projection inverse with the exact ellipsoidal ray–ellipsoid intersection, and completing the Satpy/pyspectral parity details for the Rayleigh correction.

Root cause

rusty_sat_modifiers computed per-pixel lon/lat with a flat-Earth approximation (lat = atan(tan θ_y · cos θ_x), lon = λ₀ + θ_x), exact only at the sub-satellite point. At the AHI full-disk limb (scan angle ~8.4°) it underestimates ground offsets by the h/a curvature factor (~5.6×), giving latitudes/longitudes off by ~57°. Consequences for the true-color pipeline:

pixel true sza / corr old (flat) sza / corr
bottom limb 84.5° / 10.5 29.4° / 1.15
left limb 76.1° / 4.2 26.8° / 1.12
top limb 46.7° / 1.46 15.0° / 1.04

The limb was not sun-zenith-corrected and the Rayleigh LUT path lengths were ~9× too small → dark, flat edges.

Changes

  • geos.rs: exact ellipsoidal geos inverse — ray–ellipsoid intersection (quadratic; A separable, B/C constants) + geodetic-latitude Newton iteration; matches PROJ geos / pyresample get_lonlats, verified to 1e-9 against limb reference values (horizon ±81.28°).
  • angles.rs: rad-based sun_angles_from_lonlat / sat_angles_from_lonlat; AngleSet and the strip loops route through the exact inverse; deleted precompute_columns/rows + solar/satellite_single_precomputed.
  • sun_zenith.rs / rayleigh.rs: strip loops use the exact inverse (space pixels → NaN via discriminant).
  • Rayleigh pyspectral parity:
    • LUT boundary: clamp sun/sat zenith angles to the LUT max secant before computing 1/cos (_clip_angles_inside_coordinate_range) — limb/night pixels evaluate at the LUT edge instead of extrapolating past the last grid point.
    • Correction order: interp → relax(red) → clip [0,100] → reduce → vis − corr (pyspectral np.clip semantics; dropped the result max(0)).
    • Red-band cloud relaxation: new RedBandSource (None / Dataset / SunZenithCorrectedVis), wiring Satpy rayleigh_corrected semantics (sun-zenith-corrected B03 as the red band; B03 itself uses the in-place post-sunz values with no extra memory).
  • true_color_reproduction.rs: wires the B03 red band (0.5 km self + 1 km nearest-resampled via the existing pipeline), replaces the pre-fix artifact assertion (Rayleigh monotonicity does not hold under sunz amplification), and adds edge-contrast regression assertions on the in-memory u8 image (bottom/left limb means > 120, left-limb std > 25, top std > 20, center mean/std unchanged).

Measured impact (true_color_05km.png)

region old mean old std new mean new std
center 112/116/124 54/51/46 111/115/123 55/52/46
bottom limb 47/54/57 18/19/22 169/179/201 22/17/16
left limb 76/87/95 19/16/14 170/181/206 47/41/26
top limb 129/143/151 21/14/12 143/152/162 32/26/23

Tests

  • geos parity table (limb pixels, 1e-9), spherical/ellipsoidal horizon checks
  • limb sunz/correction-factor tests (curved-Earth references)
  • LUT secant clamping, correction clip order, red relaxation (incl. red > 100 negative factor parity, SunZenithCorrectedVis path gating)
  • full workspace: 660/660 tests pass (release), cargo fmt + cargo clippy -D warnings clean

…olor

The geos projection inverse in rusty_sat_modifiers used a flat tangent-plane
approximation (lat = atan(tan theta_y * cos theta_x), lon = lon_0 + theta_x)
that is exact only at the sub-satellite point. At the AHI full-disk limb
(scan angle ~8.4 deg) it underreported ground offsets by the h/a curvature
factor (~5.6x), giving latitudes off by ~57 deg. Sun-zenith correction
factors at the limb were ~1.04-1.15 instead of the true 1.5-10.5, and the
Rayleigh LUT path lengths were ~9x too small, leaving the true-color output
edges flat and dark (bottom limb mean ~50/255).

Replace with the exact ellipsoidal ray-ellipsoid intersection (PROJ geos /
pyresample get_lonlats parity, verified to 1e-9 against limb reference
values), route the sun_zenith/rayleigh strip loops through it, delete the
precompute machinery, and add pyspectral parity details:

- LUT boundary: clamp sun/sat zenith angles to the LUT max secant before
  computing 1/cos (pyspectral _clip_angles_inside_coordinate_range), so
  limb/night pixels evaluate at the LUT edge instead of extrapolating.
- Correction order: relax(red) -> clip [0,100] -> reduce -> vis - corr
  (pyspectral np.clip semantics; drop the result max(0)).
- Red-band cloud relaxation: RedBandSource { None, Dataset, SunZenithCorrectedVis }
  wiring Satpy's rayleigh_corrected semantics (B03 sunz-corrected red; B03
  itself uses the in-place post-sunz values with no extra memory).
- true_color_reproduction: wire the B03 red band (0.5 km self + 1 km
  nearest-resampled), replace the artifact assertion (Rayleigh monotonicity
  does not hold under sunz amplification), and add edge-contrast regression
  assertions (bottom/left limb means > 120, left-limb std > 25, center
  unchanged). Measured: bottom limb mean 47-57 -> 168-200, left-limb std
  19 -> 47.

All 660 workspace tests pass (release), fmt and clippy clean.
…ght blend

- sun_zenith: replace the inverted cosine-domain 88-95 deg gradient with the
  Satpy angle-domain formula (continuous at 88 deg, zero at max_sza), verified
  against _sunzen_corr_cos_ndarray reference values; add daynight_blend_weights
  (Satpy DayNightCompositor cos-zenith weights)
- rayleigh: combined sunz+rayleigh path applies the sunz amplification for
  wavelengths outside the Rayleigh LUT range (e.g. AHI B04 at 0.86 um gets the
  Satpy [sunz_corrected]-only semantics instead of the raw array)
- composites: new DayNightCompositor (fill.rs) blending corrected/uncorrected
  band-major RGB with per-pixel weights, mask OR propagation, owned/borrowed
- true_color_reproduction test: three outputs - standard true_color (hybrid
  green, CIRA), JMA reproduced-green corrected render, and the full JMA
  true_color_reproduction day/night blend (lim 73/85) with black background
  (masked outside-scan pixels become NaN in the uncorrected composite)
…land-color parity

The JMA true_color_reproduction renders applied the standard CIRA stretch, but
Satpy's true_color_reproduction_color_stretch chain is a per-pixel color
conversion matrix (satpy/enhancements/ahi.py, Himawari-8/9) followed by a log
stretch (min 3 / max 150). The missing matrix left land too green (Australia
R-G +8.7 vs +25.0 after the fix).

- rusty_sat_image: JMA_CCM_HIMAWARI_8/9 constants, jma_ccm_for_platform,
  jma_true_color_reproduction_value, log_stretch_value (trollimage
  stretch_logarithmic parity), and the fused finalize_rgb_jma_u8 finalizer
  with hand-computed reference tests
- true_color_reproduction test: the reproduced-green render and the
  uncorrected composite now use finalize_rgb_jma_u8 (Himawari-9); the
  standard true_color keeps the CIRA stretch
…eline

Measured 178s -> ~84s standalone (2.1x), 305s -> 119s under full-suite
concurrency, with byte-identical output PNGs.

- modifiers: apply_corrections_with_sun_zenith_batch corrects several bands
  on one grid in a single pass (angles computed once per pixel; per-band
  in-place red prerequisites; bit-identical to sequential corrections), and
  apply_correction_with_sun_zenith_and_weights emits the DayNightCompositor
  weights during the correction (no second full-disk angle pass); trilinear
  interval lookup switched to binary search and secant clip angles hoisted
- true_color_reproduction test: five 1 km corrections run as one batch
  (33.9s -> 9.2s), the 0.5 km correction fuses the day/night weights
  (24s weights pass removed), and the uncorrected JMA composite reuses raw
  band clones instead of a second Scene load (21.9s saved); per-phase timing
  instrumentation added; all three output PNGs byte-identical to before
@LeonJian

LeonJian commented Aug 8, 2026

Copy link
Copy Markdown
Owner Author

Follow-up: scientific corrections, JMA TCR color parity, and performance (3 new commits)

This PR was opened with the exact-geos-inverse fix (5c7a36c). The branch has since grown three commits that close the remaining scientific gaps found by review, add the missing JMA True Color Reproduction enhancement, and cut the integration-test runtime by ~2.1× while keeping the output bytes identical. All work was verified against the inspected Satpy/pyspectral/trollimage reference code.

1. 3585c4e — SZA gradient fix, out-of-LUT sun-zenith fix, JMA day/night blend

Root cause A (whole-image red cast): the 88°–95° sun-zenith falloff was inverted. sza_correction_factor used a cosine-domain gradient with reversed endpoints: 0 at 88° (should be ~28.65, continuous with the 1/cos branch) and 28.65 at 95° (should be 0). Replaced with the exact Satpy _sunzen_corr_cos_ndarray angle-domain formula; reference checks now lock 88°/90°/94.5° → 28.65/18.26/1.50.

Root cause B (edges redder than the interior): the combined correction skipped the sun-zenith amplification for wavelengths outside the Rayleigh LUT (400–800 nm). AHI B04 (0.86 µm) — the 15% NIR term of hybrid_green — was returned raw, so the green channel was too dark at high SZA and the red cast grew toward the limb. The combined path now applies the sunz amplification unconditionally and the Rayleigh subtraction only for in-LUT wavelengths (Satpy [sunz_corrected]-only semantics for B04).

Also in this commit:

  • daynight_blend_weights (Satpy DayNightCompositor._get_coszen_blending_weights parity) and a new DayNightCompositor in rusty_sat_composites (band-major RGB blend, mask OR propagation, owned/borrowed).
  • The integration test now renders three outputs: standard true_color (hybrid green, CIRA), JMA true_color_reproduction_corr (reproduced green), and the full JMA true_color_reproduction day/night blend (lim_low 73°, lim_high 85°) that replaces the bright gray limb with the natural uncorrected colors.
  • White-background fix for the uncorrected composite: masked outside-scan pixels become NaN (the reader stores bright calibrated values there under the validity mask).

2. 2547b1a — JMA TCR enhancement (color conversion matrix + log stretch)

Root cause (land too green vs the JMA product): the JMA renders used the standard CIRA stretch. Satpy's true_color_reproduction_color_stretch chain is a per-pixel color conversion matrix (Satpy enhancements/ahi.py, Himawari-8/9) followed by a log stretch (min 3 / max 150, trollimage stretch_logarithmic). Added the fused finalize_rgb_jma_u8 finalizer plus JMA_CCM_HIMAWARI_8/9, jma_ccm_for_platform, jma_true_color_reproduction_value and log_stretch_value (trollimage parity), with hand-computed reference tests. Measured effect on Australian land: R−G +8.7 → +25.0 (warm red-brown, matching the JMA TCR look).

3. 6a497f5 — performance (byte-identical output)

Per-phase timing instrumentation showed the runtime split: corrections ~64 s (36%), a separate full-disk day/night-weight pass ~29 s, and two scene loads ~39 s. Optimizations, all verified byte-identical on the three output PNGs:

change before after
apply_corrections_with_sun_zenith_batch — five 1 km corrections share one per-pixel angle pass (bit-identical to sequential) 33.9 s 9.2 s
apply_correction_with_sun_zenith_and_weights — day/night weights emitted during the 0.5 km correction 24–29 s 0
uncorrected composite reuses raw-band clones instead of a second Scene load 21.9 s 0
trilinear interval lookup → binary search; secant clip angles hoisted per correction — —

Total: ~178 s → ~84 s standalone (2.1×), 305 s → 119 s under full-suite concurrency. The three output PNGs are byte-identical to the pre-optimization runs (verified with cmp), and the batch API has a bit-for-bit parity test against the sequential path.

Verification

  • 679/679 workspace tests pass (release); cargo fmt --all -- --check and cargo clippy --workspace --all-features --all-targets -- -D warnings clean.
  • New tests: Satpy angle-domain SZA reference values, out-of-LUT combined sunz-only parity, trilinear bit-parity vs rustyspectral, batch-vs-sequential bit-parity, JMA CCM/log-stretch hand-computed references, day/night blend semantics (day side == corrected, night side == uncorrected, gray-limb reduction, deep-night rescue).

@LeonJian
LeonJian merged commit 93730af into master Aug 15, 2026
9 checks passed
@LeonJian
LeonJian deleted the fix-modifiers-geos-red-parity branch August 15, 2026 14:20
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant