User guide

Install

The core install runs the entire algorithm on fields already in memory, and pulls in no I/O or plotting stack:

pip install shachen

Optional extras:

pip install "shachen[satellite]"   # satpy: read and calibrate ABI/AHI L1b
pip install "shachen[data]"        # earthaccess/s3fs: fetch MERRA-2, CAMEL, L1b
pip install "shachen[render]"      # cartopy/matplotlib: georeferenced PNG
pip install "shachen[all]"         # everything, for the reproduction scripts

import shachen pulls in none of the extras; CI asserts that.

The scene contract

Every stage works on the sensor’s 2 km fixed grid. A scene is an xarray.Dataset whose variables are named by DEBRA role rather than by sensor channel, so ABI and AHI share one code path:

Variable	Unit	Role
refl_vis_064	%	daytime baseline image (Eq. 23)
refl_nir_160	%	reserved for the daytime cloud test
bt_swir_39	K	night thin-cirrus test CM4 (Eq. 6)
bt_wv_62	K	deep-convection test CM2 (Eq. 4)
bt_tir_86	K	dust test DT2 (Eq. 14)
bt_tir_104	K	clean window reference
bt_tir_123	K	dirty window, RSW / DT1 (Eq. 13)

Two attributes must be present:

scene.attrs["area"]

a pyresample.geometry.AreaDefinition — the projection and grid every ancillary field is interpolated onto, and the map projection render draws in.

scene.attrs["start_time"]

a naive-UTC datetime.datetime, the scan start. Only the solar zenith angle uses it, and the terminator blend is ~30° wide, so scan start is accurate enough.

shachen.io.satellite.load_scene() builds exactly this from ABI netCDF or AHI HSD files. The mapping from role to channel lives in shachen.constants.ABI_BANDS and shachen.constants.AHI_BANDS; if you read L1b some other way, produce the same variable names and you can skip the satellite extra entirely.

Running the algorithm

shachen.pipeline.run_debra() chains background → cloud mask → dust tests → confidence and returns every intermediate field alongside the answer:

import shachen

result = shachen.run_debra(scene, skin_temperature=merra_ts, emissivity=camel)

result["cf_comb"]  # combined dust confidence, 0–1 (Eq. 22)
result["dt1"]  # split-window test against the dynamic background
result["cm_norm_day"]  # normalised daytime cloud mask
result["zenith_deg"]  # per-pixel solar zenith

skin_temperature is MERRA-2 TS (K) on its native lat/lon grid; the pipeline regrids it via shachen.geo.regrid_latlon(). Pixels with NaN inputs (off-disk, bad pixels) carry NaN confidence throughout.

Running ZHOUYE

shachen.pipeline.run_zhouye() (昼夜, “day and night”) runs the same chain with the shachen.constants.ZHOUYE preset. It changes three things, all in the terminator and night confidence sums (Eqs. 17–18) and none in the daytime sum: the 8.6 − 10.4 µm test is read on a fixed interval instead of Eq. 14 (dt2_fixed, emitted next to Eq. 14’s dt2), a pixel lit by the thermal-contrast test alone gets no confidence, and the night and terminator normalisation intervals are fitted so that dust reads the same confidence as by day. Below a solar zenith of 75° the output equals run_debra with ABI_TUNED exactly; beyond it the reading on dust stays at the day’s level. The scheme, the ground truth behind it and what it costs are described in a paper in preparation.

result = shachen.run_zhouye(scene, skin_temperature=merra_ts, emissivity=camel)
result.attrs["scheme"]  # "zhouye"

The three hooks are off in DEFAULTS and ABI_TUNED, so DEBRA’s output is untouched.

Choosing a background

The clear-sky background is estimated per pixel, by one of two schemes. Exactly one must be given: passing both, or neither, raises ValueError.

Scheme A — semi-analytic (paper §3.2, Option A). Per-band surface emissivity modifies the Planck radiance of the MERRA-2 skin temperature, and that is inverted back to a brightness temperature. Pass emissivity=, a CAMEL Dataset on its native grid:

result = shachen.run_debra(scene, skin_temperature=merra_ts, emissivity=camel)

Scheme B — cloud-cleared composite (Option B). Stack the same time-of-day scenes from the preceding ~14 days and keep, per pixel, the bands from the day with the warmest bt_tir_104. Build it yourself and pass background=, already on the scene grid:

from shachen.composite import composite_background

bg = composite_background([day_1, day_2, ...])  # same grid, same hour
result = shachen.run_debra(scene, skin_temperature=merra_ts, background=bg)

Scheme B is built from real observations, so it carries the split-window water-vapour depression that the atmosphere-free scheme A lacks: the ~1 K high bias that zeroes DT1/DT2 on transparent winter plumes. Its per-pixel candidate count is passed through to the output as n_valid.

Any Dataset carrying rsw_bg, btd_bg and bt_bg_tir_86/104/123 on the scene grid works as background=; composite_background and shachen.background.background_signals() both produce that contract.

Tuning

shachen.constants.DEFAULTS is the paper-tuned constant set, and it is the single source of every bound, offset and weight — each one appears in Equations in the formula it belongs to. Pass your own shachen.constants.DebraConstants to retune:

from shachen.constants import ABI_TUNED

result = shachen.run_debra(..., constants=ABI_TUNED)

shachen.constants.ABI_TUNED is the only shipped retune. It raises the Eq. 19 lower bound from 0.25 to 0.40 to suppress a clear-sky DT3 floor specific to the ABI + MERRA-2 + CAMEL stack. The reasoning and the numbers are in Deviations.

Eq. 19’s interval is per branch: cf_norm_day and cf_norm_ngt, with the terminator interpolating between them. The night interval is the day one halved, because Eq. 18’s raw sum peaks at 1.5 against Eq. 16’s 3.0; with one shared interval cf_ngt could never exceed 0.556 and dust faded at dusk. Passing cf_norm= sets both at once and restores the printed single-interval behaviour:

from shachen.constants import Bounds, ConfidenceConstants, DebraConstants

published = DebraConstants(
    confidence=ConfidenceConstants(cf_norm=Bounds(0.25, 2.50)),
)

Enhanced imagery

The confidence field is rendered in two steps (Eqs. 23–29): a day/night blended greyscale baseline, then a per-gun colour modulation by cf_comb.

from shachen.imagery import debra_imagery, to_uint8
from shachen.render import render_debra_png

imagery = debra_imagery(scene, result)  # vis_bg, ir_bg, b_bg, bi, rgb
render_debra_png(
    to_uint8(imagery["rgb"]),
    scene.attrs["area"],
    scene.attrs["start_time"],
    Path("dust.png"),
)

Dust is yellow by default. The paper’s other presets live in shachen.constants.COLOR_DIMMING and are selected with the dimming= argument:

imagery = debra_imagery(scene, result, dimming=COLOR_DIMMING["pink"])

render_debra_png is the only function that needs cartopy and matplotlib; both are imported inside it, so the rest of the package stays importable without the render extra.

Reproducing a reference case

End to end reproduction needs [all] plus an Earthdata login in ~/.netrc: MERRA-2 and CAMEL are authenticated downloads, while GOES L1b on AWS S3 is anonymous:

python scripts/fetch_case.py 2017-03-23-swus   # the paper's Figure 6 case
python scripts/run_case.py   2017-03-23-swus   # → netCDF + PNG