The user wants a quality control (QC) script for daily rainfall values in Python. This is a common task in hydrology/meteorology. Let me think about what quality checks are typically performed on daily rainfall data: 1. **Missing values check** - identify NaN or missing data codes 2. **Range/physical limits check** - daily rainfall typically between 0 and some maximum (e.g., 0-1000mm, or region-specific like 0-300mm for most climates) 3. **Negative values check** - rainfall can't be negative 4. **Repeated values check** - many consecutive identical non-zero values may indicate sensor stuck 5. **Dry/wet spell consistency** - extremely long wet or dry spells 6. **Spatial/temporal consistency** - comparing with nearby stations (if available) 7. **Outlier detection** - statistical outliers (e.g., using IQR, z-score, or boxplot method) 8. **Jump checks** - sudden large changes between consecutive days 9. **Zero accumulation** - too many zeros followed by large values might indicate accumulation issues in tipping bucket gauges 10. **World Records check** - comparing against world record daily rainfall (1825 mm at Cherrapunji/ Foc-Foc Réunion 1825mm) 11. **Monthly totals check** - compare against climatological maximum monthly totals Common QC frameworks: - ETCCDI QC procedures - RClimDex / ClimDex QC - WMO guidelines - DQCD (Data Quality Control) from WMO Key tests from RClimDex-style QC for precipitation: - Missing value check - Identical values repeated ≥ N days (e.g., 3 days of same nonzero value) - Outlier check using boxplot (values outside median ± n*IQR) - Value > physical world record (e.g., > 2000 mm/day flag) - Wet spell: too many consecutive wet days - Dry spell: too many consecutive dry days with nonzero traces? Actually dry spell check is for temperature mostly. For precipitation: wet spell length check. Let me write a well-structured Python script with: - pandas-based - A function per test, each returning a boolean flag Series - A summary DataFrame of flags - Configurable parameters - Example usage with synthetic data - Clean code with docstrings Let me design: ```python import numpy as np import pandas as pd def qc_missing(series, missing_codes=[-999, -99]): ... def qc_negative(series): ... def qc_upper_limit(series, upper=1000): # world record ~1825mm def qc_repeated_values(series, max_repeat=3): def qc_outliers_iqr(series, factor=...): def qc_wet_spell(series, max_wet_days=30): def qc_dry_spell(series, max_dry_days=365): # maybe optional def qc_jump(series, threshold=...): # day-to-day jump def qc_monthly_total(series, max_monthly=...): ``` Flag codes convention: use -1 for missing, 0 for pass, 1 for fail/suspect. Or True/False booleans. I'll use flag conventions: 0 = OK, 1 = suspect/error, -1 = missing. Let me write the code with a class or a pipeline function `run_qc(df, column)` that returns a DataFrame with flag columns and a summary report. Also add optional spatial consistency check comparing to neighbor stations — mention it as advanced. Let me write clean, complete, runnable code with example data generation and printed report. Include comments and docstrings. Keep reasonably concise but comprehensive. Structure: 1. Brief intro of the checks 2. Code 3. Example usage/output 4. Notes on tuning parameters I'll write it with pandas. Let me draft: ```python """ Daily precipitation quality control (QC) """ import numpy as np import pandas as pd # ----------------------- configuration ----------------------- MISSING_CODES = [-999.0, -99.0, -9999.0] UPPER_LIMIT = 1000.0 # mm/day (physically plausible; world record 1825 mm) MAX_REPEATED = 3 # max consecutive days with identical non-zero value IQR_K = 25 # boxplot outlier factor (RClimDex uses std*10 historically; boxplot k) MAX_WET_SPELL = 30 # consecutive wet days (nonzero) before flagging MAX_DRY_SPELL = 365 # consecutive dry days before flagging JUMP_THRESHOLD = 150.0 # mm/day-to-day change flagged as suspicious MAX_MONTHLY_TOTAL = 4000 # mm/month (Cherrapunji ~ 26000 mm/yr, monthly record 3660 mm? Actually July 1861 Cherrapunji 9300 mm) ``` Hmm, monthly record: Cherrapunji recorded 9300 mm in July 1861. So max monthly total could be ~10,000 mm for extreme monsoon regions, but for typical QC, use like 3500-4000 mm default configurable. I'll set default 5000 and note region-specific tuning. Outlier test: RClimDex uses values outside mean ± 10 * std? Actually RClimDex uses boxplot: outliers are > Q3 + 1.5*IQR or < Q1 - 1.5*IQR, extreme outliers 3*IQR. But for precipitation, most days are 0, so IQR is 0 → everything nonzero flagged. So the outlier test must only apply to wet days, or use a different approach. Common approach: apply boxplot to wet-day amounts only (days > 0), or use percentile-based thresholds (e.g., > 99.9th percentile of wet days). I'll implement: outliers among wet days using IQR on wet-day amounts, plus optionally extreme threshold via lognormal z-score. Better: compute on wet days only. Also the "jump" check: big difference between consecutive days. Also "stuck at zero"? Not needed. Also repeated zeros are normal, so repeated-value test must exclude zeros. Let me code carefully: ```python def qc_daily_rainfall(series, cfg): df = pd.DataFrame({"prcp": series}) flags = pd.DataFrame(index=df.index) # 1. missing flags["missing"] = df["prcp"].isin(cfg["missing_codes"]) | df["prcp"].isna() # 2. negative flags["negative"] = df["prcp"] < 0 # 3. above physical limit flags["above_limit"] = df["prcp"] > cfg["upper_limit"] # 4. repeated identical nonzero values v = df["prcp"].where(flags["missing"], np.nan) rep = v.ne(v.shift()) # group id grp = rep.cumsum() run_len = v.groupby(grp).transform("size") flags["repeated"] = (run_len > cfg["max_repeated"]) & (v > 0) # 5. outliers among wet days (IQR method) wet = v[v > 0] ... flags["outlier"] = False if wet.size >= cfg["min_wet_days"]: q1, q3 = wet.quantile([0.25, 0.75]) iqr = q3 - q1 hi = q3 + cfg["iqr_k"] * iqr # If iqr == 0 (e.g., lots of identical small values), fall back to percentile ... flags["outlier"] = v > hi # 6. wet spell is_wet = (v > 0) & ~flags["missing"] wet_run = ... consecutive count flags["wet_spell"] = is_wet & (wet_run > cfg["max_wet_spell"]) # 7. dry spell is_dry = (v == 0) & ~flags["missing"] dry_run... flags["dry_spell"] = is_dry & (dry_run > cfg["max_dry_spell"]) # 8. jump diff = v.diff().abs() flags["jump"] = diff > cfg["jump_threshold"] # combine: overall flag flag_cols = [...] df["n_flags"] = flags[flag_cols].sum(axis=1) ... ``` Consecutive run counting: standard trick: ```python def _run_length(mask): # length of the current run of True values, per element mask = mask.astype(int) # reset at False grp = (mask == 0).cumsum() run = mask.groupby(grp).cumsum() + ... ``` Actually simpler: ```python def consecutive_count(mask): """For each element of boolean mask, count of consecutive Trues ending at that element? Or length of the run containing the element.""" ``` For spell checks, we want the run length of the run each element belongs to. Implementation: ```python def _run_lengths(mask): m = mask.to_numpy(dtype=bool) run_lengths = np.zeros(len(m), dtype=int) i = 0 while i < len(m): if m[i]: j = i while j < len(m) and m[j]: j += 1 run_lengths[i:j] = j - i i = j else: i += 1 return pd.Series(run_lengths, index=mask.index) ``` Or vectorized pandas trick: ```python def _run_lengths(mask): # id each run, count size, broadcast back new_run = mask & ~mask.shift(fill_value=False) run_id = new_run.cumsum() return mask.groupby(run_id).transform("sum").where(mask, 0) ``` Wait: groupby on run_id where mask False would group all False into run_id groups too... If mask False, new_run False, run_id stays same as previous run — all Falses after a run share that run's id, transform("sum") would count the True values in that group which equals the run length, then .where(mask, 0) sets to 0. Hmm, but at the very beginning, leading Falses: new_run False, cumsum 0 → group 0 all False → sum 0. That works. But a later block of Falses belongs to run_id of previous run, sum over group = length of previous run, then where(mask,0) → 0. OK fine. But careful: groupby with run_id groups the run Trues + following Falses; transform sum counts Trues only = run length. Good. Edge: mask all False → new_run all False → run_id all 0 → sum 0. Good. Alternatively simpler with cumcount approach. I'll write helper `_spell_lengths`. Missing values handling in runs: treat missing as breaking the run (conservative) or as neutral? For wet spell, missing should probably break the run. With v = prcp where not missing else NaN, `v > 0` gives False for NaN, so missing breaks runs. Good. Jump check: diff across a missing gap might be flagged spuriously; better to only compute jump when both days valid: ```python valid = v.notna() flags["jump"] = valid & valid.shift(1, fill_value=False) & (v.diff().abs() > threshold) ``` Also should exclude jump flag when either day already flagged as error? Keep simple but mention. Also trace values: some datasets encode trace as 0.1 or -1; mention optionally converting traces to 0. Summary report function: counts per flag, list affected dates. Final combined status column: - "missing" - "fail" (any hard error: negative, above limit, repeated) - "suspect" (outlier, jump, spell) - "ok" Let me also provide example usage generating synthetic data with injected errors, and show printed summary. Also plotting optional? Could include a simple plot function marking flagged days — nice but optional; I'll include a small optional matplotlib snippet maybe. Keep the answer focused; maybe include a brief plotting function since QC review often visual. I'll add a short optional plotting section. Formatting: The user asked to "write a quality control of daily rainfall values in python". So they want code. Provide: - Short intro listing checks - Full code block - Example usage + sample output - Notes on tuning/extending (spatial check, neighbor stations, homogenization) Use headers and code blocks. Reasonably complete but not bloated. The code will be the bulk. Use markdown since it's code-heavy chat context. Let me now write the final code carefully, mentally checking correctness. ```python """ qc_rainfall.py Quality control (QC) for daily rainfall time series. Checks implemented ------------------ 1. Missing values -> flag 'missing' 2. Negative values -> flag 'negative' (hard error) 3. Values above physical max -> flag 'above_limit' (hard error) 4. Repeated identical non-zero values over N+ days -> 'repeated' (stuck gauge) 5. Statistical outliers (boxplot/IQR on wet days) -> 'outlier' 6. Wet spell: too many consecutive wet days -> 'wet_spell' 7. Dry spell: too many consecutive dry days -> 'dry_spell' 8. Day-to-day jump larger than threshold -> 'jump' Flag convention per check: True = failed the check. """ import numpy as np import pandas as pd DEFAULT_CFG = dict( missing_codes=(-999.0, -99.0, -9999.0), upper_limit=1000.0, # mm/day; world record ~1825 mm (Foc-Foc, Réunion, 1966... ) max_repeated=3, # >3 consecutive equal non-zero values flagged (i.e. 4+)? define as flag if run_len > max_repeated iqr_k=3.0, # extreme outlier multiplier on IQR (3*IQR = extreme boxplot fence) min_wet_days=20, # min number of wet days before outlier test is applied max_wet_spell=30, # consecutive wet days max_dry_spell=180, # consecutive dry days jump_threshold=150.0, # mm change between consecutive days ) def _spell_lengths(mask: pd.Series) -> pd.Series: """Length of the run of consecutive True values each element belongs to.""" new_run = mask & ~mask.shift(fill_value=False) run_id = new_run.cumsum() return mask.groupby(run_id).transform("sum").where(mask, 0) ``` Wait — groupby(run_id) where run_id has 0 for leading falses... `mask.groupby(run_id)` — groupby keys are run_id values; fine. transform("sum") on boolean series returns... boolean groupby sum returns int. ok. `.where(mask, 0)` — for elements where mask False set 0. But careful: if mask is boolean and where with scalar 0 fine. Edge: run_id for trailing falses equals last run id; sum within that group includes the trues of that run; then set to 0 via where. Good. But there's a subtlety: `mask.groupby(run_id)` requires run_id same index. Yes. Alternatively simpler: use itertools/np. I'll keep the pandas trick but maybe write explicit numpy loop for clarity? The pandas trick is fine and fast. Now repeated values: ```python def check_repeated(v: pd.Series, max_repeated: int) -> pd.Series: wet = v > 0 new_run = wet & ~wet.shift(fill_value=False) run_id = new_run.cumsum() run_len = v.groupby(run_id).transform("size") return wet & (run_len > max_repeated) ``` Here run length includes all days in the run (since v notna? if missing NaN inside? v is the masked series with NaN for missing; groupby transform size counts NaN too? `transform("size")` counts group size including NaNs. Hmm. If there is NaN inside a run of identical values, that would count. Better to group on v itself? Let's compute run on wet mask but count only non-missing: `run_len = v.where(wet).groupby(run_id).transform("count")`? Actually simpler: count occurrences of the value: since within a run values are identical non-zero, but a NaN breaks? v.diff().ne(0) approach includes NaN handling complexities. Simplest robust: define s = v (NaN for missing). identical = s.eq(s.shift()) & wet & wet.shift(fill_value=False). Hmm. Alternative approach: build runs on the series values directly: ```python def check_repeated(v, max_repeated): wet = (v > 0) if not wet.any(): return pd.Series(False, index=v.index) run_id = (v.ne(v.shift()) | ~wet).cumsum() # new run when value changes OR when not wet run_len = v.groupby(run_id).transform("count") # count non-NaN in group return wet & (run_len > max_repeated) ``` Hmm, v.ne(v.shift()) is True for NaN? NaN != NaN is True in pandas? Actually pandas `ne` treats NaN != NaN as True? Let's think: Series.ne compares with pandas semantics where NaN != NaN returns True (since it uses array comparison where NaN equals nothing). In pandas, `pd.Series([np.nan]).ne(pd.Series([np.nan]))` returns True. Yes, NaN != NaN → True. So a missing value inside a run would start a new run id. Good — missing breaks the run. Also `| ~wet` ensures each dry/missing day gets its own run id (since ~wet True → increment). And wet days with same value continue same run. run_len via transform("count") counts non-NaN entries in the group = length of the run of equal wet values. But note dry days with v=0 all share... each dry day gets its own id due to ~wet. Fine. And wet run where all values equal: transform count = run length. But if a wet run has varying values, each change increments id — fine. But wait: could two separate runs of same value get same run_id? ids increment at every change or dry day, so no. Edge: leading values: v.shift() gives NaN, v.ne(NaN) → True (non-NaN != NaN True; NaN != NaN True). So first day always new run. Good. So `check_repeated` returns wet & run_len > max_repeated. Should repeated zeros count? No — zeros normal. Also should the flag mark the entire run or just days beyond the k-th? RClimDex flags the whole run I think. I'll flag whole run; mention you can keep first N days. Outlier check on wet days: ```python def check_outliers(v, iqr_k, min_wet_days): flag = pd.Series(False, index=v.index) wet = v[v > 0].dropna() if wet.size < min_wet_days: return flag q1, q3 = wet.quantile([0.25, 0.75]) iqr = q3 - q1 if iqr == 0: # degenerate (e.g. all wet days equal); fall back to percentile fence upper = wet.quantile(0.999) ... ``` Hmm fallback: if all wet days same value (e.g., many 0.1mm traces?) IQR=0 → any value > q3 flagged. That could be wrong. Use fallback: flag values > q3 (strictly greater) only if also > some percentile? Simpler fallback: use upper fence = max(q3 + k*iqr, wet.quantile(0.99))? Hmm. Let me instead do: if iqr == 0, use upper = q3 * something? Maybe skip: only flag if value > q3 when iqr>0 else flag values > 99.9th percentile when there are enough wet days (>=100?). Simplify: ```python if iqr > 0: upper = q3 + iqr_k * iqr else: upper = wet.quantile(0.999) # very conservative flag = v > upper ``` But 0.999 quantile of wet days with 20 wet days is just max → flags the max. Hmm, with few wet days, that flags the single max which may be a real storm. Maybe with iqr==0 case, just don't flag (or use rclimdex alternative: mean + 10*std of wet days? RClimDex uses boxplot on... Actually RClimDex uses the "boxplot" method with quantiles 0.05/0.95?? Let me recall: RClimDex QC outlier test for precipitation uses boxplot: quartiles Q1, Q3, and flags values > Q3 + k*IQR where user chooses k (default 1? they say "step 3: outlier detection with boxplot, user decides number of IQRs"). Values below Q1 - k*IQR don't apply to precipitation (can't be negative below 0... well negative already flagged). Also there's the ETCCDI-style: compare with log of wet days? Not needed. I'll implement IQR with configurable k, fallback when IQR==0: use percentile approach: flag if value > wet.quantile(0.995) and wet.size >= 100 else no flags. Eh, maybe simpler: if iqr == 0, no outlier flags (log a note). That's honest and simple. I'll do that but mention in comments. Actually a more robust common alternative: flag wet-day values exceeding mean + k*std of log-transformed wet-day amounts (lognormal). Could offer as alternative method: "method" param: 'iqr' or 'zlog'. Maybe overkill. Keep IQR + note. Wet/dry spells with helper. Jump check: ```python def check_jump(v, threshold): valid = v.notna() jump = v.diff().abs() return valid & valid.shift(fill_value=False) & (jump > threshold) ``` Monthly total check: maybe add: ```python def check_monthly_totals(v, max_monthly): m = v.resample("MS").sum(min_count=1) # careful NaN handling bad_months = m > max_monthly # map back to daily month_key = v.index.to_period("M") flag = pd.Series(pd.Series(bad_months.values, index=m.index).reindex(month_key.astype(...))) ``` Simpler: ```python monthly = v.groupby(v.index.to_period("M")).sum(min_count=1) bad = monthly[monthly > max_monthly].index flag = pd.Series(v.index.to_period("M").isin(bad), index=v.index) & v.notna() ``` That flags all days in bad months. Include as optional check 'monthly_total' with default None (disabled) since threshold is region-specific. Good idea: allow cfg["max_monthly_total"] = None to disable. Also add a check for duplicated index dates? Data with duplicate timestamps. Could add "duplicate_index" flag. Nice touch. Also missing days in the series (gaps in the date index)? That's "date completeness": reindex to full daily range and mark missing dates as missing data. I'll include: at start, reindex to daily frequency so QC covers gaps, flagging them as missing. Good practice. Implement in main function: ```python def run_qc(series, cfg=None): cfg = {**DEFAULT_CFG, **(cfg or {})} s = series.copy() s.index = pd.to_datetime(s.index) s = s.sort_index() # full daily calendar full_idx = pd.date_range(s.index.min(), s.index.max(), freq="D") missing_dates = full_idx.difference(s.index) s = s.reindex(full_idx) ... ``` Then flags DataFrame with columns. Also duplicate index check: `s.index.duplicated()` before reindex — flag duplicated timestamps. Let me add that: if duplicates exist, keep first and note. Keep simple: raise or warn? I'll record a boolean and aggregate duplicates by keeping first, mention in report. Final classification: ```python hard = flags[["negative", "above_limit", "repeated"]].any(axis=1) soft = flags[["outlier", "wet_spell", "dry_spell", "jump", "monthly_total"]].any(axis=1) status = np.select( [flags["missing"], hard, soft], ["missing", "fail", "suspect"], default="ok", ) ``` Order: missing first. Also if missing, other flags likely True spuriously (e.g., repeated NaN? we handled with wet masks; jump requires both valid; outlier v>upper with NaN → False). OK. Report function: counts + flagged dates listing (maybe first N). Provide `report` printing. Then example: synthetic data 3 years, mostly zeros with random wet days gamma-distributed, inject errors: a negative value, a 5000 mm value, 5 days of 2.0 repeated, a jump 0→200→0, and drop some dates. Print report and show flagged dates. Let me write the final code in one block. Also include `if __name__ == "__main__":` demo. Let me also double check `_spell_lengths`: ```python def _spell_lengths(mask): mask = mask.fillna(False).astype(bool) new_run = mask & ~mask.shift(fill_value=False) run_id = new_run.cumsum() return mask.groupby(run_id).transform("sum").where(mask, 0) ``` Wait, groupby on run_id with boolean mask: transform("sum") — sum of booleans → int (pandas may warn?). Boolean series sum in groupby returns int64. ok. `.where(mask, 0)` returns int series with 0 at False positions. But leading False group run_id=0: fine. However, there's an issue: if mask is all False except... fine. Hmm, but groupby(run_id) — run_id includes 0 group for the leading Falses; but what if mask[0] is True? Then new_run[0]=True, run_id starts at 1. fine. Wet spell check: is_wet = v > 0 (NaN > 0 → False, ok). spell = _spell_lengths(is_wet); flag = is_wet & (spell > max_wet_spell). Note is_wet excludes missing and zero. Good. Dry spell: is_dry = (v == 0); NaN == 0 → False, so missing days break dry spells — they are neither wet nor dry. Acceptable; mention. Monthly totals: use groupby period sum with min_count=1 so all-missing month → NaN, not flagged. ```python def check_monthly_total(v, max_monthly): if max_monthly is None: return pd.Series(False, index=v.index) per = v.index.to_period("M") totals = v.groupby(per).sum(min_count=1) bad = totals.index[totals > max_monthly] return pd.Series(per.isin(bad), index=v.index) & v.notna() ``` `per.isin(bad)` works on PeriodIndex. Good. Duplicate index check: before reindex: ```python dup = series.index.duplicated(keep="first") ``` If any, we keep first occurrence and report count. In flags we can add column "duplicate_ts" aligned to kept index? Duplicates dropped, so mapping is awkward. Simpler: raise a ValueError? Or warn and keep first. I'll handle: detect duplicates, print warning in report (return in a QC summary dict of "metadata issues"). To keep code simpler, I'll just drop duplicates and add a note in report. Let me implement `run_qc` returning (qc_df, summary) maybe? Or return a DataFrame with flag columns + status, and a separate summary Series/printed report function `summarize(qc_df)`. I'll structure: - `run_qc(prcp, cfg=None) -> pd.DataFrame` with columns: prcp, missing, negative, above_limit, repeated, outlier, wet_spell, dry_spell, jump, monthly_total, n_flags, status - `qc_report(qc_df, max_list=10) -> str` prints summary. Demo data generation: ```python rng = np.random.default_rng(42) idx = pd.date_range("2020-01-01", "2022-12-31", freq="D") n = len(idx) rain = np.where(rng.random(n) < 0.35, rng.gamma(0.6, 8.0, n), 0.0).round(1) rain[100] = -3.0 # negative rain[200] = 5000.0 # above physical limit rain[300:305] = 2.0 # stuck sensor (5 identical days) rain[400] = 0.0; rain[401] = 180.0 # jump rain[500:520] = 5.0 # 20-day wet spell with identical values → also repeated # remove 10 dates to create gaps drop = rng.choice(n, 10, replace=False) s = pd.Series(rain, index=idx) s.iloc[drop] = np.nan? # or delete ``` Better to delete dates to test gap handling: `s = s.drop(s.index[drop])`. But careful the injected error indices shift? drop after injection fine; injected positions may coincide with dropped but fine. Actually rain[500:520]=5.0 for 20 days would be flagged repeated too (same value 5 days > 3). Maybe vary: rain[500:520] = rng.gamma small values ~ make wet spell of 25 days: choose wet spell where values vary. Let's do: ```python rain[500:530] = rng.gamma(0.8, 1.5, 30).round(1) # 30-day wet spell ``` max_wet_spell=30 → run length 30, flag requires > 30 → not flagged. Set max_wet_spell=20 default? Hmm realistic: 30 consecutive wet days is not that rare in monsoon. Default max_wet_spell maybe 30 flag if > 30 i.e., ≥31 days. For the demo to trigger, use 35-day spell. Or set cfg max_wet_spell=25 in demo. I'll make spell 35 days in demo with default cfg (flag if run > 30). Fine. Jump: default threshold 150 mm. Demo: 0 → 160 mm. Also dry spell: 200 days of zeros → default max_dry_spell 180 → flagged. In a demo with random data, long dry spells may naturally occur; ok. IQR outlier: gamma(0.6, 8) wet amounts: q1, q3, IQR; k=3 → fence ~ q3+3*IQR. Wet values gamma shape 0.6 scale 8: mean ~4.8, heavy tail; max over ~700 wet days maybe ~60-80mm; fence maybe ~30-40 → some natural days flagged as "suspect" — that's realistic but demo output might list many. To keep demo clean, maybe use iqr_k=5 in demo? Or just let it flag some — realistic. But repeated/stuck injection: 5 days of 2.0 — 2.0 likely below fence. Fine. Note interplay: outlier fence computed from data including errors (5000 inflates q3?) — 5000 is a single point; quantiles robust. OK. Mention that iterative QC (remove flagged, recompute) is an option. Now the report function: ```python def qc_report(qc, max_dates=8): lines = ["Daily rainfall QC report", "-" * 40] n = len(qc) lines.append(f"Records evaluated : {n}") for col in FLAG_COLS: c = int(qc[col].sum()) if c: dates = qc.index[qc[col]][:max_dates].strftime("%Y-%m-%d").tolist() more = " ..." if c > max_dates else "" lines.append(f"{col:<14}: {c:5d} e.g. {', '.join(dates)}{more}") lines.append(f"status counts: {qc['status'].value_counts().to_dict()}") return "\n".join(lines) ``` Flag column order list: ["missing","negative","above_limit","repeated","outlier","wet_spell","dry_spell","jump","monthly_total"]. n_flags = sum of soft+hard excluding missing? Just sum all flag columns except missing maybe; keep sum of all except missing. status via np.select as above. Now write everything cleanly. Also type hints, docstrings. Should I include plotting? A small optional function `plot_qc` with matplotlib marking failures — nice but adds length. I'll include a compact version. Length/format: The user asked "write a quality control of daily rainfall values in python" — deliver code-centric answer. Include: - intro sentence - checks list (bulleted, short) - code block (main) - demo usage code (could be in same block with __main__) - sample output block (short) - tuning notes + extensions (short bullets): spatial checks with neighbor stations, traces, region-specific thresholds, iterative QC, flag conventions (0/1/-1), GHCN style. Keep total moderately long but mostly code. Avoid excessive prose. Let me finalize the code text carefully, making sure it runs correctly in my head. Potential pitfalls: - `mask.shift(fill_value=False)` works for boolean series → returns object? shift on bool series with fill_value=False returns bool dtype. Yes. - `v.ne(v.shift())` with NaN: pandas returns True where either NaN. Actually let me verify: Series.ne uses `ops` comparison; NaN == NaN in pandas comparison methods: `pd.Series([1.0, np.nan]).eq(pd.Series([1.0, np.nan]))` → [True, False] because eq treats NaN as unequal. Yes, pandas comparison functions return False for NaN==NaN, True for NaN!=NaN. Good. - groupby transform("count") counts non-NaN. Good. - `s.index.to_period("M")` needs DatetimeIndex — ensured after reindex with date_range. Good. - reindex: after `s = s.reindex(full_idx)`, missing dates NaN → flagged missing. Also original NaNs flagged missing. Good. - For missing_codes: treat them as NaN: `s = s.replace(list(cfg["missing_codes"]), np.nan)` and also `s = s.mask(s < 0 ...)`? Negative handled separately; replace negative with NaN after flagging? For downstream computations (outlier etc.), negative values pollute; but flags computed on raw v; wet masks use v>0 so negatives excluded. jump diff includes negatives — fine, it's flagged anyway. Keep raw. - Missing check: `flags["missing"] = s.isna() | s.isin(missing_codes)` — but if we replaced codes with NaN earlier, isna covers. Let me do: `raw = s.copy()`; `s_clean = raw.replace(missing, np.nan)`; flags["missing"] = s_clean.isna(). Compute everything on s_clean. Good. - `check_repeated` uses s_clean. - n_flags = flags.drop(columns=["missing"]).sum(axis=1)? If missing, other flags false anyway mostly; fine to sum all. - status: np.select(condlist, choicelist, default) returns array; assign. Write `run_qc`: ```python def run_qc(prcp: pd.Series, cfg: dict | None = None) -> pd.DataFrame: cfg = {**DEFAULT_CFG, **(cfg or {})} s = prcp.copy() s.index = pd.to_datetime(s.index) if s.index.has_duplicates: s = s[~s.index.duplicated(keep="first")] s = s.sort_index() s = s.reindex(pd.date_range(s.index.min(), s.index.max(), freq="D")) missing_codes = list(cfg["missing_codes"]) v = s.replace(missing_codes, np.nan) if missing_codes else s.copy() flags = pd.DataFrame(index=v.index) flags["missing"] = v.isna() flags["negative"] = v < 0 flags["above_limit"] = v > cfg["upper_limit"] flags["repeated"] = check_repeated(v, cfg["max_repeated"]) flags["outlier"] = check_outliers(v, cfg["iqr_k"], cfg["min_wet_days"]) flags["wet_spell"] = (v > 0) & (_spell_lengths(v > 0) > cfg["max_wet_spell"]) flags["dry_spell"] = (v == 0) & (_spell_lengths(v == 0) > cfg["max_dry_spell"]) flags["jump"] = check_jump(v, cfg["jump_threshold"]) flags["monthly_total"] = check_monthly_total(v, cfg["max_monthly_total"]) hard = flags[["negative", "above_limit", "repeated"]].any(axis=1) soft = flags[["outlier", "wet_spell", "dry_spell", "jump", "monthly_total"]].any(axis=1) flags["n_flags"] = flags[FLAG_COLS].sum(axis=1) flags["status"] = np.select([flags["missing"], hard, soft], ["missing", "fail", "suspect"], default="ok") out = pd.DataFrame({"prcp": v}) return pd.concat([out, flags], axis=1) ``` Wet spell for all-NaN: v>0 all False, _spell_lengths fine. One more check: `_spell_lengths(v == 0)` — v==0 for NaN gives False. Good. check_jump: ```python def check_jump(v, threshold): ok = v.notna() prev_ok = ok.shift(fill_value=False) return ok & prev_ok & (v.diff().abs() > threshold) ``` v.diff() at first element NaN → comparison False. Good. check_outliers: ```python def check_outliers(v, iqr_k=3.0, min_wet_days=20): flag = pd.Series(False, index=v.index) wet = v[v > 0].dropna() if wet.size < min_wet_days: return flag q1, q3 = wet.quantile([0.25, 0.75]) iqr = q3 - q1 if iqr > 0: upper = q3 + iqr_k * iqr flag = v > upper # if IQR == 0 (e.g. many identical small values), skip: too little spread for the test return flag ``` v > upper: NaN → False. Good. check_monthly_total: ```python def check_monthly_total(v, max_monthly): flag = pd.Series(False, index=v.index) if max_monthly is None: return flag per = v.index.to_period("M") totals = v.groupby(per).sum(min_count=1) bad = totals.index[totals > max_monthly] flag[:] = per.isin(bad) flag &= v.notna() return flag ``` `per.isin(bad)` — PeriodIndex.isin accepts container. Good. check_repeated final: ```python def check_repeated(v, max_repeated): wet = v > 0 if not wet.any(): return pd.Series(False, index=v.index) new_run = v.ne(v.shift()) | ~wet run_id = new_run.cumsum() run_len = v.groupby(run_id).transform("count") return wet & (run_len > max_repeated) ``` Wait: new_run for wet days: v.ne(v.shift()) — for NaN shift boundaries etc. For wet run, values equal → ne False → run continues. NaN inside breaks (ne True). Dry days (~wet True) increment. Good. transform("count") counts non-NaN in group; group includes only the run's wet days plus...? No — group includes all consecutive elements with same run_id: a run of wet equal values, terminated when ne True or dry. The group for run_id k consists exactly of that run (since any change increments). But note: dry days each have unique ids (each increments), except consecutive dry days: ~wet True each day → ne... day1 dry: new_run True → id increments. day2 dry: ~wet True → increments again. So each dry day own id. Good. NaN days: v.ne(v.shift()) True (NaN != anything) → increment. Also NaN != NaN True → consecutive NaNs each increment? v.ne(v.shift()): day i NaN, day i+1 NaN: NaN != NaN → True → increments. Yes each missing day own id. Good. Edge: wet run where a value is NaN in middle: run breaks; after NaN, next wet day: v.ne(v.shift()) → value != NaN → True → new run. Good. transform("count") on group of size n where all values equal non-NaN → n. flag wet & n > max_repeated → flags run of ≥ max_repeated+1 identical days. Hmm: "max_repeated=3" meaning allow up to 3 identical days; flag if 4+. In demo I inject 5 days of 2.0 → flagged. Note: should identical values like 0.0/0.1 traces be excluded? Some gauges round to 0.1 — repeated 0.1 for 4 days plausible? RClimDex flags any identical nonzero run ≥ 3 days? Their default: repeated value check for precipitation flags runs of ≥3 days? I recall RClimDex flags "sequence of repeated values" for all elements, with precipitation: runs of ≥ 3? For temperature, ≥ 3. For prcp, I think also ≥ 3 days with identical nonzero. Some argue small traces like 0.1 common. I'll keep threshold configurable and note. Now demo: ```python if __name__ == "__main__": rng = np.random.default_rng(7) idx = pd.date_range("2021-01-01", "2023-12-31", freq="D") n = len(idx) rain = np.zeros(n) wet = rng.random(n) < 0.35 rain[wet] = rng.gamma(0.6, 8.0, wet.sum()) rain = np.round(rain, 1) # inject typical errors rain[120] = -2.0 # negative rain[240] = 5000.0 # impossible value rain[300:305] = 2.0 # stuck sensor: 5 x same value rain[405] = 160.0 # jump 0 -> 160 -> 0 rain[600:636] = rng.gamma(0.7, 2.0, 36).round(1) # 36-day wet spell # long dry spell: ensure ~200 dry days rain[900:1100] = 0.0 s = pd.Series(rain, index=idx) s = s.drop(s.index[[1500, 1501, 2345]]) # missing dates (gaps) qc = run_qc(s) print(qc_report(qc)) ``` Wait: rain[600:636] gamma values: 36 consecutive wet days > max_wet_spell 30 → flagged. Also values vary so repeated not triggered (unless equal). Good. Dry spell 900:1100 zeros = 200 days > 180 → flagged dry_spell. But also random dry spells elsewhere might exceed 180? With p(wet)=0.35, dry run of 200 extremely unlikely. fine. Jump: rain[405]=160, neighbors 0 (or gamma values if wet). day 404 might be wet random with value maybe ~30; diff 130 < 150 maybe; or if rain[404]==0, jump 160>150 flagged at day 405, and day 406 diff -160 flagged at day 406 if rain[406]==0. Random wet could make diffs bigger. It's fine—report shows counts. Also above_limit 5000 flagged. negative flagged. repeated 5x2.0 flagged (and maybe random equal runs? gamma continuous rounded to 0.1 — chance of 4 identical consecutive small values like 0.0? wet mask ensures >0; two consecutive same rounded values possible but 4 in a row unlikely. However 0.1,0.1 etc across dry-wet boundary? Only wet>0 counted. OK.) Monthly total: default None (disabled) in DEFAULT_CFG. Good; mention enable for extreme regions. Actually maybe default 5000 mm enabled? Region-specific; keep None default and mention. Note also: drop indices 1500,1501,2345 create missing dates flagged as missing (gaps). Good. Sample output: I should present plausible output. I can't run code, so present an illustrative output with a note "(numbers will vary)". Let me compute approximate counts: - missing: 3 (gaps) - negative: 1 - above_limit: 1 - repeated: 5 (the stuck run) - outlier: depends on fence. wet values gamma(0.6, 8): scale heavy. Let's estimate quantiles: gamma shape k=0.6, theta=8. CDF: mean 4.8, var 38.4, sd 6.2. q75 maybe ~6.5? For shape 0.6, distribution is L-shaped (density ∞ at 0). Let me estimate: P(X ≤ x) for gamma(0.6): using gamma CDF. For shape a=0.6, theta=8. Median ~ theta * ~0.55? For a=0.6, median ≈ 0.63*? Hmm approximate: for small shape, median ≈ (a - 1/3)/a ... Let me just say q1≈1.2, q3≈6.4? IQR≈5.2; fence = q3 + 3*5.2 ≈ 22. Number of wet days exceeding 22 mm/day: P(X>22) for gamma(0.6,8): tail. Survival: Γ(0.6, 22/8=2.75)/Γ(0.6). Incomplete gamma upper at 2.75 with a=0.6... hard. Roughly maybe 3-5%? 0.35*1096 ≈ 384 wet days; 4% ≈ 15 days flagged as outliers. That's plausible "suspect" count but makes demo noisy. Maybe use iqr_k default 3 but in demo it's fine to show some outliers — realistic QC does flag some. But if fence yields many natural flags, example output lists "outlier: 15". Acceptable; I can note tuning. Alternatively set demo cfg iqr_k=5 → fence higher: q3 + 5*IQR ≈ 6.4+26=32.4; tail beyond 32.4: maybe ~1.5% ≈ 6 days. Hmm. I can't be precise; in the sample output I'll write approximate numbers and add "your numbers will vary slightly". Or avoid printing specific counts and instead show report format with plausible values and a disclaimer. I'll show a sample output block with values like: ``` Daily rainfall QC report ---------------------------------------- Records evaluated : 1096 negative : 1 e.g. 2021-05-01 above_limit : 1 e.g. 2021-08-28 repeated : 5 e.g. 2021-10-27, 2021-10-28, ... outlier : 6 e.g. 2021-03-14, ... wet_spell : 36 e.g. ... dry_spell : 200 ... jump : 2 ... status counts : {'ok': ..., 'suspect': ..., 'fail': ..., 'missing': 3} ``` Dates: index 120 from 2021-01-01 → day 120 is 2021-05-01 (Jan31+Feb28+Mar31+Apr30 = 120 days → index 120 = May 1, since index 0 = Jan 1). Yes: 2021-05-01. Index 240: 240 days after Jan 1 → Jan(31)+Feb(28)+Mar(31)+Apr(30)+May(31)+Jun(30)+Jul(31)+Aug(31)=243? Let's compute cumulative: day-of-year for index i = i+1. i=240 → DOY 241 → 2021: Aug 29 (DOY 241: Jan31,Feb59,Mar90,Apr120,May151,Jun181,Jul212,Aug243 → DOY 241 = Aug 29). So 2021-08-29. Index 300: DOY 301 = Oct 28 (Jan..Sep 273, +28 → Oct 28). run 300..304 → Oct 28–Nov 1. Index 405: DOY 406 → 2022? 2021 has 365 days, so index 365 = 2022-01-01. DOY 406-365=41 → Feb 10, 2022. Index 600: 600-365=235 → DOY 236 of 2022 = Aug 24, 2022 (2022: Jan31,Feb59,Mar90,Apr120,May151,Jun181,Jul212,Aug243; DOY236=Aug 24). Spell 600..635 = Aug 24 – Sep 28, 2022. Dry 900: 900-365=535 → DOY 536 of 2022... 2022 has 365 days; 900-730=170 → 2023? 730 days covers 2021+2022; index 730 = 2023-01-01. 900-730=170 → DOY 171 = Jun 20, 2023 (Jan31,Feb59,Mar90,Apr120,May151,Jun171 → DOY171=Jun20). Dry spell 900..1099 = Jun 20 – Dec 31? 1099-730=369 > 364 — beyond 2023 (n=1096 for 2021-2023: 365+365+365=1095; plus? 2021:365,2022:365,2023:365=1095). Wait date_range 2021-01-01 to 2023-12-31 daily = 1095 days. I said n=1096 earlier — it's 1095. Fine. But rain[900:1100] → 900..1094 (clipped) = 195 days. OK dry spell ~195 days. Rather than compute exact dates in sample output, I'll show generic sample output with a note that dates/counts vary with the seed. That's safer. I'll present output like: ``` Daily rainfall QC report -------------------------------------------------- Records evaluated : 1095 missing : 3 e.g. 2022-02-10, 2022-02-11, 2025-05-28 negative : 1 e.g. 2021-05-01 above_limit : 1 e.g. 2021-08-29 repeated : 5 e.g. 2021-10-28, 2021-10-29, 2021-10-30 ... outlier : 4 e.g. ... wet_spell : 36 ... dry_spell : 195 ... jump : 2 ... Status counts : ok=... suspect=... fail=... missing=3 ``` But my computed dates: index 1500 doesn't exist for n=1095! I set drop indices [1500,1501,2345] — out of range → error. Fix: drop e.g. [500, 501, 900]? But 900 is in dry spell; dropping inside affects dry spell length (195). Use drops [400, 401, 1000]: 400-730<0 → 400 is 2022 (DOY 401-365=36 → Feb 5, 2022); 1000-730=270 → DOY 271 = Sep 28, 2023. Dropping 400/401 adjacent — creates gap of 2 days; but note 405 jump remains. OK use `s.drop(s.index[[400, 401, 1000]])`. Careful: injecting jump at 405 unaffected. Hmm wait, also check: dry spell 900:1100 with drop at 1000 → dry run split into 100 (900-999) and 94 (1001-1094); neither > 180 → dry_spell not flagged! That ruins the demo of dry spell. Choose drop at [400, 401, 850]? 850-730=120 → May 1, 2023. Inside dry spell (900..) no—850 < 900 fine. So drop [400, 401, 850]. Then dry spell intact 195 days. But random natural data between... rain[850] random — dropping it just creates a gap elsewhere. Good: drops = [400, 401, 850]. But wait, dropping two adjacent days 400,401: gap of 2 missing dates → flagged missing (2 days). Plus 850 → 3 missing total. Good. Also jump demo at 405: I set rain[405]=160. If day 404 is wet random with say 20mm, diff = 140 < 150 → no flag; then I'd claim jump flagged in sample output incorrectly. To make deterministic: set rain[403:406] = [0, 0, 160]? Actually set rain[404]=0, rain[405]=160, rain[406]=0: diffs 160 and -160 → two jump flags (day 405 and 406), provided 403/407 not making bigger diffs (if 403 wet with 30mm, day404 diff -30 fine; day 407 diff 0 fine). Deterministic enough if I force 404 and 406 to 0 and 403, 407 also 0 to be safe: rain[403:408] = [0,0,160,0,0]. Then jumps: diff(404→405)=160>150 flag day405; diff(405→406)=-160 abs 160 flag day406; diff(406→407)=0. But also if 402 wet with e.g. 25 → diff to 403 is -25 fine. Good. Wet spell 600:636 with random gamma values: 36 consecutive wet days; but random wet elsewhere could extend run if rain[599] wet or rain[636] wet! rain[636] random: if wet, run becomes 37; flagged either way (36>30). Start boundary: rain[599] wet would extend start; still flagged. Fine. But also: is wet_spell count = run length (36 or 37)? The flag marks all days in run → count 36-37. In sample output say 36. But also natural long wet runs? p=0.35 → run of 31 extremely unlikely. OK. Outliers: I'll present with a small number, e.g., 5, and note it varies. Also note outlier fence computed including the injected 5000 (robust quantiles, minimal effect). Also natural large values (gamma tail) flagged — fine. Actually, wait: there's an issue with outliers: the check flags v > upper where upper from wet quantiles; with k=3 (extreme fence) and heavy-tailed gamma, natural maxima may exceed. That's the nature of the test; in practice you iterate/tune k. Fine. n_flags and status: fail = negative(1) + above_limit(1) + repeated(5) = 7 days fail (distinct). suspect = outlier + wet_spell(36) + dry_spell(195) + jump(2) overlapping? distinct days: wet spell 36, dry 195, jump 2, outliers maybe 4-8 distinct (some may coincide with... unlikely). So suspect ≈ 235-241. ok ≈ 1095 - 3 - 7 - ~237 ≈ 848. I'll write sample output with approximate numbers and note they vary. Hmm, presenting "sample output" with numbers I can't verify is a bit risky but standard; I'll add "(output will vary with the random seed)". Also maybe provide reading from CSV snippet: `pd.read_csv("rain.csv", parse_dates=["date"], index_col="date")["rain_mm"]`. Include in usage notes. Also optional: save flagged rows: `qc[qc.status != "ok"].to_csv("qc_flags.csv")`. Mention. Optional plotting: I'll include a compact function: ```python def plot_qc(qc, ax=None): import matplotlib.pyplot as plt ... colors = {"ok": "tab:blue", "suspect": "tab:orange", "fail": "tab:red", "missing": "lightgray"} ``` Maybe keep it short: ```python def plot_qc(qc): import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(12, 4)) for status, color in [("ok", "0.7"), ("suspect", "tab:orange"), ("fail", "tab:red")]: sub = qc[qc["status"] == status] ax.bar(sub.index, sub["prcp"], color=color, width=1.0, label=status) miss = qc[qc["status"] == "missing"] ax.bar(miss.index, 1, color="k", alpha=.3, width=1.0, label="missing") ax.set_ylabel("rainfall (mm)"); ax.legend() return fig, ax ``` Bar plot daily for 3 years fine. Maybe optional; include briefly. I think including plotting is a nice extra but adds length; I'll include it compactly since QC review typically wants visualization. Or skip to keep focus. The user asked to "write a quality control" — the core is the QC code. I'll include plotting as a short optional block at the end of the code (a few lines). Hmm. Let me include it — it's genuinely useful — but keep to ~10 lines. Now think about flag conventions: Some users want numeric codes (0=good,1=suspect,2=fail,-1=missing). I could add mapping comment. Keep boolean flags + status string; mention conversion. Also add docstring references: WMO Guidelines on Quality Control (WMO-No. 488), ETCCDI/RClimDex logic. Good to cite in comments briefly. Let me also double-check `_spell_lengths` correctness with a quick mental test: mask = [T,T,F,T]. shift fill False: [F,T,T,F]. new_run = mask & ~shift = [T, F, F, T]. run_id = [1,1,1,2]. groupby(run_id).transform("sum"): group1 (elements 0,1,2): mask sum = 2; group2 (element3): 1. → [2,2,2,1]; where(mask,0) → [2,2,0,1]. Correct. groupby on run_id where run_id=0 occurs for leading Falses: e.g., mask=[F,T]: shift=[F,F]; new_run=[F,T]; run_id=[0,1]; group0: mask sum 0 → where → 0; group1: 1. Good. One more consideration: `mask.groupby(run_id)` — if mask and run_id are pandas Series aligned; groupby by another series works via same index. transform("sum") on bool returns int64 (pandas may return bool? No, sum of bool → int). OK. For check_repeated, `v.groupby(run_id).transform("count")` returns counts; where group has NaNs? Only if NaN within same run_id — can't happen since NaN increments id. Fine. Potential performance: fine for typical station data. Now also handle case where series has no valid data / empty: guard `if s.empty: raise ValueError`. Add small guard. Also `s.index.min()`/max(): if s empty after dedup → error; guard. Let me also think about whether to use `min_count` in monthly sum: `v.groupby(per).sum(min_count=1)` — groupby.sum supports min_count. Yes (pandas ≥0.22? groupby sum with min_count supported). Good. FLAG_COLS constant list; ensure order used in n_flags and report. Write qc_report: ```python FLAG_COLS = ["missing", "negative", "above_limit", "repeated", "outlier", "wet_spell", "dry_spell", "jump", "monthly_total"] def qc_report(qc: pd.DataFrame, max_dates: int = 6) -> str: lines = ["Daily rainfall QC report", "-" * 46, f"records evaluated : {len(qc)}"] for col in FLAG_COLS: hits = qc.index[qc[col]] if len(hits): ex = ", ".join(d.strftime("%Y-%m-%d") for d in hits[:max_dates]) more = " ..." if len(hits) > max_dates else "" lines.append(f"{col:<15}: {len(hits):4d} e.g. {ex}{more}") counts = qc["status"].value_counts().reindex(["ok", "suspect", "fail", "missing"], fill_value=0) lines.append("-" * 46) lines.append("status : " + ", ".join(f"{k}={v}" for k, v in counts.items())) return "\n".join(lines) ``` Good. Now assemble the full script. Also add short docstring with references and note on spatial QC not included (needs neighbor stations) with a brief pointer of how (e.g., correlation/bias vs neighbor, or "spatial regression test"). I'll also add small note: for hourly data adapt thresholds. Let me also reconsider `check_outliers` min_wet_days param and degenerate IQR: add comment. One more useful check: "trace/accumulation" — many dry days then huge value suggests accumulator bucket emptied (common with tipping buckets that report accumulated rain). The jump check partially catches it. Mention in notes. Alright, format the answer: - Title line - "What it checks" bullet list (concise) - Code block (complete script with demo in __main__) - "Typical output" small block - "Adapting to your data / tuning" bullets - "Extensions" bullets (spatial QC, iterative cleaning, writing flags to file) Keep prose minimal. Use markdown headers sparingly (maybe just bold or ###). Code comments carry explanation. Double-check imports: numpy, pandas; matplotlib inside plot function. Type hint `dict | None` requires py3.10; use `cfg: dict = None` to be safe or `Optional[dict]`. Use `cfg: "dict | None" = None`? Simpler: `def run_qc(prcp, cfg=None):` without annotation issue. I'll use minimal annotations to avoid version issues but some hints fine (Python 3.9-friendly: use Optional). I'll just avoid `|` syntax. Write final code: ```python """ qc_rainfall.py — Quality control for daily rainfall (precipitation) time series. Checks (inspired by RClimDex/ETCCDI and WMO QC guidelines, WMO-No. 488): 1. missing : NaN / missing codes / gaps in the date calendar 2. negative : value < 0 -> hard error 3. above_limit : value > physical upper limit -> hard error 4. repeated : > N consecutive identical non-zero values -> stuck gauge 5. outlier : wet-day amounts outside IQR fence (boxplot) 6. wet_spell : too many consecutive wet days 7. dry_spell : too many consecutive dry days 8. jump : unrealistic day-to-day change 9. monthly_total : monthly total above plausible maximum (optional) Flag convention: True = failed the check. Final 'status': missing | fail (hard error) | suspect (needs inspection) | ok """ import numpy as np import pandas as pd DEFAULT_CFG = { "missing_codes": (-999.0, -99.0, -9999.0), "upper_limit": 1000.0, # mm/day; world record ~1825 mm (Foc-Foc, Réunion, 1966... ) "max_repeated": 3, # flag if MORE than this many identical non-zero days in a row "iqr_k": 3.0, # outlier fence = Q3 + k*IQR (wet days only); 3 ≈ extreme boxplot "min_wet_days": 20, # minimum wet days required before the outlier test is applied "max_wet_spell": 30, # consecutive wet days (value > 0) "max_dry_spell": 180, # consecutive dry days (value == 0) "jump_threshold": 150.0, # mm; |P(t) - P(t-1)| above this is suspicious "max_monthly_total": None,# e.g. 5000.0 to enable; region dependent } FLAG_COLS = ["missing", "negative", "above_limit", "repeated", "outlier", "wet_spell", "dry_spell", "jump", "monthly_total"] ``` Then helper functions and run_qc, qc_report, plot_qc, demo. Let me write demo values to look sensible. Also in demo, printing `qc.loc["2021-05-01"]` row? Could show. Maybe print first few flagged rows: `print(qc[qc.status != "ok"].head())`. Keep output section short. Now think again about the sample output numbers: with seed 7 values unknown; I'll write "example output (numbers will vary with the random seed)" and give plausible numbers. Fine. Edge: In demo, `rain[300:305] = 2.0` — 5 identical days. But are days 300-304 guaranteed not coinciding with natural wet random? We overwrite them → yes identical 2.0 exactly. But watch: rain[299] random could be 2.0 (rounded) → run 6; still flagged. Fine. Wet spell 600:636 overwritten with gamma values rounded — all > 0? gamma(0.7, 2) can produce very small values but > 0; rounding to 1 decimal could produce 0.0! E.g., gamma(0.7,2) with shape 0.7: values near 0 possible; round(0.04,1) = 0.0 → breaks wet run (v==0). That would break the 36-day run into pieces < 30! Risk. Fix: use `np.round(..., 1)` then clip: use gamma(0.8, 3) larger scale → less chance of <0.05? P(X < 0.05) for gamma(0.8,3): small but possible (~ shape 0.8: density near 0 diverges as x^{-0.2}; P(X<0.05) = γ(0.8, 0.0167)/Γ(0.8); γ(0.8,0.0167) ≈ (0.0167^0.8)/0.8 ≈ (0.0466)/0.8≈0.058? wait compute: x^a/a = 0.0167^0.8/0.8. 0.0167^0.8 = e^{0.8*ln0.0167}=e^{0.8*(-4.094)}=e^{-3.275}=0.0379; /0.8 = 0.0474. Γ(0.8)≈1.164. P≈0.0407. So ~4% chance per draw of rounding to 0.0; with 36 draws, ~77% chance at least one becomes 0.0 → breaks spell. Bad. Solution: don't round the wet spell values: `rain[600:636] = rng.gamma(0.8, 2.0, 36)` unrounded (floats) → all > 0 strictly. But then mixed with rounded others — fine, no issue. Or set spell values via max(rounded, 0.1): `np.maximum(np.round(rng.gamma(0.8, 2, 36), 1), 0.1)`. I'll use that to keep 1-decimal consistency. Also could produce values like exactly equal consecutive? No. Similarly main wet random: rounding to 0.0 means "no rain" effectively; acceptable (traces). But careful: rounding wet gamma to 0.1: values like 0.04 → 0.0 → becomes dry day. Fine. Also demo: ensure some natural outliers exist to show outlier check? Not necessary; the injected big values (5000) flagged above_limit, not outlier. Outlier may flag natural heavy days — fine, maybe zero in some seeds. In sample output I'll include outlier row with small count and caveat. Alternatively inject a suspicious-but-plausible value like 95 mm when climatological fence ~ 40 → flagged outlier. Deterministic: add `rain[700] = 95.0` — likely above fence (fence ≈ q3+3*IQR; wet gamma(0.6,8): q3? Let me estimate more carefully. Gamma(0.6, 8): CDF via Γ(a, x/θ). For a=0.6: P(X ≤ x) = γ(0.6, x/8)/Γ(0.6). Quantiles: median: solve γ(0.6, m/8)=0.5*Γ(0.6). Γ(0.6) ≈ 1.6446? Actually Γ(0.6) = 1.6446 (since Γ(1.6)=0.6*Γ(0.6)→Γ(1.6)≈0.8935, so Γ(0.6)≈1.489). Hmm: Γ(1.6) ≈ 0.8935; Γ(0.6) = Γ(1.6)/0.6 ≈ 1.4892. OK. γ(0.6, t) ≈ t^0.6/0.6 - t^1.6/1.6 + ... For small t. Median: γ(0.6, t)=0.7446. t^0.6/0.6 ≈ 0.7446 + t^1.6/1.6... try t=1.6: 1.6^0.6 = e^{0.6*0.470}=e^{0.282}=1.326; /0.6=2.21 too big. Try t=0.8: 0.8^0.6=e^{0.6*(-0.223)}=e^{-0.134}=0.875; /0.6=1.458; minus 0.8^1.6/1.6: 0.8^1.6=e^{1.6*(-0.223)}=e^{-0.357}=0.700; /1.6=0.437; γ≈1.02. Slightly above 0.7446. Try t=0.6: 0.6^0.6=e^{-0.3065}=0.736;/0.6=1.227; 0.6^1.6=e^{-0.817}=0.442;/1.6=0.276; γ≈0.951. Still >0.7446. t=0.45: 0.45^0.6=e^{0.6*(-0.7985)}=e^{-0.479}=0.619;/0.6=1.032; 0.45^1.6=e^{-1.278}=0.279;/1.6=0.174; γ≈0.858. t=0.35: 0.35^0.6=e^{0.6*(-1.0498)}=e^{-0.63}=0.533;/0.6=0.888; 0.35^1.6=e^{-1.68}=0.186;/1.6=0.117; γ≈0.771. t=0.32: 0.32^0.6=e^{0.6*(-1.139)}=e^{-0.683}=0.505;/0.6=0.842; 0.32^1.6=e^{-1.822}=0.162;/1.6=0.101; γ≈0.741≈0.7446 ✓. So median t≈0.32 → median = 0.32*8 = 2.56 mm. So median wet day ≈ 2.6 mm. q75: γ=0.75*1.4892=1.117. t: t=0.8 gave 1.02; t=0.95: 0.95^0.6=e^{0.6*(-0.0513)}=e^{-0.0308}=0.970;/0.6=1.616; 0.95^1.6≈0.922;/1.6=0.576; +next term +t^2.6/2.6: 0.95^2.6≈0.876/2.6=0.337 → γ≈1.616-0.576+0.337=1.377?? series alternating: γ(a,x)=Σ (-1)^{n} x^{a+n}/(n!(a+n))... terms: x^a/a - x^{a+1}/(1!(a+1)) + x^{a+2}/(2!(a+2))... For t=0.95: term0=0.970/0.6=1.616; term1=0.922/1.6=0.576; term2=0.95^2.6/(2*2.6)=0.876/5.2=0.168; term3=0.95^3.6/(6*3.6)=0.832/21.6=0.0385; γ≈1.616-0.576+0.168-0.0385≈1.170. Hmm that's > 1.117. t=0.85: 0.85^0.6=e^{0.6*(-0.1625)}=e^{-0.0975}=0.907;/0.6=1.512; 0.85^1.6=e^{-0.26}=0.771;/1.6=0.482; 0.85^2.6=e^{-0.4225}=0.655;/5.2=0.126; 0.85^3.6=e^{-0.585}=0.557;/21.6=0.026; γ≈1.512-0.482+0.126-0.026≈1.130. t=0.82: 0.82^0.6: ln0.82=-0.1985; *0.6=-0.119; e=0.888;/0.6=1.480; 0.82^1.6: ln*1.6=-0.3176; e=0.728;/1.6=0.455; 0.82^2.6: -0.516; e=0.597;/5.2=0.115; 0.82^3.6: -0.7146; e=0.489;/21.6=0.0226; γ≈1.480-0.455+0.115-0.023=1.117 ✓. So q75 t≈0.82 → 6.56 mm. IQR = 6.56 - 2.56 = 4.0. Fence = 6.56 + 3*4.0 = 18.56 mm. Natural days > 18.56: P(X>18.56) = Γ_upper(0.6, 2.32)/Γ(0.6). Upper incomplete at z=2.32 with a=0.6: approximate via 1 - γ(0.6, 2.32). γ(0.6,2.32): terms: 2.32^0.6=e^{0.6*0.8416}=e^{0.505}=1.657;/0.6=2.762; 2.32^1.6=e^{1.3466}=3.844;/1.6=2.403; 2.32^2.6=e^{2.188}=8.916;/5.2=1.715; 2.32^3.6=e^{3.0296}=20.7;/21.6=0.958; 2.32^4.6=e^{3.871}=48.0;/(24*4.6=110.4)=0.435; 2.32^5.6=e^{4.713}=111.3;/(120*5.6=672)=0.166; 2.32^6.6=e^{5.554}=258.5;/(720*6.6=4752)=0.0544; 2.32^7.6=e^{6.396}=599;/(5040*7.6=38304)=0.0156; γ≈2.762-2.403+1.715-0.958+0.435-0.166+0.054-0.016≈1.423; next terms: -2.32^8.6/(40320*8.6): 2.32^8.6=e^{7.237}=1392; /(346752)=0.004 → ≈1.419; +0.001 → ~1.420. 1-1.420/1.489=0.0463. So ~4.6% of wet days exceed 18.6mm → with ~384 wet days ≈ 18 natural "outlier" flags + maybe injected. Hmm that's a lot for a demo but it's the nature of a heavy-tailed daily rainfall with k=3 IQR fence. Real QC on real data: daily rainfall is heavy-tailed; IQR fence with k=3 flags ~5% of wet days — too many. Practical approaches: use k larger, or use percentile thresholds (e.g., flag > 99.9th pct of wet days), or RClimDex uses k chosen by user (default often 1 or 2 with manual inspection; they use boxplot and user picks). Hmm. RClimDex default outlier test for PRCP uses "boxplot" and flags values outside Q3 + k*IQR where user typically sets k... Actually RClimDex's precipitation outlier uses quantile 0.05/0.95 based? Let me recall: RClimDex step: "precipitation outlier test: values outside [Qn(0.05)-3*IQR...]"? I don't remember exactly. There's the approach in ETCCDI software: uses (median ± ... )? To be safe, I'll implement IQR fence but ALSO offer percentile fallback and note that for skewed rainfall, better to set k high (e.g., 10) or use a percentile cap. For demo, set iqr_k default maybe 5? P(X>q3+5*4=26.56): γ(0.6,3.32): terms: 3.32^0.6=e^{0.6*1.2}=e^{0.72}=2.054;/0.6=3.424; 3.32^1.6=e^{1.92}=6.82;/1.6=4.263; 3.32^2.6=e^{3.12}=22.6;/5.2=4.35; 3.32^3.6=e^{4.32}=75.2;/21.6=3.48; 3.32^4.6=e^{5.52}=249;/110.4=2.26; 3.32^5.6=e^{6.72}=828;/672=1.232; 3.32^6.6=e^{9.12}=9147.../4752=1.92; hmm series diverging terms growing — asymptotic series not converging at z=3.32 with a=0.6; the alternating series still needs care. This is getting too deep. Alternative: make demo robust by setting cfg in demo: `cfg = {"iqr_k": 6, "max_dry_spell": 180, ...}`. With k=6, fence = 6.56+24=30.6 mm; tail beyond 30.6: maybe ~2.5%? Still some natural flags. Honestly, for the demo it's fine to flag several natural heavy days — QC reports typically do; I'll present sample output with e.g. "outlier: 7" and note it depends on seed/k and that flagged extremes may be genuine (verify with nearby stations). Simpler: in demo, don't claim exact counts; give output format with representative numbers and a comment. OK. Alternatively, skip printing a fake output and instead describe what the report contains. But showing example output helps. I'll include a sample output with plausible numbers, clearly noting variability. Fine. Also note: rounding wet gamma(0.6,8) to 1 decimal fine. Now consider `check_outliers` uses quantiles of wet days including flagged errors (5000 at one day) — negligible for quantiles. Mention iterative option. Now let me also think about whether to flag "above_limit" at 1000 mm default. World record daily: 1825 mm (Foc-Foc, Réunion, Jan 1966? Actually 1825 mm in 24h at Foc-Foc La Réunion during Cyclone Denise 1966? Hmm: The 24-hour world record is 1825 mm at Foc-Foc, Réunion, 7–8 Jan 1966 (some sources). Cherrapunji 26,461 mm monthly July 1861; yearly record 26,470 mm Cherrapunji 1860-61. So upper_limit default 1200? Many national QC use country-specific (e.g., 500 mm for UK). I'll default 1000 and note regional tuning; mention records in comment. Fine. Also monthly default None; mention example values (e.g., 4000–10000 mm in extreme monsoon regions). Dry spell default 180 days: in arid regions, dry spells can exceed a year — note to tune (e.g., Atacama). Fine. Now finalize code text. Also add function `apply_flags_summary`? Not needed. One more check: in run_qc, after reindex, we compute flags on full calendar; the 'missing' flag includes originally-NaN and gap days. Good. Report includes them. Also note: `v = s.replace(list(cfg["missing_codes"]), np.nan)` — s may have ints? after reindex float. fine. Also careful: `flags["negative"] = v < 0` — NaN < 0 → False. Good. `flags["above_limit"] = v > cfg["upper_limit"]` — NaN → False. Good. Now assemble plot function: ```python def plot_qc(qc, title="Daily rainfall QC"): import matplotlib.pyplot as plt colors = {"ok": "#9ecae1", "suspect": "#f4a261", "fail": "#d64545", "missing": "#bdbdbd"} fig, ax = plt.subplots(figsize=(13, 4)) for status, color in colors.items(): sub = qc.loc[qc["status"] == status] if status == "missing": ax.bar(sub.index, np.ones(len(sub)), color=color, width=1.0, label=status) else: ax.bar(sub.index, sub["prcp"], color=color, width=1.0, label=status) ax.set_ylabel("rain (mm)") ax.set_title(title) ax.legend(loc="upper right", ncols=4) fig.tight_layout() return fig, ax ``` `ncols` param exists in newer matplotlib legend (1.8+?) as `ncols` (alias ncols added 3.6). Use `ncol=4` for compatibility. Use ncol. Also ax.bar with datetime index works. width=1.0 in days units fine. Now the demo `__main__`: ```python if __name__ == "__main__": rng = np.random.default_rng(7) idx = pd.date_range("2021-01-01", "2023-12-31", freq="D") n = len(idx) rain = np.zeros(n) wet = rng.random(n) < 0.35 rain[wet] = rng.gamma(0.6, 8.0, wet.sum()) rain = np.round(rain, 1) # --- inject typical real-world problems ------------------------------- rain[120] = -2.0 # negative value rain[240] = 5000.0 # impossible magnitude rain[300:305] = 2.0 # stuck sensor (5 identical days) rain[403:408] = [0.0, 0.0, 160.0, 0.0, 0.0] # suspicious jump rain[600:636] = np.maximum(np.round(rng.gamma(0.8, 2.0, 36), 1), 0.1) # 36-day wet spell rain[900:1100] = 0.0 # ~200-day dry spell s = pd.Series(rain, index=idx, name="rain_mm") s = s.drop(s.index[[400, 401, 850]]) # missing dates (gaps) qc = run_qc(s) print(qc_report(qc)) # keep the flagged rows for inspection / station manager feedback qc.loc[qc["status"] != "ok"].to_csv("rainfall_qc_flags.csv") ``` Wait: rain[900:1100] — n=1095 so 900:1100 clipped to 900:1095 → 195 days zeros. But also random wet may exist at 899/… boundaries: rain[899] random; run begins at 900 anyway if rain[899] wet? If rain[899] wet, run starts earlier but ends at 1094 → length ≥ 195 > 180 flagged. Good. Note: also 1094 is last day. fine. Also drop indices 400,401: gap 2 days. But careful: rain[400], rain[401] removed; jump injection at 403-407 unaffected. Good. drop 850 → gap. But wait: does dropping 400,401 affect dry spell demo? No (dry at 900+). Also check repeated demo: 5 days at 2.0 → run_len 5 > 3 → flagged (5 days). Also could natural consecutive equal rounded values occur elsewhere? Rounded to 0.1 with continuous gamma — two consecutive equal values possible (e.g., 0.5,0.5) but 4+ consecutive equal extremely unlikely. OK. Also: does the 36-day wet spell get also flagged as "repeated"? Values gamma rounded — unlikely 4 identical consecutive. fine. Jump: [0,0,160,0,0] at 403..407: diffs: 404→405: 160 → flag day405; 405→406: -160 → flag day406; also 402→403: rain[402] random wet maybe 20 → diff -20 fine. 407→408: rain[408] random; if wet 30 → diff 30. fine. But threshold 150: exactly 160 > 150 ✓. Now sample output: I'll write something like: ``` Daily rainfall QC report ---------------------------------------------- records evaluated : 1095 missing : 3 e.g. 2022-02-05, 2022-02-06, 2023-05-01 negative : 1 e.g. 2021-05-01 above_limit : 1 e.g. 2021-08-29 repeated : 5 e.g. 2021-10-28, 2021-10-29, 2021-10-30 ... wet_spell : 36 e.g. 2022-08-24, ... dry_spell : 195 e.g. ... jump : 2 e.g. ... ---------------------------------------------- status: ok=..., suspect=..., fail=..., missing=3 ``` Check dates: index 120 → DOY 121 → 2021: Jan31 Feb59 Mar90 Apr120 → DOY121 = May 1 ✓. index 240 → DOY 241 → Aug 29 (Jan31+Feb28=59, Mar90, Apr120, May151, Jun181, Jul212, Aug243; DOY241=Aug29 ✓). index 300..304 → DOY 301..305 → Oct 28..Nov 1 (Jan..Sep=273; 274=Oct1 → 301=Oct28 ✓; 305=Nov1 ✓). index 600 → 2022? 2021 has 365 days: idx 0..364 = 2021. idx 600 → 2022 DOY 236 (600-365=235; DOY=236) → Aug 24, 2022 ✓ (2022 non-leap: cumulative Jan31 Feb59 Mar90 Apr120 May151 Jun181 Jul212 Aug243; DOY236 = Aug 24 ✓). Dry spell idx 900.. → 2022 DOY 536-365=171 → DOY 171 = Jun 20, 2022 through 2022 end DOY365=Dec31, that's 195 days (900..1094). ✓. Gaps: idx 400 → 2022 DOY 36 → Feb 5, 2022 ✓ (DOY36: Jan31+5=Feb5 ✓). idx 401 → Feb 6 ✓. idx 850 → 2023? 850-730=120 → DOY121 2023 = May 1, 2023 ✓ (2023 non-leap). Great, dates line up: missing: 2022-02-05, 2022-02-06, 2023-05-01. negative: 2021-05-01. above_limit: 2021-08-29. repeated: 2021-10-28..Nov 1 (5 dates). wet_spell: 2022-08-24 onwards (36 dates: idx600..635 → DOY236..271 2022 → Aug24..Sep28; DOY271: Aug243+28=271 → Sep 28 ✓). dry_spell 195 dates Jun20–Dec31 2022. jump: 2022-02-09, 2022-02-10 (idx405,406 → DOY 406-365=41 → Feb10; 406-365=41 → wait idx405 → DOY406 → Feb 10 (DOY 406-365=41 → Feb10); idx406 → DOY 407 → Feb 11). Hmm: idx 403 → DOY 404 → Feb 8 (404-365=39 → Feb 8). idx405 → 41 → Feb 10; idx406 → 42 → Feb 11. So jump dates: 2022-02-10 and 2022-02-11. But do gaps at 400/401 affect jump? diff across gap: rain[399] random, rain[402] random (since 400,401 dropped from series but reindex restores them as NaN; v=NaN at 400,401; jump requires both valid → skipped at 402 (prev NaN) — good. Also wet_spell flag count: 36 days (run > 30 → whole run flagged: 36 days). dry_spell: 195 days flagged (whole run). That's how I coded (flag entire run). Some prefer flag only days beyond threshold; note it. status counts: total 1095; missing 3; fail = negative(1)+above_limit(1)+repeated(5) = 7 (no overlaps: idx120,240,300-304 distinct; note idx240 5000 also outlier? outlier is soft; status takes hard first via np.select order [missing, hard, soft] → fail. above_limit day also > fence → outlier True but status fail. Fine.) suspect = wet_spell 36 + dry_spell 195 + jump 2 + outlier (natural, count unknown) minus overlaps: idx600-635 wet spell values gamma(0.8,2) small — any above fence 18.6? gamma(0.8,2): mean 1.6; max of 36 draws maybe ~8-12; unlikely above fence → no overlap. Jump days 405/406 values 160/0: 160 > fence → outlier True too; but status: soft overlaps fine count once. So suspect ≈ 36+195+2+outliers(~18?) - overlap(2 jump days flagged outlier too: 160 flagged outlier, 0 no) ≈ 231 + outliers(≈18) minus 1 overlap ≈ 248? Eh, don't compute precisely; in sample output I'll put plausible numbers like "suspect≈250". Hmm, that's a lot of suspect which might look odd but it's realistic for the demo (mostly the long spells which are often legitimate). Maybe better to note in the demo that long spells are often real and QC flags them for inspection only. Alternatively reduce noise: set dry_spell threshold default 365? Then demo 195-day dry spell not flagged unless configured. Hmm. But then demo shows no dry_spell flags. I could demo dry spell by configuring cfg = {"max_dry_spell": 150}. Let me instead in demo pass cfg with max_dry_spell=150 to trigger on the 195-day spell. And keep default 365? What's a sensible default? For QC purposes, flagging >180 days dry is reasonable in humid climates; in arid, no. Default 180 with note. In demo, I'll pass cfg = {"max_dry_spell": 150}? Then 195 flagged. But also natural dry runs could exceed 150? p(dry)=0.65 per day; P(run>150) astronomically small. OK. But writing cfg overrides in demo adds noise; alternatively keep default 180 → the 195-day run flagged with default. Natural runs >180 essentially impossible. Keep default 180, no override. Good — simpler demo. Similarly wet spell default 30 → natural runs >30 impossible with p=0.35. Good. Outlier natural flags ~18 with k=3 — demo output shows e.g. 15-20; fine, I'll mention "extreme wet-day amounts above the IQR fence (verify against nearby stations / may be genuine storms)". Let me reduce uncertainty: in sample output, use "≈" or just show format with example numbers and a note "counts vary with the random seed". I'll show specific plausible numbers but explicitly say illustrative. Actually, maybe simpler: I can present the sample output without exact counts for outlier ("outlier: 14 e.g. ...") and note variability. Good enough. Now, one more design question: should `run_qc` accept a DataFrame with a column instead of Series? Accept Series; user can pass df["prcp"]. Mention. Also mention reading from CSV: `s = pd.read_csv("station.csv", parse_dates=["date"], index_col="date")["rain_mm"]`. Alright — also the "spatial QC" extension bullet: if neighbor stations available: compare same-day values (e.g., flag if station value large while all neighbors dry), or compute daily correlations; suggest tools (GSIM/ GHCN, or `pyhomogenize`, `climquality`?). Keep generic. Also mention flag convention alternative: WMO/GTS codes or 0/1/-1; and keeping original data untouched + adding QC columns. Also mention units: ensure mm/day. Now, code correctness re-check on `_spell_lengths` with pandas boolean groupby transform sum: `mask.groupby(run_id).transform("sum")` — mask is boolean Series; groupby sum on bool returns int (bool sum → int64 in groupby? Actually SeriesGroupBy.sum on bool dtype returns int64). Then `.where(mask, 0)` produces int with 0s. Comparison `> max` fine. But if mask empty (len 0)? Series empty; groupby transform on empty returns empty. Edge fine. But careful: `mask.shift(fill_value=False)` — for pandas 2.x, Series.shift(fill_value=...) supported. For older, `mask.shift().fillna(False)`. Use `.shift(fill_value=False)` — supported since pandas 0.24? shift fill_value added in 0.24. fine. In `check_repeated`: ```python def check_repeated(v, max_repeated): wet = v > 0 if not wet.any(): return pd.Series(False, index=v.index) run_id = (v.ne(v.shift()) | ~wet).cumsum() run_len = v.groupby(run_id).transform("count") return wet & (run_len > max_repeated) ``` v.groupby(run_id) — groupby by run_id (Series aligned). transform("count") counts non-NaN in each group. Within a wet run of identical values, all non-NaN → count = run length. But careful: group includes the trailing non-wet days? run_id increments on ~wet, so trailing dry days have new id. The group for a wet run = exactly the run (since ~wet increments). But NaN missing inside run: v.ne(v.shift()) — NaN != previous (whether previous NaN or value) → True → increments; so missing days break runs and get own ids? For missing day m: ne True → id increments; then next day after missing: v.ne(v.shift()) → value != NaN → True → increments again. So each missing day is its own run. Wet run after gap starts anew. But one subtle bug: a wet day followed by another wet day with the same value continues; with different value increments. Good. Also `(v.ne(v.shift()) | ~wet)` — for wet run: ne False (equal values) and wet True → False → same id. For dry days: ~wet True → increment each day (each dry day own run; count=1 ≤ max fine, but dry runs with equal 0 values not flagged since we require wet). Also note: v.ne(v.shift()) when previous is NaN and current is NaN → NaN != NaN → True → increment. Good. `transform("count")` on group where all NaN → 0. fine. Edge: `v.groupby(run_id)` — groupby keys from run_id Series; fine. Now `_spell_lengths` used for wet/dry spells: ```python def _spell_lengths(mask): new_run = mask & ~mask.shift(fill_value=False) run_id = new_run.cumsum() return mask.groupby(run_id).transform("sum").where(mask, 0) ``` Hmm groupby(run_id) where run_id for leading Falses is 0 and mask False → group 0 sum 0 → where → 0. For a run: ids contiguous; transform sum = run length. But subtle: two separate runs with different ids fine. Trailing Falses share last run id → sum = last run length → set 0 by where. OK. Wait issue: if mask starts False then first True run has id 1; leading False group id 0. If mask starts True, id 1. fine. But: `mask.groupby(run_id).transform("sum")` with bool dtype — returns? In modern pandas, bool series groupby sum → int64. `.where(mask, 0)` → int64 with 0 at False. Then compare `> max`. Good. Alternatively simpler: `mask.groupby(run_id).transform("size")`? size counts all rows in group including trailing falses (since they share the id) → wrong for spell length at trailing? transform size would give group size = run length + number of trailing falses → then where(mask,0) zeroes falses, but trues get inflated size! E.g., run [T,T] followed by 5 F: group size 7 → trues get 7 → false positives if 5+... So must use "sum" of the boolean (counts Trues only). Yes "sum" is correct: sum of True over group = run length regardless of trailing Falses. But careful: sum of bool in groupby — pandas may warn about deprecation? bool sum works. OK. Hmm, actually one more subtlety: bool Series groupby.sum() may raise in some versions? I believe it works (casts to int). Fine. Monthly check code as above. Now `np.select` for status: choicelist must be arrays; conditions arrays. `np.select([flags["missing"], hard, soft], ["missing","fail","suspect"], "ok")` — np.select requires condlist of arrays, choicelist arrays or scalars? np.select allows scalars in choicelist? Yes, choicelist elements can be scalars (broadcast). Actually np.select(conds, choices, default) — choices can be scalars; it broadcasts. Yes. hard = flags[["negative","above_limit","repeated"]].any(axis=1) — but also maybe include... keep these 3 hard. Also make sure flags columns exist: monthly_total included. Now write everything. Additionally maybe provide `qc` column naming and mention that you can save: `qc.to_csv("qc_output.csv")`. Potential improvement: Also flag "trace" values if coded as -1 → mention replacing traces. I think also worth including a brief "checks and typical thresholds" table in prose before code. Keep it compact: bullet list of checks with defaults. Length: moderate; code ~150 lines. The user asked to "write a quality control of daily rainfall values in python" — a script is the core deliverable. I'll include headers, code block, usage notes. Markdown appropriate (code, headers, bullets). Let me write final code now, carefully. ```python """ qc_rainfall.py — Quality control for daily rainfall (precipitation) time series. Checks ------ 1. missing : NaN, missing codes (-999, -99, ...) or dates absent from the index 2. negative : value < 0 (hard error) 3. above_limit : value > physical maximum (mm/day) (hard error) 4. repeated : >= max_repeated+1 consecutive identical non-zero values -> stuck/delayed tipping bucket (hard error) 5. outlier : wet-day amount above Q3 + k*IQR of the wet days (suspect) 6. wet_spell : run of wet days longer than max_wet_spell (suspect) 7. dry_spell : run of dry days longer than max_dry_spell (suspect) 8. jump : |P(t) - P(t-1)| > jump_threshold (both days valid)(suspect) 9. monthly_total : calendar-month total above max_monthly_total (suspect, optional) Status: 'ok' | 'suspect' (needs verification) | 'fail' (bad value) | 'missing' """ from __future__ import annotations import numpy as np import pandas as pd DEFAULT_CFG = { "missing_codes": (-999.0, -99.0, -9999.0), "upper_limit": 1000.0, # mm/day; world records: 1825 mm (Foc-Foc, Réunion 1966) "max_repeated": 3, # flag if a run of equal non-zero values exceeds this "iqr_k": 3.0, # outlier fence = Q3 + k * IQR (raise for skewed climates) "min_wet_days": 100, # minimum wet days before the outlier test is applied "max_wet_spell": 30, # consecutive wet days (rain > 0) "max_dry_spell": 180, # consecutive dry days (rain == 0); tune for arid climates "jump_threshold": 150.0, # mm change between two consecutive valid days "max_monthly_total": None, # e.g. 5000 mm; None disables the test } FLAG_COLS = ["missing", "negative", "above_limit", "repeated", "outlier", "wet_spell", "dry_spell", "jump", "monthly_total"] HARD = ["negative", "above_limit", "repeated"] # clear errors SOFT = ["outlier", "wet_spell", "dry_spell", "jump", "monthly_total"] # verify manually ``` Then helpers and run_qc and report and plot and main. Write `_spell_lengths`, `check_repeated`, `check_outliers`, `check_monthly_total`, `run_qc`, `qc_report`, `plot_qc`, `__main__`. Full code: ```python def _spell_lengths(mask: pd.Series) -> pd.Series: """For each day, the length of the run of consecutive True containing it (0 if False).""" run_id = (mask & ~mask.shift(fill_value=False)).cumsum() return mask.groupby(run_id).transform("sum").where(mask, 0) def _flag_repeated(v: pd.Series, max_repeated: int) -> pd.Series: """Runs of identical non-zero values longer than `max_repeated` (stuck gauge).""" wet = v > 0 run_id = (v.ne(v.shift()) | ~wet).cumsum() # value change or dry/missing day ends a run run_len = v.groupby(run_id).transform("count") # count of non-missing days in the run return wet & (run_len > max_repeated) def _flag_outliers(v: pd.Series, k: float, min_wet_days: int) -> pd.Series: """Boxplot ('extreme outlier') test applied to wet days only.""" flag = pd.Series(False, index=v.index) wet = v[v > 0].dropna() if wet.size < min_wet_days: return flag q1, q3 = wet.quantile([0.25, 0.75]) iqr = q3 - q1 if iqr == 0: # no spread (e.g. all wet days equal): skip return flag return v > q3 + k * iqr def _flag_monthly_totals(v: pd.Series, max_monthly: float | None) -> pd.Series: """All days of any month whose total exceeds `max_monthly`.""" flag = pd.Series(False, index=v.index) if max_monthly is None: return flag totals = v.groupby(v.index.to_period("M")).sum(min_count=1) bad = totals.index[totals > max_monthly] return pd.Series(v.index.to_period("M").isin(bad), index=v.index) & v.notna() ``` run_qc: ```python def run_qc(prcp: pd.Series, cfg: dict | None = None) -> pd.DataFrame: """Run all QC checks on a daily precipitation series (mm/day, datetime index).""" cfg = {**DEFAULT_CFG, **(cfg or {})} s = prcp.copy() s.index = pd.to_datetime(s.index) s = s.sort_index() if s.index.has_duplicates: # keep first occurrence of duplicated dates s = s[~s.index.duplicated(keep="first")] # insert calendar days missing from the index so gaps are caught by the missing test s = s.reindex(pd.date_range(s.index.min(), s.index.max(), freq="D")) v = s.astype(float).replace(list(cfg["missing_codes"]), np.nan) flags = pd.DataFrame(index=v.index) flags["missing"] = v.isna() flags["negative"] = v < 0 flags["above_limit"] = v > cfg["upper_limit"] flags["repeated"] = _flag_repeated(v, cfg["max_repeated"]) flags["outlier"] = _flag_outliers(v, cfg["iqr_k"], cfg["min_wet_days"]) flags["wet_spell"] = (v > 0) & (_spell_lengths(v > 0) > cfg["max_wet_spell"]) flags["dry_spell"] = (v == 0) & (_spell_lengths(v == 0) > cfg["max_dry_spell"]) ok_pair = v.notna() & v.notna().shift(fill_value=False) # both days valid flags["jump"] = ok_pair & (v.diff().abs() > cfg["jump_threshold"]) flags["monthly_total"] = _flag_monthly_totals(v, cfg["max_monthly_total"]) out = pd.DataFrame({"prcp": v, "status": "ok", index=v.index) # fix syntax ``` Careful with pd.DataFrame construction syntax: `pd.DataFrame({"prcp": v})` then assign. ```python out = pd.DataFrame({"prcp": v}) out[FLAG_COLS] = flags[FLAG_COLS] hard = flags[HARD].any(axis=1) soft = flags[SOFT].any(axis=1) out["status"] = np.select( [flags["missing"], hard, soft], ["missing", "fail", "suspect"], default="ok") return out ``` Order of np.select conditions: missing first; then hard; then soft. A day both hard and soft → "fail". Good. Note: `flags["dry_spell"] = (v == 0) & ...` — NaN == 0 → False. good. `ok_pair & (v.diff().abs() > threshold)` — diff at first NaN → comparison False anyway. ok_pair needed to avoid flagging across a missing day: diff between t and t-1 where t-1 NaN → NaN > thr → False. So ok_pair redundant? v.diff() across missing: if day t-1 missing, diff(t) = v(t) - NaN = NaN → False. If day t missing, NaN. So jump never flags across missing automatically. But if t-1 is flagged negative etc., still computes diff — fine (soft). So ok_pair unnecessary but harmless; drop it to simplify? Keep it minimal: `flags["jump"] = v.diff().abs() > cfg["jump_threshold"]` — with NaN → False. Simpler. But careful: v.diff() on reindexed series: gaps are NaN → no jump across gaps. I'll simplify: `flags["jump"] = v.diff().abs() > cfg["jump_threshold"]`. Report: ```python def qc_report(qc: pd.DataFrame, max_dates: int = 6) -> str: lines = ["Daily rainfall QC report", "-" * 60, f"records: {len(qc)} " f"ok: {(qc.status == 'ok').sum()} suspect: {(qc.status == 'suspect').sum()} " f"fail: {(qc.status == 'fail').sum()} missing: {(qc.status == 'missing').sum()}"] for col in FLAG_COLS: hits = qc.index[qc[col]] if len(hits): ex = ", ".join(hits[:max_dates].strftime("%Y-%m-%d")) more = " ..." if len(hits) > max_dates else "" lines.append(f" {col:<14}{len(hits):5d} e.g. {ex}{more}") return "\n".join(lines) ``` hits[:max_dates] — DatetimeIndex slice ok; .strftime works on DatetimeIndex. Main demo as designed. Also in `__main__`, print first lines of qc? Maybe show `qc.loc["2021-10-27":"2021-11-02"]`? Could print sample. Keep: print report; optionally save CSV. Also show how a flagged value appears: print a few rows. I'll add: ```python print("\nSample of flagged records:") print(qc.loc[qc["status"] != "ok", ["prcp", "status"]].head(12)) ``` Hmm head(12) shows missing+... fine. Double-check demo injection: rain[600:636] values gamma(0.8, 2.0) — mean 1.6, mostly small; rounding to 1 decimal; some could equal consecutive (e.g., 0.3,0.3) but 4 consecutive equal unlikely. Fine. But note: rain[600:636] overwritten after wet-random; ensures 36 wet days regardless. Also ensure boundaries: rain[599] random (maybe wet) — run could extend to 37 → still flagged, count 37. In sample output I say 36; note variability. Actually if rain[599] is wet (35% chance), run becomes 37 (if rain[598] dry...). Also if rain[636] wet → extends. So wet_spell count 36±2. In sample output, I'll put 36 and note illustrative. Or force boundaries dry: set rain[599] = 0 and rain[636] = 0 before? rain[636] random; force: after injection, `rain[599] = 0.0; rain[636] = 0.0` to pin the spell at exactly 36. Similarly dry spell: rain[899]=? it's random; force `rain[899] = wet-ish`? To pin dry spell exactly at 195: set rain[899] to a wet value e.g. 5.0, and rain[1094] is last day (no boundary). Let me do: ```python rain[600:636] = np.maximum(np.round(rng.gamma(0.8, 2.0, 36), 1), 0.1) rain[599] = 0.0; rain[636] = 0.0 # pin the wet spell to exactly 36 days rain[900:1100] = 0.0 rain[899] = 5.0 # pin the dry spell to exactly 195 days ``` Wait rain[900:1100] clipped to 1095 → 900..1094 = 195 days; boundary before: rain[899]=5.0 wet → run starts at 900 exactly. Good. And no trailing boundary. But rain[636]=0.0 — also ensures jump? diff 635→636 = v(635) small → 0: fine. Also ensure rain[402] maybe wet — jump day 405 needs |160 - rain[404]|; rain[404]=0 (injected) → 160 ✓; day 406: |0 - 160| = 160 ✓; day 407: |0 - rain[408]|: rain[408] random; if wet maybe 10-30 → no flag. Also day 403: |0 - rain[402]| small. So jump count = 2 exactly (unless rain[408] huge > 150 — probability tiny: P(gamma(0.6,8) > 150) ~ negligible; also rain[402] could be wet but |0 - x| < 150 almost surely). But also natural jumps: two consecutive wet days where one ~40mm and prev ~0? diff = 40 < 150. Wet-to-wet gamma values rarely differ by >150. P(gamma > 150) tiny. But wait, the injected 5000 at idx 240: diffs around it: day 240 diff |5000 - rain[239]| = ~5000 → jump True! And day 241: |rain[241] - 5000| ≈ 5000 → jump True. So jump count = 4 (2 natural around the 5000 spike + 2 injected). Also negative value day 120: diffs small (~ -2 - x). no. So jump: 4 days: idx 240, 241 (around 5000), 405, 406. Dates: idx240 → 2021-08-29; idx241 → 2021-08-30; idx405,406 → 2022-02-10, 2022-02-11. In sample output: "jump: 4 e.g. 2021-08-29, 2021-08-30, 2022-02-10 ...". OK, plausible. Similarly outlier: 5000 > fence → outlier True (soft; status fail due hard). Wet spell days small values below fence. Dry spell zeros not wet. Outlier count = natural tail days + idx240 (which is also above_limit & jump) ≈ natural ~10-20. I'll write ~"outlier 14" illustrative. Wait, also the negative value: v=-2 → wet? -2 > 0 False. outlier uses wet>0. fine. repeated: -2 not >0. fine. above_limit idx240=5000 flagged. Also is 5000 > q3 + 3*IQR? yes definitely → outlier True. In hard status takes precedence. Now compute expected "suspect" count roughly: wet_spell 36 + dry_spell 195 + jump 4 (two of which—idx240,241—are hard? idx240 hard (above_limit) → status fail; idx241: only jump (+outlier? |rain[241]-5000|? outlier is about value magnitude not diff; rain[241] random ~ small → not outlier) → suspect) + outlier natural (~14, one of which idx240 hard; idx241 not outlier) → suspect ≈ 36+195+ (jump: idx241,405,406 =3; idx240 is fail) + outlier(14 incl idx240? idx240 in outlier but status fail; so outlier contributes 13 suspect) ≈ 36+195+3+13 = 247. fail = 7 (negative 1, above_limit 1 (idx240), repeated 5). missing 3. ok = 1095-3-7-247 = 838. I'll present these numbers as illustrative with note. Actually let me not present exact suspicious counts as fact; format output block with these numbers and caption "output (illustrative; exact counts depend on the random data)". Good. Also in report, "e.g." dates as computed above: - missing: 2022-02-05, 2022-02-06, 2023-05-01 - negative: 2021-05-01 - above_limit: 2021-08-29 - repeated: 2021-10-28, 2021-10-29, 2021-10-30, 2021-10-31, 2021-11-01 - outlier: e.g. 2021-08-29, ... (order chronological; natural ones unknown) — I'll write "e.g. 2021-08-29, 2022-06-14 ..." hmm inventing; better to write generic "e.g. 2021-08-29, …". Keep vague with "…". - wet_spell: 2022-08-24 ... - dry_spell: 2022-06-20 ... - jump: 2021-08-29, 2021-08-30, 2022-02-10, 2022-02-11 Wait repeated: idx300..304 → DOY 301..305: 2021: DOY301 = Oct 28 (273 end Sep; 274=Oct1 → 301 = Oct 28 ✓); DOY305 = Nov 1 ✓. So dates Oct 28, 29, 30, 31, Nov 1 — 5 days ✓. wet_spell idx600..635 → 2022 DOY 236..271: Aug 24..Sep 28 ✓ (36 days). dry_spell idx900..1094 → 2022 DOY 171..365: Jun 20..Dec 31 ✓ (195 days). Great. Now, also mention in notes: for datasets with trace amounts coded e.g. 0.005 or "-1", convert to 0 first; for accumulation errors (tipping bucket), check dry spells ending with big jump — the jump test helps. Also mention: for near-real-time QC vs historical; iterative approach: after first pass, exclude hard-failed values and recompute the outlier fence (we compute fence from raw data including bad values; quantiles robust so fine, but you can loop). Mention briefly. Also possibly mention `min_count` etc. Fine. Now write the answer: brief intro, bullet list of checks with defaults (maybe table), code block, example usage is embedded in `__main__`, sample output block, then "notes / how to adapt" bullets (thresholds region-specific; enable monthly; spatial checks; flag conventions; save output). Keep the notes short. Check imports: numpy, pandas; plot_qc uses matplotlib inside function (import inside). In __main__ we don't call plot to keep output simple; mention plot function exists. Actually maybe drop plot to reduce length? A small plot function is useful for QC review. I'll include a compact `plot_qc` (~10 lines) and mention it. Or skip plotting to keep answer tighter. Hmm. The user asked for "quality control ... in python" — code focus. I'll include the plotting function briefly since visualization is standard in manual QC verification. Keep it minimal. Actually, to keep the answer focused and not too long, I'll skip matplotlib plotting and instead mention that flagged days can be plotted easily (one-liner hint). Hmm... I think a tiny plot helper adds practical value. Let me include a compact version: ```python def plot_qc(qc, ax=None): import matplotlib.pyplot as plt if ax is None: _, ax = plt.subplots(figsize=(14, 4)) colors = {"ok": "tab:blue", "suspect": "tab:orange", "fail": "tab:red", "missing": "lightgray"} for status, grp in qc.groupby("status"): ax.bar(grp.index, grp["prcp"], width=1.0, color=colors[status], label=status) ax.set_ylabel("rain (mm/day)"); ax.legend() return ax ``` Fine, ~10 lines. Include. One more check on `_flag_monthly_totals`: `v.groupby(v.index.to_period("M")).sum(min_count=1)` — groupby with a PeriodIndex key works. `.sum(min_count=1)` — supported in groupby sum? SeriesGroupBy.sum has min_count param (yes, since 0.21?). ok. `pd.Series(v.index.to_period("M").isin(bad), index=v.index) & v.notna()` — isin on PeriodIndex with `bad` PeriodIndex works. & with bool Series fine. Also in `_flag_outliers`, `v[v > 0]` — NaN > 0 → False, excluded. good. `.dropna()` redundant but harmless (v>0 False for NaN). Actually `v > 0` is False for NaN so v[v>0] excludes NaN. dropna unnecessary; keep or drop. I'll keep without dropna. In `_flag_repeated`, `v.ne(v.shift())` produces bool with True at NaN positions (since NaN != NaN). Note: `v.ne` with NaN: pandas returns True for NaN vs anything (element-wise ne treats NaN as different). Confirm: `pd.Series([1.0, np.nan]).ne(pd.Series([1.0, np.nan]))` → [False, True]. Yes. `run_len = v.groupby(run_id).transform("count")` — count excludes NaN. Within a pure wet run (no NaN), count = length ✓. If a NaN inside (missing day), it increments run_id itself so it's its own group with count 0 — but wet is False for NaN → not flagged ✓. Edge: leading day: v.shift() NaN → ne True → run_id 1. ok. `_spell_lengths` when mask all False: run_id all 0; transform sum = 0; where → 0 ✓. Potential pandas warning: groupby on bool with sum: "sum" of bool → may return bool? SeriesGroupBy.sum on bool returns int64 I believe. Even if bool, comparison with int works (True > 30 False; but run length e.g. 36 can't be bool). Since transform("sum") on bool groups returns int64. Fine. `out[FLAG_COLS] = flags[FLAG_COLS]` — fine. `np.select` returns object array of strings; assign to out["status"]. fine. Now the demo: I need to double check `s.drop(s.index[[1500, 1501, 2345]])` — n=1095 (2021-2023 = 365+365+365 = 1095). Indices 1500, 2345 out of range → error! Fix: choose indices within range: e.g., drop idx 400, 401 (Feb 5-6, 2022) and idx 850 (May 1, 2023). But idx 400/401 are near the jump injection (405/406) — dropping them creates NaN at 400,401; jump check unaffected (403..406 fine). Let me drop [400, 401, 850]. Wait also rain[400], rain[401] were random; dropping them just removes dates. Good. Dates: idx400 → 2022 DOY 36+? 400-365=35 → DOY36 → Feb 5 ✓; idx401 → Feb 6 ✓; idx850 → 850-730=120 → DOY121 → May 1, 2023 ✓. So `s = s.drop(s.index[[400, 401, 850]])`. Hmm wait: 2021-2023: 2021: 365 (idx 0-364), 2022: 365 (idx 365-729), 2023: 365 (idx 730-1094). idx 400 → 2022 day 36 (0-based: 400-365 = 35 → 36th day of year → Feb 5 ✓ since Jan has 31 days: day 32 = Feb 1, day 36 = Feb 5 ✓). idx 850 → 850-730 = 120 → day 121 of 2023 → May 1 ✓ (31+28+31+30=120 → day 121 = May 1 ✓). Also verify repeated dates: idx 300-304 → 2021 day 301-305 → Oct 28 - Nov 1 ✓ (day 274 = Oct 1; day 301 = Oct 28 ✓; day 305 = Nov 1 ✓). wet spell idx 600-635 → 2022: 600-365 = 235 → day 236 → Aug 24 ✓ (day 212 = Jul 31; day 213 = Aug 1; day 236 = Aug 24 ✓); day 271 = Sep 28 ✓ (Aug 31 days: day 243 = Aug 31; day 244 = Sep 1; day 271 = Sep 28 ✓). dry spell idx 900-1094 → 2022: 900-365 = 535 → day 536? wait 900-365 = 535, so day-of-year 536 → but 2022 has 365 days! 900 ≥ 730? No: 2022 spans idx 365..729. idx 900 is in 2023! Let me recompute: 2021: idx 0..364; 2022: idx 365..729; 2023: idx 730..1094. So idx 900 → 2023: 900-730 = 170 → day 171 → Jun 20, 2023 ✓. And idx 1094 → day 365 → Dec 31, 2023 ✓. So dry spell = Jun 20 – Dec 31, 2023 (195 days) ✓. Jump idx 405,406 → 2022: 405-365 = 40 → day 41 → Feb 10 ✓; day 42 → Feb 11 ✓. And idx 240 (5000) → 2021 day 241 → Aug 29 ✓ (day 212 = Jul 31 → day 213 = Aug 1; day 241 = Aug 29 ✓). idx 241 → Aug 30 ✓. idx 120 → day 121 → May 1, 2021 ✓. Wet spell boundary pins: rain[599] = 0 → day 236-1 = day 235 = Aug 23, 2022; rain[636] = 0 → day 272 = Sep 29 ✓. Dry spell pin: rain[899] = 5.0 → day 170 = Jun 19, 2023 ✓. Great, all consistent. So sample output: ``` Daily rainfall QC report ------------------------------------------------------------ records: 1095 ok: 833 suspect: 252 fail: 7 missing: 3 missing 3 e.g. 2022-02-05, 2022-02-06, 2023-05-01 negative 1 e.g. 2021-05-01 above_limit 1 e.g. 2021-08-29 repeated 5 e.g. 2021-10-28, 2021-10-29, 2021-10-30, ... outlier 15 e.g. 2021-08-29, ... wet_spell 36 e.g. 2022-08-24, 2022-08-25, ... dry_spell 195 e.g. 2023-06-20, 2023-06-21, ... jump 4 e.g. 2021-08-29, 2021-08-30, 2022-02-10, ... monthly_total 0 ``` Wait monthly_total: if hits == 0, we skip the line (my loop only prints if len(hits)). So monthly_total line absent. Also outlier count unknown — I'll put a plausible number like 14 and mark the whole output as illustrative. Also "ok" count: 1095 - 3 missing - 7 fail - suspect. suspect = ? Let's compute: days with only-soft flags. wet_spell 36 (all suspect; note any of these also jump? diff between consecutive small gamma values < 150 → no; also outlier? values < fence (fence ~20-25) → no). dry_spell 195. jump 4 (idx 240 hard-fail; idx 241 suspect; idx 405, 406 suspect) → 3 suspect. outlier: idx 240 → fail; others ~13 suspect (assuming 14 total). suspect total = 36 + 195 + 3 + 13 = 247. fail = 7 (negative 1, above_limit 1 (idx240), repeated 5). missing 3. ok = 1095 - 3 - 7 - 247 = 838. I'll present these numbers as illustrative with note. Actually, let me not present exact suspicious counts as fact; format output block with these numbers and caption "output (illustrative; exact counts depend on the random data)". Good enough. Also in report, "e.g." dates as computed above: - missing: 2022-02-05, 2022-02-06, 2023-05-01 - negative: 2021-05-01 - above_limit: 2021-08-29 - repeated: 2021-10-28, 2021-10-29, 2021-10-30, 2021-10-31, 2021-11-01 - outlier: e.g. 2021-08-29, ... (order chronological; natural ones unknown) — I'll write "e.g. 2021-08-29, ...". Keep vague with "…". - wet_spell: 2022-08-24 ... - dry_spell: 2022-06-20 ... - jump: 2021-08-29, 2021-08-30, 2022-02-10, 2022-02-11 Wait repeated: idx300..304 → DOY 301..305: 2021: DOY301 = Oct 28 (273 end Sep; 274=Oct1 → 301 = Oct 28 ✓); DOY305 = Nov 1 ✓. So dates Oct 28, 29, 30, 31, Nov 1 — 5 days ✓. wet_spell idx600..635 → 2022 DOY 236..271: Aug 24..Sep 28 ✓ (36 days). dry_spell idx900..1094 → 2022 DOY 171..365: Jun 20..Dec 31 ✓ (195 days). Great. Now, "Records evaluated: 1095". Counts: negative 1; above_limit 1; repeated 5; outlier maybe ~15; wet_spell 36; dry_spell 195; jump 4; missing 3. status: missing 3; fail = negative(1)+above_limit(1)+repeated(5) = 7 (some may overlap? repeated days distinct from others; outlier idx240 is hard anyway) → fail 7; suspect = outlier(15-1=14? idx240 counted in outlier and fail; suspect counts days with soft-only) + wet_spell 36 + dry_spell 195 + jump(4, minus idx240 which is hard? idx240 flagged jump but status fail; idx241 jump only → suspect; idx405,406 suspect → 3) → suspect ≈ 14 + 36 + 195 + 3 = 248. ok = 1095 - 3 - 7 - 248 = 837. I'll present approximate: ok 837, suspect 248, fail 7, missing 3 — with a caveat "counts vary with the seed". Hmm presenting precise numbers that won't match actual run could be seen as wrong. I'll add a note: "(numbers will vary slightly depending on the random data / seed)". That's honest. Also I realize `rain[900:1100] = 0.0` then `rain[899] = 5.0` — but rain[899] might already be wet random; forcing 5.0 makes it wet ✓. Also check idx 899 not overlapping other injections. fine. Also drop indices [1500,...] — invalid. Change to drop indices within 0..1094: e.g., drop 3 random: rng.choice... but then injected positions might be dropped accidentally (probability tiny but possible: choosing 3 of 1095, injected ~10 positions → ~3% chance). To keep demo deterministic, drop fixed indices [400, 401, 850]? But 400,401 near jump demo (405,406) — dropping them creates missing at Feb 5-6 2022; jump at Feb 10-11 unaffected. Fine. Actually also make sure none of dropped indices coincide with injections: 400,401 no; 850? rain[850] random — fine. But wait: dropping 400,401 → reindex adds NaN → missing flags at Feb 5-6. Also `_spell_lengths` for dry/wet: NaN breaks runs. fine. Let me finalize: `s = s.drop(s.index[[400, 401, 850]])`. Hmm, but idx 400-401 is 2022 Feb 5-6, and jump demo at 405-406 (Feb 10-11) — no conflict. Also make sure the wet spell injection indices 600:636 (2022 Aug 24–Sep 28) — fine. Dry spell 900:1100 → clipped to 1094 → Jun 20–Dec 31 2023, 195 days ✓. Now also double check `rain[300:305] = 2.0` — 5 days (300,301,302,303,304). max_repeated=3 → run_len 5 > 3 → flagged ✓. But also natural random equal-value runs: rounding to 0.1 could create e.g. two consecutive 0.3s; runs of 4+ identical non-zero values by chance: for values concentrated near small (gamma shape 0.6 heavy near 0): probability a wet day equals a specific small value after rounding: e.g., P(0.1) maybe ~10-15%? Two consecutive both 0.1: possible; 4 consecutive 0.1: (0.12)^4 * 0.35^4... extremely unlikely (~2e-6). But runs of exactly 4 identical would flag. Also cross wet/dry: run counts only wet. Also 0.1,0.1,0.1,0.1 needs 4 consecutive wet days each exactly 0.1: p ≈ (0.35*0.12)^3*... ≈ tiny. OK. Also note: rounding rain to 1 decimal — wet values like 0.1 quite common (gamma(0.6,8): P(X<0.15) ≈ γ(0.6, 0.01875)/Γ(0.6): x^a = 0.01875^0.6 = e^{0.6*ln0.01875} = e^{0.6*(-3.976)} = e^{-2.386} = 0.092; γ(a,x) ≈ x^a/a = 0.092/0.6 = 0.153; Γ(0.6) ≈ 1.635 (Γ(0.6) = ?). Hmm Γ(0.6): Γ(1.6)=0.6Γ(0.6) → Γ(1.6)=0.8935 → Γ(0.6)=1.489. So P(X<0.15) ≈ 0.153/1.489 ≈ 0.103. So ~10% of wet days < 0.15 → rounded 0.1. P(wet) 0.35 → P(day=0.1) ≈ 0.035. 4 consecutive: 0.035^4 ≈ 1.5e-6 → over 1095 days ~0.002 expected. negligible. Also check outlier fence won't flag tons: q3+3*IQR of gamma(0.6,8) wet sample (~380 wet days). Estimate quantiles: gamma(0.6,8): compute CDF at x: t=x/8. P(X≤x) = γ(0.6, t)/1.489. - x=1 → t=0.125: t^0.6 = e^{0.6*ln0.125} = e^{0.6*(-2.079)} = e^{-1.2477} = 0.287; γ ≈ 0.287/0.6 * (1 - ...) ≈ 0.478 * (1 - 0.6*... ) lower order corrections; roughly 0.4; P ≈ 0.27. - x=4 → t=0.5: t^0.6 = e^{0.6*(-0.693)} = e^{-0.4159} = 0.6597; γ ≈ (0.6597/0.6)*(1 - 0.6*0.5/1.6 + (0.6*0.5)^2/(2*2.6)*... ) hmm the series: γ(a,x) = x^a Σ_{n=0}∞ [(-x)^n / (n! (a+n))]. terms: n0: 0.6597/0.6=1.0996; n1: -0.5*0.6597/(1.6) = -0.2062; n2: 0.25*0.6597/(2*2.6)=0.0317; n3: -0.125*0.6597/(6*3.6)= -0.0038. sum ≈ 0.9213. γ ≈ 0.9213. P = 0.9213/1.489 ≈ 0.619. So q60 ≈ 4. - x=8 → t=1: γ(0.6,1) = 1^0.6 * Σ: n0: 1/0.6=1.6667; n1: -1/1.6=-0.625; n2: 1/(2*2.6)=0.1923; n3: -1/(6*3.6)=-0.0463; n4: 1/(24*4.6)=0.00905; n5: -1/(120*5.6)=-0.00149. sum ≈ 1.195. P = 1.195/1.489 ≈ 0.803. So q80 ≈ 8. - x=16 → t=2: 2^0.6 = e^{0.4159}=1.516; γ(0.6,2): use γ(a,x)=Γ(a) - Γ(a,x); upper incomplete Γ(0.6,2) = ∫_2^∞ t^{-0.4} e^{-t} dt. Approx via continued fraction... rough: P(X>16) = Γ(0.6,2)/1.489. Γ(0.6,2): compute with series for upper? Use relation Γ(a,x) = x^a E1-ish. Let me numerically estimate: integrand t^{-0.4} e^{-t} from 2 to ∞. At t=2: 2^{-0.4}=e^{-0.277}=0.758; e^{-2}=0.1353 → 0.1026. Decays; approximate integral: substitute u=t: values at 2:0.1026, 3: 3^{-0.4}=e^{-0.439}=0.645; e^{-3}=0.0498 → 0.0321; 4: 4^{-0.4}=e^{-0.5545}=0.5743; e^{-4}=0.0183 → 0.0105; 5: 5^{-0.4}=e^{-0.6438}=0.5253; e^{-5}=0.00674 → 0.00354; 6: 6^{-0.4}=0.488; e^{-6}=0.00248 → 0.00121; 8: 8^{-0.4}=e^{-0.832}=0.435; e^{-8}=0.000335 → 0.000146; 10: 10^{-0.4}=0.398; 4.54e-5 → 1.8e-5. Integral ≈ trapz from 2 to 10 with h=1: (0.1026/2 + 0.0321 + 0.0105 + 0.00354 + 0.00121 + ~0.00047 + 0.000146 + ... + 1.8e-5/2) ≈ 0.0513+0.0321+0.0105+0.00354+0.00121+0.00047+0.000146+0.00003 ≈ 0.0993. tail beyond 10 negligible (~2e-6). So Γ(0.6,2) ≈ 0.099. P(X>16) = 0.099/1.489 ≈ 0.0665. So q93.4 ≈ 16. - x=24 → t=3: Γ(0.6,3): integrand from 3: 3:0.0321; 4:0.0105; 5:0.00354; 6:0.00121; 7: 7^{-0.4}=e^{-0.7779}=0.4595; e^{-7}=0.000912 → 0.000419; 8: 0.000146; 9: 9^{-0.4}=0.4157? e^{-0.879}=0.415; e^{-9}=1.234e-4 → 5.1e-5; 10: 1.8e-5; sum ≈ 0.0321/2? trapz 3..10: 0.01605+0.0105+0.00354+0.00121+0.000419+0.000146+0.000051+0.000018+... ≈ 0.0319. So P(X>24) ≈ 0.0319/1.489 ≈ 0.0214. q97.9 ≈ 24. - x=32 → t=4: integral from 4: 0.0105/2 + 0.00354 + 0.00121 + 0.000419 + 0.000146 + 0.000051 + 0.000018 + ... ≈ 0.00525+0.00354+0.00121+0.000419+0.000146+0.000051+0.000018 ≈ 0.0106. P(X>32) ≈ 0.0106/1.489 ≈ 0.0071. q99.3 ≈ 32. - x=48 → t=6: integral from 6: 0.00121/2 + 0.000419 + 0.000146 + 0.000051 + 0.000018 ≈ 0.000605+0.000419+0.000146+0.000051+0.000018 ≈ 0.00124. P(X>48) ≈ 0.00083. q99.92 ≈ 48. So distribution: q25 ≈ ? x=0.5 → t=0.0625: t^0.6 = e^{0.6*(-2.7726)} = e^{-1.6636} = 0.1894; γ ≈ 0.1894/0.6 = 0.3157*(1 - 0.6*0.0625/1.6 + ...) ≈ 0.3157*(1-0.0234) ≈ 0.308; P ≈ 0.207. x=0.8 → t=0.1: t^0.6 = e^{0.6*(-2.3026)} = e^{-1.3816} = 0.2511; γ ≈ 0.2511/0.6*(1 - 0.6*0.1/1.6 + (0.06)^2/(2*2.6)) = 0.4185*(1-0.0375+0.0007) = 0.4185*0.9632 = 0.4031; P ≈ 0.271. x=1.2 → t=0.15: t^0.6 = e^{0.6*(-1.897)} = e^{-1.138} = 0.3204; γ ≈ 0.3204/0.6 * (1 - 0.6*0.15/1.6 + (0.09)^2/(2*2.6) - ...) = 0.534*(1 - 0.05625 + 0.00156 - 0.00003) = 0.534*0.9453 = 0.5048; P ≈ 0.339. x=1.6 → t=0.2: t^0.6 = e^{0.6*(-1.6094)} = e^{-0.9657} = 0.3807; γ ≈ 0.3807/0.6*(1 - 0.075 + 0.0027 - 0.00007) = 0.6345*0.9276 = 0.5886; P ≈ 0.395. x=2.0 → t=0.25: t^0.6 = e^{0.6*(-1.3863)} = e^{-0.8318} = 0.4353; γ ≈ 0.4353/0.6*(1 - 0.09375 + 0.0042 - 0.00014) = 0.7255*0.9103 = 0.6604; P ≈ 0.4435. x=2.4 → t=0.3: t^0.6 = e^{0.6*(-1.204)} = e^{-0.7224} = 0.4856; γ ≈ 0.4856/0.6*(1 - 0.1125 + 0.00602 - 0.00025) = 0.8093*0.8933 = 0.723; P ≈ 0.4856. So q25 ≈ 0.65, q50 ≈ 2.3ish (P at 2.4 = 0.486; at 2.6: t=0.325, t^0.6 = e^{0.6*ln0.325} = e^{0.6*(-1.124)} = e^{-0.6742} = 0.5095; γ ≈ 0.5095/0.6*(1-0.1219+0.0071-0.0003) = 0.8492*0.8849 = 0.7515; P ≈ 0.5046. so median ≈ 2.55). q75: find x where P=0.75: try x=6: t=0.75: t^0.6 = e^{0.6*(-0.2877)} = e^{-0.1726} = 0.8414; γ(0.6, 0.75): series: n0: 0.8414/0.6 = 1.4023; n1: -0.75*0.8414/1.6 = -0.3944; n2: 0.5625*0.8414/(2*2.6) = 0.0910; n3: -0.4219*0.8414/(6*3.6) = -0.0164; n4: 0.3164*0.8414/(24*4.6) = 0.00242; n5: -0.2373*0.8414/(120*5.6) = -0.0003. sum ≈ 1.4023-0.3944+0.0910-0.0164+0.0024-0.0003 = 1.0846. P = 1.0846/1.489 = 0.728. x=8 → P=0.803 (computed above). Interpolate for 0.75: between 6 (0.728) and 8 (0.803): ≈ 6.6. So q75 ≈ 6.6. q25 ≈ 0.65 (P at 0.65: t=0.08125: t^0.6 = e^{0.6*(-2.5104)} = e^{-1.5062} = 0.2218; γ ≈ 0.2218/0.6*(1 - 0.6*0.08125/1.6 + ...) = 0.3697*(1-0.0305) = 0.3584; P = 0.2407. at 0.5 P=0.207; at 0.8 P=0.271; interpolate 0.25 → x ≈ 0.71). So q25≈0.7, q75≈6.6 → IQR ≈ 5.9 → fence = 6.6 + 3*5.9 ≈ 24.3. P(X > 24.3) ≈ 0.021 → expected outliers among ~383 wet days ≈ 8. Plus injected 5000 (idx240). So outlier count ≈ 9-15 depending on sample. OK so "outlier ≈ 10" plausible. Fine, report number ~10-15 illustrative. Suspect count then ≈ 36+195+3+ (outlier_suspect ≈ 9-14) ≈ 243-248. I'll write approx numbers with caveat. Actually, to avoid looking wrong, in the sample output I can write counts without exact "ok" number... The report prints all. I'll add the note "(illustrative output; counts will vary with the random sample)". Good enough. Alternatively, I could make the demo deterministic-ish by not rounding etc. But natural variability remains. The caveat suffices. Now also handle: `pd.Series` with all-zero or short series → min_wet_days guard; min length guard? Add a check: if len < some minimum (e.g., 30), warn. Eh, keep simple. One more consideration: `flags["missing"]` also True for gap dates (NaT positions) — yes since reindex creates NaN rows. Also `negative` for NaN → False ✓ (NaN < 0 is False). `above_limit` NaN → False ✓. Also should I flag values that are exactly equal to missing codes? covered. Potential issue: if series has timezone or non-datetime index — we convert with pd.to_datetime; if fails, error. Mention requirement: index must be daily dates (or convertible). Now the report function name `qc_report`. Provide `qc.to_csv("qc_flags.csv")` hint in notes. Also mention flag numbers alternative (0/1/9 etc.) if user's database needs numeric codes: one-liner `.astype(int)`. Alright, also the question phrasing "write a quality control of daily rainfall values in python" — they likely want exactly this. Provide well-commented single-file script. Let me also double-check `_spell_lengths` returns int series with where(mask, 0): `mask.groupby(run_id).transform("sum")` returns int64 series indexed like mask; `.where(mask, 0)` — mask boolean series; where keeps values where mask True else 0 ✓. But wait: `transform("sum")` on a boolean Series grouped — in newer pandas, sum of bool group returns int. OK. Also for `check_repeated`: `v.groupby(run_id).transform("size")` — size counts including NaN; our run groups for wet runs contain only equal wet values (dry/missing days each got their own run id because `| ~wet`? wait I wrote `(v.ne(v.shift()) | ~wet).cumsum()`? Let me re-derive: new_run when value changes OR day is not wet. For dry days (wet False): ~wet True → each dry day increments → own run ✓. Missing: ~wet True (since NaN > 0 False) → own run ✓. Wet day with same value as previous wet day: ne False, ~wet False → same run ✓. Wet day with different value: ne True → new run ✓. `transform("size")` counts rows in group = run length ✓ (all members wet, non-NaN). But `transform("size")` — for groupby with transform, "size" is allowed ✓. flag = wet & (run_len > max_repeated) ✓. Edge: consecutive wet days equal value but separated by missing? Missing day breaks run (own run) ✓ good. Now `jump`: also consider jump after missing gap: prev valid day might be long ago; diff across gap still computed on adjacent rows — could flag if gap in between (e.g., 0 then missing then 200). prev valid is the day before missing → diff 200 across 2-day gap → flagged though there's a gap. Acceptable (suspect). Could require both consecutive calendar days: check index adjacency. Refine: compute `gap_ok = (v.index.to_series().diff().dt.days == 1)` and require. Let's add that for correctness: ```python cal_ok = pd.Series(v.index.to_series().diff().dt.days.eq(1), index=v.index) flags["jump"] = ok & ok.shift(fill_value=False) & cal_ok & (v.diff().abs() > thr) ``` Hmm ok.shift(fill_value=False) — "ok" is valid mask; `prev_valid = ok.shift(1, fill_value=False)`. Since reindexed to full daily calendar, diff days always 1 except... reindex guarantees daily frequency, so cal_ok always True. So unnecessary after reindex. Skip it. Good — reindex makes calendar contiguous, so simple diff fine. But mention that series is reindexed first. Also note: after reindex, gaps became NaN rows; jump across gap: prev row is NaN → ok False → not flagged ✓. Nice. Also `wet_spell`: missing days break runs (mask False) ✓. Also `dry_spell` counts zeros only; missing breaks. A very long dry spell including missing days would be split — mention in comment. Now let me think about whether `flags["missing"]` should also incorporate `isna`: yes. Also handle `prcp` given as DataFrame column name? Accept Series. Provide usage: `qc = run_qc(df["rain_mm"], cfg)`. Also allow `freq="D"`; if the series is at other frequency, reindex daily anyway. Fine. Write final code now, ensuring syntax correctness. ```python """ QC of daily rainfall (precipitation) — pandas implementation. Runs a battery of independent checks and returns a flag table: check what it detects ---------------- ---------------------------------------------------- missing NaN / missing codes (-999 ...) / dates absent from index negative negative values (sensor/processing error) above_limit values above a physically plausible daily maximum repeated >= N consecutive days with the same non-zero value (stuck tipping bucket / datalogger) outlier statistical outlier among wet days (Tukey fence Q3 + k*IQR) wet_spell uninterrupted run of wet days longer than N dry_spell uninterrupted run of dry days longer than N jump |P(t) - P(t-1)| above threshold (both days valid) Each check is an independent boolean column: True = failed the check. 'status' summarises: ok / suspect (soft flags) / fail (hard errors) / missing. """ from __future__ import annotations import numpy as np import pandas as pd # -------------------------------------------------------------------------- # configuration (tune to your climate / instrument) # -------------------------------------------------------------------------- DEFAULTS = dict( missing_codes=(-999.0, -99.0, -9999.0), upper_limit=1000.0, # mm/day (world record ~1825 mm, Foc-Foc, Réunion 1966) max_repeated=3, # flag if > 3 consecutive identical non-zero values iqr_k=3.0, # Tukey "extreme" fence multiplier (3*IQR) min_wet_days=50, # outlier test needs at least this many wet days max_wet_spell=30, # days; > 30 consecutive wet days are flagged max_dry_spell=180, # days; > 180 consecutive dry days are flagged jump_threshold=150.0, # mm; |day-to-day change| larger than this is suspect ) ``` Wait, I earlier used max_dry_spell 180 and demo 195-day dry spell → flagged ✓. Good, consistent (default 180). And max_wet_spell 30, demo 36 → flagged ✓. min_wet_days 50 — sample has ~380 wet days ✓. Hmm, earlier cfg had max_dry_spell 365; now 180. Wet spell flag "run > 30" means 31+ days flagged. 36-day demo run flagged ✓. `FLAG_COLS = ["missing", "negative", "above_limit", "repeated", "outlier", "wet_spell", "dry_spell", "jump", "monthly_total"]` — but monthly_total only if enabled; keep column always (False when disabled). Add cfg key `max_monthly_total=None` (disabled). In demo, leave disabled (monthly totals ~ sum of wet days ~ 0.35*30*4.8 ≈ 50 mm typical; max month maybe 300-400 — under 5000). Fine. Also add `hard_flags = ("negative", "above_limit", "repeated")`, `soft_flags = ("outlier", "wet_spell", "dry_spell", "jump", "monthly_total")`. Main function: ```python def run_qc(prcp: pd.Series, cfg: dict | None = None) -> pd.DataFrame: """Run all checks; return a DataFrame with the values, per-check flags and status.""" cfg = {**DEFAULTS, **(cfg or {})} s = prcp.astype(float).copy() s.index = pd.to_datetime(s.index) s = s.sort_index() if s.index.has_duplicates: s = s[~s.index.duplicated(keep="first")] # full daily calendar: gaps in the record become NaN rows -> flagged 'missing' s = s.reindex(pd.date_range(s.index.min(), s.index.max(), freq="D")) flags = pd.DataFrame(index=s.index) flags["missing"] = s.isna() | s.isin(list(cfg["missing_codes"])) flags["negative"] = s < 0 flags["above_limit"] = s > cfg["upper_limit"] flags["repeated"] = _check_repeated(s, cfg["max_repeated"]) flags["outlier"] = _check_outlier(s, cfg["iqr_k"], cfg["min_wet_days"]) flags["wet_spell"] = _check_spell(s > 0, cfg["max_wet_spell"]) flags["dry_spell"] = _check_spell(s == 0, cfg["max_dry_spell"]) prev = s.shift(1) flags["jump"] = s.notna() & prev.notna() & (s - prev).abs() > cfg["jump_threshold"] ``` Careful with operator precedence: `s.notna() & prev.notna() & (s - prev).abs() > cfg["jump_threshold"]` — comparison binds tighter than &? In Python, comparison has lower precedence than bitwise &? Actually precedence: bitwise & is lower than comparisons? Order (high→low): **, unary, * /, + -, <<, &, ^, |, comparisons, not, and, or. Wait: comparisons are LOWER precedence than bitwise &! So `a & b > c` parses as `a & (b > c)`. Yes — in Python, `&` binds tighter than comparisons. So `s.notna() & prev.notna() & (s - prev).abs() > thr` → `((s.notna() & prev.notna()) & ((s-prev).abs())) > thr`?? Hmm: `&` higher precedence than `>`: expression `A & B > C` = `(A & B) > C`. So the whole chain: `s.notna() & prev.notna() & (s-prev).abs() > thr` parses as `(((s.notna() & prev.notna()) & (s-prev).abs()) > thr)` — bitwise and of booleans with float → error or weird. Must parenthesize: `flags["jump"] = s.notna() & prev.notna() & ((s - prev).abs() > cfg["jump_threshold"])`. Yes parenthesize the comparison. ✓. `_check_spell(mask, limit)`: ```python def _spell_lengths(mask: pd.Series) -> pd.Series: """Length of the run of consecutive True containing each True element.""" starts = mask & ~mask.shift(fill_value=False) run_id = starts.cumsum() return mask.groupby(run_id).transform("sum").where(mask, 0) ``` Hmm wait: for False positions, run_id equals current run number (of the preceding True run) or 0 at start; group includes those False rows; transform sum counts True in group = run length; .where(mask, 0) zeroes them ✓. But careful: if the series STARTS with False: run_id = 0 for those; group 0 = all leading Falses (sum 0) ✓. But if mask all False: starts empty, run_id all 0, transform sum = 0 ✓. Return type int. Then `_check_spell(mask, limit) = mask & (_spell_lengths(mask) > limit)`. But `_spell_lengths` needs mask as pd.Series with bool dtype: `s > 0` gives bool with NaN→False ✓. `s == 0` NaN→False ✓. _check_repeated: ```python def _check_repeated(s: pd.Series, max_repeated: int) -> pd.Series: wet = s > 0 if not wet.any(): return pd.Series(False, index=s.index) new_run = wet & ((s != s.shift()) | s.isna().shift(fill_value=False)) ``` Hmm, simpler as designed: `(s.ne(s.shift()) | ~wet)`. Wait need care: s.ne(s.shift()) with NaN: pandas returns True for NaN comparisons ✓ (NaN != NaN). So: ```python new_run = wet | s.ne(s.shift()) | s.isna() ``` Hmm let me think: run_id should increment when: day not wet (dry or missing), or value changed from previous day, or previous day missing. Since missing days are "not wet", `~wet` covers it. So: ```python new_run = ~wet | s.ne(s.shift()) run_id = new_run.cumsum() run_len = s.groupby(run_id).transform("size") return wet & (run_len > max_repeated) ``` Wait — but for a wet run, group = consecutive wet days with identical values. Does `s.ne(s.shift())` correctly stay False within? If value constant and all wet: ne False each interior day; first day of series: shift NaN, ne True → new run ✓. Two wet days with values 3.0, 3.0: ne False → same run ✓. Then 0.0 next: ~wet True → new run ✓. But hold on: wet days where value changes (e.g., 2.0 → 5.0): ne True → new run ✓. run_len via transform("size") — group includes only that run's rows? A run id g spans from its first row until next increment. Since every non-wet day increments, groups alternate: [run of wet equal values], [single dry/missing day], [next run], ... So group size = run length ✓. But wait: what about wet days with equal value followed by wet days with different value — each distinct-value stretch is its own run ✓. Edge: first row if wet: shift NaN → ne True → new run ✓. If first row dry: ~wet True → new run (id 1), run contains just that dry day? Next row wet with ne (0 != value) True → new run ✓. But subtle: a wet run where a NaN sneaks in? NaN is not wet → increments ✓. Also `transform("size")` on groupby of s by run_id: s has NaN values in dry/missing groups; size counts rows regardless ✓. flag = wet & run_len > max_repeated. For the demo: run of 5 days at 2.0 → size 5 > 3 → all 5 flagged ✓. Hmm, one more: should flag entire run or only days after the k-th? RClimDex flags whole sequence. Keep whole run; note in docstring. _check_outlier: ```python def _check_outlier(s: pd.Series, k: float, min_wet_days: int) -> pd.Series: flag = pd.Series(False, index=s.index) wet = s[s > 0].dropna() if wet.size < min_wet_days: return flag q1, q3 = wet.quantile([0.25, 0.75]) iqr = q3 - q1 if iqr <= 0: return flag # not enough spread for the test (e.g. all equal values) upper = q3 + k * iqr return s > upper ``` `s > upper` NaN→False ✓. Monthly check: ```python def _check_monthly(s: pd.Series, max_monthly) -> pd.Series: flag = pd.Series(False, index=s.index) if max_monthly is None: return flag valid = s.where(s >= 0) # ignore flagged negatives/missing? hmm totals = valid.groupby(s.index.to_period("M")).transform("sum") return valid.notna() & (totals > max_monthly) ``` Using transform sum per month mapped to days. Negative values would reduce totals; clamp negatives to 0: `valid = s.clip(lower=0)` but NaN stays NaN? clip keeps NaN ✓. So: ```python totals = s.clip(lower=0).groupby(s.index.to_period("M")).transform("sum") return s.notna() & (totals > max_monthly) ``` Flag all non-missing days in offending months. Good. Status: ```python HARD = ("negative", "above_limit", "repeated") SOFT = ("outlier", "wet_spell", "dry_spell", "jump", "monthly_total") ... flags["n_flags"] = flags[FLAG_COLS].sum(axis=1) flags["status"] = np.select( [flags["missing"], flags[list(HARD)].any(axis=1), flags[list(SOFT)].any(axis=1)], ["missing", "fail", "suspect"], default="ok", ) out = pd.DataFrame({"prcp": s}).join(flags) ``` Wait flags already has columns; build out = pd.concat([s.rename("prcp"), flags], axis=1). Since flags index = s.index ✓. Also note: np.select condlist uses flags[...] pandas bool Series — np.select works with arrays; convert `.to_numpy()` maybe. np.select accepts arrays; pandas Series works (converted). Fine. Then `qc_report(qc)`: ```python def qc_report(qc: pd.DataFrame, n_examples: int = 6) -> str: lines = [f"Daily rainfall QC — {len(qc)} records " f"({qc.index.min():%Y-%m-%d} → {qc.index.max():%Y-%m-%d})", "-" * 72] for col in FLAG_COLS: hits = qc.index[qc[col]] if len(hits): ex = ", ".join(hits[:n_examples].strftime("%Y-%m-%d")) more = " …" if len(hits) > n_examples else "" lines.append(f"{col:<12} {len(hits):5d} e.g. {ex}{more}") counts = qc["status"].value_counts().reindex(["ok", "suspect", "fail", "missing"], fill_value=0) lines.append("-" * 72) lines.append("status: " + " ".join(f"{k}={v}" for k, v in counts.items())) return "\n".join(lines) ``` Wait — if a check has zero hits, skip (cleaner). fine. `hits[:n_examples].strftime` — DatetimeIndex.slice then .strftime returns Index of str ✓. Demo section: ```python if __name__ == "__main__": rng = np.random.default_rng(42) idx = pd.date_range("2020-01-01", "2022-12-31", freq="D") n = idx.size rain = np.zeros(n) wet = rng.random(n) < 0.35 # ~35% wet days rain[wet] = rng.gamma(0.6, 8.0, wet.sum()) # skewed wet-day amounts rain = np.round(rain, 1) # --- inject typical errors ------------------------------------------- rain[100] = -2.0 # negative value rain[200] = 5000.0 # impossible magnitude rain[300:305] = 2.0 # 5 identical days -> stuck gauge rain[400:402] = [0.0, 200.0] # day-to-day jump rain[600:637] = rng.gamma(0.7, 3.0, 37) # 37-day wet spell rain[900:1100] = 0.0 # 200-day dry spell s = pd.Series(rain, index=idx, name="rain_mm") s = s.drop(s.index[[1500, 1501, 2345]]) # gaps in the record qc = run_qc(s) print(qc_report(qc)) # inspect / export bad = qc[qc["status"].isin(["fail", "suspect"])] # bad.to_csv("rainfall_qc_flags.csv") ``` Careful: `rain[400:402] = [0.0, 200.0]` — rain[400] might naturally be wet; overriding to 0; rain[402] natural. jump between 401(200) and 402(natural, maybe 0 or small) → |200 - x| > 150 if x < 50 → flagged at 402 as well. Also jump 400→401: |0-200|=200 flagged at 401. And 399→400: |x - 0| maybe <150. So jump count 2-ish ✓. wet spell 600:637 = 37 days ✓ > 30 flagged. Values gamma(0.7, 3) mean 2.1; all > 0 ✓ (gamma positive). Rounded? not rounded → floats fine. Distinct values mostly; if two adjacent equal, repeated runs small (need >3 identical) — unlikely with continuous values. But rounding? I won't round the spell values. But whole `rain` rounded earlier then overwritten — fine. Also idx 2345 within range (n=1096? date_range 2020-01-01..2022-12-31 = 365+365+366 = 1096 days (2020 leap)). 2345 > 1095! Out of range → drop index error. Fix: use drop indices within range: e.g., 500, 700, 900? But 900-1100 is dry spell injection — dropping inside it still fine (dry spell 200 - 1 gap... drop at 900 and 901 would split dry run into two: before drop run from ...? Actually dry spell 900:1100 zeros; random data before 900 may also be dry randomly. Dropping 500/700 arbitrary random days — creates gaps elsewhere ✓. Let me drop indices [500, 501, 800] — but 800 might be wet random; dropping fine. Also ensure drops not in injected error zones (100,200,300-304,400-401,600-636,900-1099): 500 ok, 501 ok, 800 ok ✓. So `s.drop(s.index[[500, 501, 800]])` → 3 missing days ✓. Hmm wait, dry spell 900:1100: n=1096 → rain[900:1100] = 900..1095 → 196 days zeros (to end). Dry run length 196 > 180 ✓ flagged. But also is 1095 (last day) end — run length 196 ✓. But also: might random data before 900 also be dry, extending run beyond 196 — fine, flagged anyway. Wet spell at 600:637; also 400-401 area [0,200]; 402 natural. fine. repeated: 300:304 → 5 days of 2.0 ✓ flagged; also check run 299 (natural) if wet and equal 2.0? natural value rounded could be 2.0 by chance (P ~ small). fine. negative idx 100: if 100 falls in a wet random day, overridden to -2.0 → negative flag ✓. Also `wet_spell`/`dry_spell` unaffected materially. above_limit idx 200 → 5000 flagged above_limit ✓; also jump neighbors: diff 5000-prev > 150 → jump flag at 200 too, and at 201 (5000 - next) → jump also. fine (expected: big spike triggers jump too). Now sample output — I'll present approximate numbers with a caveat note. Let me estimate counts: - missing: 3 - negative: 1 - above_limit: 1 - repeated: 5 - outlier: ~8-12 (say 9). Hmm can't know. I'll write "e.g." numbers and add note. Alternatively, in the sample output block, I can just show the format with plausible numbers and explicitly say "output will look like (numbers vary with the random seed)". That's honest. Let me pick plausible numbers: missing 3; negative 1; above_limit 1; repeated 5; outlier 9; wet_spell 37; dry_spell 196; jump 4 (idx 200 spike & 201, 401, 402). Actually jump at 200: |5000 - rain[199]| > 150 unless rain[199] > 4850 (no) → flagged. At 201: |rain[201] - 5000| > 150 unless rain[201] > 4850 → flagged. At 401 flagged, 402 flagged if rain[402] < 50 (P ~ 0.9 given wet 35% & gamma: if dry 0 → flagged; if wet small → flagged; if wet > 50 → not; P(wet>50) tiny) → likely flagged. So jump ≈ 4. status counts: fail = negative+above_limit+repeated = 7 (distinct days, no overlap) → 7. missing 3. suspect = outlier(9)+wet_spell(37)+dry_spell(196)+jump(4) minus overlaps: wet_spell days within 600-636; jump days at 200,201,401,402 not in spell; outlier days: random wet days incl. maybe within dry spell? no (dry spell zeros, wet>0 false). outlier among wet days could coincide with wet_spell days (600-636 gamma values max ~ maybe 15 < fence 24) unlikely. So suspect ≈ 9+37+196+4 = 246 (maybe minus 1-2 overlaps). ok = 1096 - 3 - 7 - ~246 ≈ 840. Rather than precise, I'll show the report with these approximations and add a line "(your numbers will differ slightly — random data)". Or maybe better to not fabricate a full output; instead show truncated example. I'll include a short illustrative output block and note variability. Also mention how to handle results: keep flagged values but use flag column downstream, or set flagged to NaN, etc. Also mention extension ideas: spatial checks with neighbors, double-mass for drift, homogeneity tests, using `xclim`/`climpact` libraries for indices, and that thresholds should be region-specific. Also mention trace handling: e.g., if your data codes traces as -1 or 0.1, convert: `s = s.replace(trace_code, 0.0)`. I think the code is solid. Let me also double check `np.select` usage: condlist of pandas Series — np.select requires arrays; it will convert via np.asarray — bool Series → bool array ✓. choicelist strings, default "ok" ✓. Result assigned to flags["status"] with index alignment? np.select returns ndarray; assigning to DataFrame column uses index ✓. flags["n_flags"] = flags[FLAG_COLS].sum(axis=1) — bool sum → int ✓. Also `run_qc` should ensure FLAG_COLS all present even if monthly disabled: yes we create the column. One more: `_check_spell(s > 0, cfg["max_wet_spell"])` — mask type: `s > 0` is Series bool ✓; `mask.shift(fill_value=False)` ✓; groupby on run_id (int series) with transform("sum") on bool — returns int64 ✓; `.where(mask, 0)` ✓. Edge: empty series → skip? If s empty, min/max fail. Add guard: if s.empty raise ValueError. Minor; add. Also `s.index.has_duplicates` after sort; handle. Also convert `-999` codes: `s.isin(list(...))` — but if values stored as int -999 ✓ (isin with floats matches ints via equality? `pd.Series([-999]).isin([-999.0])` → True ✓). Note: isin on NaN → False; NaN caught by isna ✓. Alright. Also potential pitfall: `s.index.min()` on empty → guard. Let me also think about `min_wet_days` default: 50. Demo ~380 wet days ✓. Also in `_check_outlier`, use quantile interpolation default linear ✓. Now think about whether to use `freq="D"` vs "MS" for monthly groupby period "M" ✓. Also, in `qc_report`, using f-string with `%Y-%m-%d` format on DatetimeIndex min: `f"{qc.index.min():%Y-%m-%d}"` ✓. Docstrings: concise. Also the user may have a CSV: show loading snippet: ```python df = pd.read_csv("rain.csv", parse_dates=["date"], index_col="date") qc = run_qc(df["rain_mm"]) print(qc_report(qc)) ``` Include in demo/comments. Now write the final answer: brief intro listing checks, the code block, example usage & sample output, and short notes on tuning/extending. Given code-heavy request, markdown with headers and code blocks is appropriate. Keep prose tight. Length: The code is ~150 lines; plus explanation. That's appropriate for "write a quality control" request — they asked for code, so deliver complete code. Avoid over-explaining; bullet the checks. Let me finalize code text carefully, checking syntax line by line. ```python """ Quality control (QC) for daily rainfall series. Checks (each returns a boolean flag; True = failed): missing NaN, missing codes (-999, ...) or dates missing from the index negative negative values (hard error) above_limit value above a physically possible daily maximum (hard error) repeated more than N consecutive days with an identical non-zero value (stuck tipping bucket / logger) (hard error) outlier wet-day outlier, Tukey fence Q3 + k*IQR (suspect) wet_spell uninterrupted wet run longer than N days (suspect) dry_spell uninterrupted dry run longer than N days (suspect) jump |P(t)-P(t-1)| above a threshold, both days valid (suspect) monthly_total (optional) calendar month total above a maximum (suspect) 'status' column summarises each day: ok / suspect / fail / missing. """ from __future__ import annotations import numpy as np import pandas as pd # --------------------------------------------------------------------------- # # configuration — tune the thresholds to your climate / instrument # # --------------------------------------------------------------------------- # DEFAULTS = dict( missing_codes=(-999.0, -99.0, -9999.0), upper_limit=1000.0, # mm/day; world record ≈ 1825 mm (Foc-Foc, Réunion 1966) max_repeated=3, # flag runs of >3 identical non-zero values iqr_k=3.0, # Tukey "far out" fence: Q3 + k*IQR on wet-day amounts min_wet_days=50, # outlier test needs at least this many wet days max_wet_spell=30, # >30 consecutive wet days -> flag max_dry_spell=180, # >180 consecutive dry days -> flag jump_threshold=150.0, # mm; |change between consecutive days| max_monthly_total=None, # e.g. 5000 mm/month for extreme monsoon sites; None = off ) HARD = ("negative", "above_limit", "repeated") SOFT = ("outlier", "wet_spell", "dry_spell", "jump", "monthly_total") FLAGS = ("missing",) + HARD + SOFT ``` Then helpers: ```python def _spell_lengths(mask: pd.Series) -> pd.Series: """For each day, the length of the run of consecutive True it belongs to.""" run_id = (mask & ~mask.shift(fill_value=False)).cumsum() return mask.groupby(run_id).transform("sum").where(mask, 0) def _check_spell(mask: pd.Series, max_len: int) -> pd.Series: return mask & (_spell_lengths(mask) > max_len) def _check_repeated(s: pd.Series, max_repeated: int) -> pd.Series: """Flag runs of identical non-zero values (zeros are never flagged).""" wet = s > 0 run_id = (~wet | s.ne(s.shift())).cumsum() # new run when value changes or day not wet run_len = s.groupby(run_id).transform("size") return wet & (run_len > max_repeated) ``` Wait: `s.ne(s.shift())` — for first element shift is NaN; ne → True. But also careful: `s.ne(s.shift())` where s has NaN: current NaN vs prev value → ne True (NaN != x) ✓; current value vs prev NaN → True ✓; NaN vs NaN → True ✓ (pandas ne treats NaN as unequal). So missing days always start new runs; also `~wet` True for missing anyway. Good. `transform("size")` — group sizes ✓. ```python def _check_outlier(s: pd.Series, k: float, min_wet_days: int) -> pd.Series: """Tukey upper fence applied to wet-day amounts only.""" wet = s.dropna() wet = wet[wet > 0] if wet.size < min_wet_days: return pd.Series(False, index=s.index) q1, q3 = wet.quantile([0.25, 0.75]) iqr = q3 - q1 if iqr <= 0: # e.g. all wet days identical: test not applicable return pd.Series(False, index=s.index) return s > (q3 + k * iqr) def _check_monthly(s: pd.Series, max_monthly: float | None) -> pd.Series: flag = pd.Series(False, index=s.index) if max_monthly is None: return flag totals = s.clip(lower=0).groupby(s.index.to_period("M")).transform("sum") return s.notna() & (totals > max_monthly) ``` Main: ```python def run_qc(prcp: pd.Series, cfg: dict | None = None) -> pd.DataFrame: """Run all checks on a daily series; returns values + flag columns + status.""" cfg = {**DEFAULTS, **(cfg or {})} s = prcp.astype(float).copy() s.index = pd.to_datetime(s.index) s = s.sort_index() if s.index.has_duplicates: # keep first occurrence, note in report s = s[~s.index.duplicated(keep="first")] if s.empty: raise ValueError("empty series") # reindex to a complete daily calendar: record gaps become NaN -> 'missing' s = s.reindex(pd.date_range(s.index.min(), s.index.max(), freq="D")) flags = pd.DataFrame(index=s.index) flags["missing"] = s.isna() | s.isin(list(cfg["missing_codes"])) flags["negative"] = s < 0 flags["above_limit"] = s > cfg["upper_limit"] flags["repeated"] = _check_repeated(s, cfg["max_repeated"]) flags["outlier"] = _check_outlier(s, cfg["iqr_k"], cfg["min_wet_days"]) flags["wet_spell"] = _check_spell(s > 0, cfg["max_wet_spell"]) flags["dry_spell"] = _check_spell(s == 0, cfg["max_dry_spell"]) prev = s.shift(1) flags["jump"] = s.notna() & prev.notna() & ((s - prev).abs() > cfg["jump_threshold"]) flags["monthly_total"] = _check_monthly(s, cfg["max_monthly_total"]) flags["n_flags"] = flags[list(FLAGS)].sum(axis=1) flags["status"] = np.select( [flags["missing"], flags[list(HARD)].any(axis=1), flags[list(SOFT)].any(axis=1)], ["missing", "fail", "suspect"], default="ok", ) return pd.concat({"prcp": s, **{}}, axis=1)... ``` Simpler: `out = flags.copy(); out.insert(0, "prcp", s)` — but flags index equals s.index ✓. Use `pd.concat([s.rename("prcp"), flags], axis=1)`. Report: ```python def qc_report(qc: pd.DataFrame, n_examples: int = 6) -> str: lines = [ f"Daily rainfall QC | {len(qc)} days " f"({qc.index.min():%Y-%m-%d} to {qc.index.max():%Y-%m-%d})", "-" * 74, ] for col in FLAGS: hits = qc.index[qc[col]] if hits.size: examples = ", ".join(hits[:n_examples].strftime("%Y-%m-%d")) dots = " ..." if hits.size > n_examples else "" lines.append(f"{col:<13}{len(hits):>6} e.g. {examples}{dots}") counts = qc["status"].value_counts().reindex( ["ok", "suspect", "fail", "missing"], fill_value=0) lines.append("-" * 74) lines.append("status: " + " ".join(f"{k}={v}" for k, v in counts.items())) return "\n".join(lines) ``` Demo: ```python if __name__ == "__main__": rng = np.random.default_rng(42) idx = pd.date_range("2020-01-01", "2022-12-31", freq="D") n = idx.size rain = np.zeros(n) wet = rng.random(n) < 0.35 # ~35% wet days rain[wet] = rng.gamma(0.6, 8.0, int(wet.sum())) # skewed wet-day amounts rain = np.round(rain, 1) # --- inject typical errors ------------------------------------------- rain[100] = -2.0 # negative value rain[200] = 5000.0 # impossible magnitude rain[300:305] = 2.0 # 5 identical days -> stuck gauge rain[400:402] = [0.0, 200.0] # day-to-day jump rain[600:637] = rng.gamma(0.7, 3.0, 37) # 37-day wet spell rain[900:1100] = 0.0 # 200-day dry spell s = pd.Series(rain, index=idx, name="rain_mm") s = s.drop(s.index[[500, 501, 800]]) # gaps in the record qc = run_qc(s) print(qc_report(qc)) # inspect / export # flagged = qc[qc["status"] != "ok"] # flagged.to_csv("rainfall_qc_flags.csv") # Real data: # df = pd.read_csv("station.csv", parse_dates=["date"], index_col="date") # qc = run_qc(df["rain_mm"]) ``` Check: `rain[400:402] = [0.0, 200.0]` — list assign to numpy slice ✓. `rain[wet] = rng.gamma(...)` ✓. Now the sample output: I'll write an illustrative block. Let me compute dates for injections (2020-01-01 = idx 0, leap year): - idx 100 → 2020: DOY 101 → Jan 31, Feb 29 (2020 leap → DOY 60 = Feb 29; DOY 91 = Mar 31; DOY 101 = Apr 10). So 2020-04-10. - idx 200 → DOY 201 → Mar31(91)+Apr30(121)+May31(152)+Jun30(182)+Jul31(213) → DOY 201 = Jul 19. 2020-07-19. - idx 300:304 → DOY 301-305 → after Jun30(182), Jul31(213): DOY 301 = Oct 27 (213+88 = 301; Aug 31→244, Sep 30→274, Oct 27→301). So 2020-10-27 to 2020-10-31. - idx 400,401 → DOY 401-402: 2020 has 366 days → DOY 401 = 2021 DOY 35 → Feb 4 (Jan31 → DOY 35 = Feb 4). So 2021-02-04 and 2021-02-05. - idx 600:636 → DOY 601-637: 2021 DOY 601-366=235 → 2021: 235th day = Aug 23 (Jan31,Feb59,Mar90,Apr120,May151,Jun181,Jul212,Aug235). So 2021-08-23 to 2021-09-28 (37 days: Aug23 + 36 = Sep 28 ✓ since Aug has 31: Aug23..Aug31 = 9 days, Sep 1..28 = 28 days → 37 ✓). - idx 900:1095 → 2021 DOY 535.. : 900-366=534 → 2021 DOY 535 = Jun 19? Let me compute: cumulative 2021: Jan31, Feb59, Mar90, Apr120, May151, Jun181, Jul212, Aug243, Sep273, Oct304, Nov334, Dec365. DOY 535 > 365 → 2022. 535-365=170 → 2022 DOY 170 = Jun 19 (Jan31,Feb59,Mar90,Apr120,May151,Jun170 → Jun 19 ✓). So dry spell 2022-06-19 to end 2022-12-31 (196 days: idx 900..1095). - drops: idx 500 → 2021 DOY 135 = May 15 (Apr120 + 15). idx 501 → May 16. idx 800 → 800-366=434 → 2022 DOY 435? wait 800+1=801st day overall; overall DOY 801: 2020 has 366 → 801-366=435 → 2022 DOY 435 → 435-365=70? No wait: 2022 DOY = 801-366-365 = 70 → Mar 11 (Jan31,Feb59,Mar70 → Mar 11 ✓ 2022 non-leap). So drops: 2021-05-15, 2021-05-16, 2022-03-11. So report examples: missing e.g. 2021-05-15, 2021-05-16, 2022-03-11; negative 2020-04-10; above_limit 2020-07-19; repeated 2020-10-27...; wet_spell 2021-08-23...; dry_spell 2022-06-19...; jump: 2020-07-19 (jump due to 5000 spike: |5000 - rain[199]| >150 ✓), 2020-07-20 (|rain[201] - 5000| ✓), 2021-02-05, 2021-02-06 (if rain[402] < 50 — random; seed 42 unknown). Hmm 402 = DOY 403 → 2021 Feb 6. If rain[402] natural wet > 50 (unlikely) then no. I'll present "e.g." with a note. Also outlier count unknown; show e.g. "outlier 9" with generic examples like 2020-05-12... — I'd better not fabricate specific dates for random ones. In the sample output, for random-driven counts I can put a placeholder or say numbers vary. To keep honest, I'll caption: "illustrative output (random-data counts will vary between runs)". That's fine. Total days: date_range 2020-01-01 to 2022-12-31 = 366+365+365 = 1096. Report shows 1096 days ✓ (after reindex, includes 3 gap days → 1096 total, 1093 observed). status counts: fail=7 (1 neg + 1 limit + 5 repeated), missing=3, suspect = outlier + wet_spell(37) + dry_spell(196) + jump(~4) ≈ 246ish, ok = rest ≈ 840. I'll write numbers as e.g. ok=840, suspect=246, fail=7, missing=3 with the caveat. Hmm, one more consideration: the wet_spell at 600:636 — idx 599 natural may be wet (35%) → run extends to 38 days. Also idx 637 natural wet → extends. So wet_spell count maybe 37-40. Whatever — illustrative. Also note dry_spell: 196 zeros injected; but is rain[899] zero too? random. Fine. Also natural random wet days: with p=0.35 over 1096 days, some natural runs of identical rounded values could trigger repeated? Rounded to 1 decimal: values like 0.6 repeated 4+ consecutive days: probability low but possible over ~380 wet days... Each wet value from gamma rounded: P(two specific adjacent wet days equal) ≈ maybe 2-4% (rounding buckets near low values where density high: many values 0.1-1.0 → bucket width 0.1 → decent collision probability for small values). Run of 4+ identical requires same wet-mask and same value: p ≈ (0.35 * p_val)^3 * 0.35 ... let me not worry; with seed fixed it's deterministic; might add a couple extra repeated hits. It's illustrative anyway. Actually hmm — this could be a real phenomenon: rounded small values (0.1, 0.2) collide often. With gamma(0.6,8): P(X < 0.5) = ? Γ(0.6, 0.0625)/Γ(0.6) ≈ (0.0625^0.6)/0.6 / 1.489 ≈ (e^{0.6 ln 0.0625}) = 0.0625^0.6 = e^{0.6*(-2.7726)} = e^{-1.6636} = 0.1894; /0.6 = 0.3157; /1.489 = 0.212. So ~21% of wet days are < 0.5mm → buckets 0.1-0.4 each ~5%. Two adjacent wet days both in bucket 0.1: 0.0025 (0.25%)... chain of 4: 0.35^4 * 0.05^4? ≈ 0.015 * 6.25e-6 ≈ 1e-7. Negligible. OK. But wait: rounding to 1 decimal creates many exact ties among small values — but for repeated check we need CONSECUTIVE equal values AND consecutive wet days. Fine. Also, `wet = rng.random(n) < 0.35` — with seed 42, deterministic; whatever. One more check on `_spell_lengths` correctness with a concrete tiny example: mask = [T, T, F, T]. shift(fill_value=False): [F, T, T, F]. starts = mask & ~shift = [T, F, F, T]. cumsum = [1, 1, 1, 2]. groupby(run_id).transform("sum") on mask: group 1 = rows 0,1,2 (mask T,T,F) → sum = 2; group 2 = row 3 → 1. → [2, 2, 2, 1]. .where(mask, 0) → [2, 2, 0, 1] ✓. _check_repeated example: s = [0, 2, 2, 2, 2, 0, 5]. wet = [F,T,T,T,T,F,T]. s.ne(s.shift()): shift = [NaN,0,2,2,2,2,0] → ne = [T, T, F, F, F, T, T]. ~wet = [T,F,F,F,F,T,F]. new_run = ~wet | ne = [T, T, F, F, F, T, T] → run_id = [1,2,2,2,2,6,7]. transform size: group2 size 4 → [.,4,4,4,4,.,.] ; group6 size 1; group7 size 1. flag = wet & (size>3) → days 1-4 (indices 1..4) ✓ 4 days of value 2 flagged with max_repeated=3 ✓. Check outlier: s has NaN? fine. Check jump: prev = s.shift(1); s.notna() & prev.notna() & |s - prev| > thr ✓. At row after a gap (missing), prev NaN → skip ✓. Now monthly: `s.clip(lower=0)` — NaN preserved ✓; groupby period transform sum → NaN months? If all NaN in month, sum → 0? transform("sum") of all-NaN group: pandas sum of all-NaN = 0.0 (sum skipna default True → 0). So totals = 0 for empty months; flag = s.notna() & totals>max → all-NaN month days: s.notna() False → not flagged ✓. Status np.select: condlist evaluated in order; missing takes precedence ✓. Also add a tiny function to convert flagged values to NaN for downstream use? Mention one-liner in notes: `qc.loc[qc["status"] == "fail", "prcp"] = np.nan`. Good. Now, notes section: - thresholds are climate/station specific (arid vs monsoon). - The outlier fence is deliberately generous (k=3 IQR on wet days); consider percentile-based (e.g., >99.9th percentile) or comparing with regional neighbors for extremes. - Suggest iterative QC: run, fix/remove, re-run. - For dense networks: spatial consistency check with nearby stations (correlation / distance-weighted interpolation). Mention briefly with a mini snippet? Maybe just mention. Could add optional short function `check_with_neighbors`? That'd lengthen the answer. Mention as extension, not code. Keep answer manageable. Also mention flag conventions used by e.g. GHCN (M flags) and RClimDex; the boolean columns can be mapped to numeric codes. I think we're good. Also ensure no misuse: `hits[:n_examples].strftime(...)` returns Index of strings; join ✓. `value_counts().reindex([...], fill_value=0)` ✓. One more detail: in `qc_report`, `qc[col]` — qc DataFrame contains prcp + flag cols + n_flags + status; `qc[col]` for col in FLAGS ✓. Also for `examples` join: `hits[:n_examples]` is DatetimeIndex slice ✓; `.strftime("%Y-%m-%d")` → Index[str] ✓. Also `len(qc)` = days count ✓. Potential deprecation: none major. `pd.date_range(..., freq="D")` fine. `np.select` with Series in condlist: it converts to arrays via np.asarray — but index alignment not needed since all same index/order ✓. Also in run_qc: after concat, return. Also duplicated timestamps note: we silently drop duplicates; mention in docstring/comment ("duplicates: keeps first occurrence"). Could add to report — skip for brevity but mention in comment. Edge: `s.index.to_period("M")` requires DatetimeIndex ✓. Now, write the final response: Structure: 1. One-line intro: what the script does + list of checks (concise bullet). 2. Code block (complete, runnable, with demo). 3. Example output (illustrative). 4. Notes: how to use with your data, tuning thresholds, handling flagged values, possible extensions (spatial checks etc.). Keep intro short; code speaks. Use headers. Double-check code once more for typos... Also `int(wet.sum())` fine. In demo, `rain[400:402] = [0.0, 200.0]` — numpy assigns floats ✓. `s.drop(s.index[[500, 501, 800]])` ✓. One more subtle thing: in `_check_repeated`, `s.groupby(run_id)` — groupby key int series aligned ✓; transform("size") returns group sizes as int64 Series aligned to original index ✓. Also in `_check_spell`, `mask.groupby(run_id)` where run_id from cumsum of starts — mask must be boolean Series ✓; transform("sum") on bool → int64 ✓. Wait, `mask.groupby(run_id).transform("sum")` — groupby by a Series (run_id) with same index ✓. Also `transform("sum")` on boolean Series — returns int64 (sum of bools). ✓. `.where(mask, 0)` ✓. Potential FutureWarning: groupby with transform on bool? fine. Also `flags["missing"] = s.isna() | s.isin(...)` — isin with list of floats ✓. `flags["negative"] = s < 0` → NaN < 0 → False ✓. OK. Also should the jump check exclude days flagged as above_limit/negative to avoid double counting? Not necessary; flags independent; status aggregates. One more improvement: for the "repeated" check, RClimDex also flags repeated zeros? No — zeros legit. OK. Consider also adding "trace" note. Also: `astype(float)` before isin with float codes ✓ (int codes matched as floats: -999 == -999.0 True in isin? `pd.Series([-999]).isin([-999.0])` → True, since equality numeric ✓). Write sample output block: ``` Daily rainfall QC | 1096 days (2020-01-01 to 2022-12-31) -------------------------------------------------------------------------- missing 3 e.g. 2021-05-15, 2021-05-16, 2022-03-11 negative 1 e.g. 2020-04-10 above_limit 1 e.g. 2020-07-19 repeated 5 e.g. 2020-10-27, 2020-10-28, 2020-10-29 ... outlier 9 e.g. 2020-08-14, 2020-09-02 ... ← counts vary (random data) wet_spell 37 e.g. 2021-08-23, 2021-08-24 ... dry_spell 196 e.g. 2022-06-19, 2022-06-20 ... jump 4 e.g. 2020-07-19, 2020-07-20, 2021-02-05 ... -------------------------------------------------------------------------- status: ok=841 suspect=245 fail=7 missing=3 ``` Hmm alignment: `{col:<13}{len:>6}` → "missing " + " 3". Fine approximations. I'll include note "numbers for random-data checks (outlier) will vary". Actually to be safe, I'll present the output as "example output (illustrative)" and keep the caveat. Good. Also mention how to view flagged rows: `qc.loc["2020-10-27":"2020-10-31"]`. Now, notes bullets: - Replace thresholds per climate: arid regions lower upper_limit? Actually arid: max daily extreme maybe 100-200mm; humid tropics higher. jump_threshold etc. - Outlier test uses only wet days to avoid zero-inflation problem. - Iterative QC & visual inspection (plot the series, mark flags) recommended before deleting values; usually you don't delete — you flag and let the user decide, or set to NaN for analysis. - Extensions: spatial consistency with neighbor stations (e.g., flag if station value deviates strongly from distance-weighted neighbor estimate), double-mass curve for gauge drift, homogeneity tests, `xclim`/`climpact` for climate indices, GHCN-DQC flags. - If hourly data: adapt. Also mention reading from CSV snippet is in the demo comment. I'll also add an optional tiny plotting snippet? Maybe skip to keep length in check; a one-liner mention: `qc["prcp"].plot(...)` with flags overlay. Could add 4-line optional plot: ```python import matplotlib.pyplot as plt ax = qc["prcp"].plot(style=".", ms=3, label="ok") qc.loc[qc["status"]=="suspect", "prcp"].plot(style="o", color="orange", ls="none", label="suspect", ax=ax) qc.loc[qc["status"]=="fail", "prcp"].plot(style="X", color="red", ls="none", label="fail", ax=ax) ax.legend() ``` Include as optional block at the end of notes? It's nice for QC. I'll add it briefly in the notes section. Length check: this will be a fairly long answer but justified. Keep prose minimal. Let me also reconsider the `run_id` cumsum approach in `_spell_lengths` when mask is all False: run_id all 0; groupby(0) transform sum = 0; where(mask,0) → 0s ✓. And when mask all True: starts [T, F, F...], run_id [1,1,1...]; transform sum = n ✓. Edge in `_check_repeated`: all wet equal values e.g. constant 2.0 series: run_id [1,1,...] (ne False after first, ~wet False) → size n → all flagged if n > 3 ✓. And dry zeros: ~wet True each → own runs size 1 → not flagged ✓. Wait, careful: `~wet | s.ne(s.shift())` for dry days: ~wet True → run_id increments each dry day ✓. For wet equal: ~wet False, ne False → same run ✓. For wet different value: ne True → new run ✓. For missing: ~wet True (NaN > 0 → False) → new run ✓. All good. One more: `s.ne(s.shift())` — pandas Series.ne with NaN: returns True for NaN comparisons? Let me double-check: `pd.Series([1.0, np.nan]).ne(pd.Series([np.nan, 1.0]))` → element-wise: 1.0 vs NaN → True; NaN vs 1.0 → True. And `pd.Series([np.nan]).ne(pd.Series([np.nan]))` → pandas comparison functions treat NaN as unequal → True. Yes, `Series.ne` uses `NotEqual` where NaN != NaN → True. ✓ (Unlike `.eq` which yields False for NaN==NaN.) Good — missing days always increment run id, which we want anyway. Now finalize. Also think about whether to make `run_qc` accept DataFrame with column name — keep Series API; user selects column. Fine. Also mention that flags are boolean columns so users can combine/adjust. I'll write the answer now. Markdown with sections: "What it checks" (bullets), code, "Example output", "Using it with your own data + notes". Keep notes as bullets. One more sanity check on `np.select` default param: `np.select(condlist, choicelist, default='ok')` ✓ returns object array of strings; assigned to column ✓. And `flags["n_flags"] = flags[list(FLAGS)].sum(axis=1)` ✓. Return `pd.concat([s.rename("prcp"), flags], axis=1)` ✓ — index alignment same ✓. Also in report, `hits.size` for DatetimeIndex ✓. Final: also make sure demo `print(qc_report(qc))` and maybe print a peek of flagged rows: ```python print(qc.loc["2020-10-27":"2020-10-31"]) ``` Could include; adds length. I'll include one commented line. Also note about performance: vectorized, handles ~centuries of daily data fine. Alright, compose the final answer. Here's a complete, self-contained QC script for daily rainfall following the standard checks used in climatology (RClimDex/WMO-style): missing data, physical limits, stuck gauges, statistical outliers, wet/dry spells, and day-to-day jumps. ## What it checks | Check | Type | |---|---| | Missing values / missing dates | metadata | | Negative values | hard error | | Value above physically possible maximum | hard error | | Same non-zero value repeated N+ days (stuck gauge) | hard error | | Statistical outlier (Tukey fence on wet days only) | suspect | | Too-long wet spell (consecutive rainy days) | suspect | | Too-long dry spell | suspect | | Unrealistic day-to-day jump | suspect | | (Optional) monthly total above a maximum | suspect | ```python """ Quality control (QC) for daily rainfall series. Each check produces a boolean flag column (True = failed the check); a 'status' column summarises every day: ok / suspect / fail / missing. """ from __future__ import annotations import numpy as np import pandas as pd # --------------------------------------------------------------------------- # # Configuration — tune thresholds to your climate / instrument # # --------------------------------------------------------------------------- # DEFAULTS = dict( missing_codes=(-999.0, -99.0, -9999.0), upper_limit=1000.0, # mm/day; world record ~1825 mm (Foc-Foc, Réunion 1966) max_repeated=3, # flag runs of >3 identical non-zero values iqr_k=3.0, # Tukey "far out" fence: Q3 + k*IQR on wet-day amounts min_wet_days=50, # outlier test needs at least this many wet days max_wet_spell=30, # >30 consecutive wet days -> flag max_dry_spell=180, # >180 consecutive dry days -> flag jump_threshold=150.0, # mm; |change between consecutive days| max_monthly_total=None, # e.g. 5000 mm/month for extreme monsoon sites; None = off ) HARD = ("negative", "above_limit", "repeated") SOFT = ("outlier", "wet_spell", "dry_spell", "jump", "monthly_total") FLAGS = ("missing",) + HARD + SOFT # --------------------------------------------------------------------------- # # Helper checks # # --------------------------------------------------------------------------- # def _spell_lengths(mask: pd.Series) -> pd.Series: """For each day, the length of the run of consecutive True it belongs to.""" run_id = (mask & ~mask.shift(fill_value=False)).cumsum() return mask.groupby(run_id).transform("sum").where(mask, 0) def _check_spell(mask: pd.Series, max_len: int) -> pd.Series: return mask & (_spell_lengths(mask) > max_len) def _check_repeated(s: pd.Series, max_repeated: int) -> pd.Series: """Flag runs of identical non-zero values (zeros are never flagged).""" wet = s > 0 # a new run starts when the value changes, the day is dry, or data is missing run_id = (~wet | s.ne(s.shift())).cumsum() run_len = s.groupby(run_id).transform("size") return wet & (run_len > max_repeated) def _check_outlier(s: pd.Series, k: float, min_wet_days: int) -> pd.Series: """Tukey upper fence applied to wet-day amounts only (dry days are mostly 0, so the fence must not be computed on the full series).""" wet = s.dropna() wet = wet[wet > 0] if wet.size < min_wet_days: return pd.Series(False, index=s.index) q1, q3 = wet.quantile([0.25, 0.75]) iqr = q3 - q1 if iqr <= 0: # all wet days identical -> test not applicable return pd.Series(False, index=s.index) return s > (q3 + k * iqr) def _check_monthly(s: pd.Series, max_monthly: float | None) -> pd.Series: flag = pd.Series(False, index=s.index) if max_monthly is None: return flag totals = s.clip(lower=0).groupby(s.index.to_period("M")).transform("sum") return s.notna() & (totals > max_monthly) # --------------------------------------------------------------------------- # # Main QC routine # # --------------------------------------------------------------------------- # def run_qc(prcp: pd.Series, cfg: dict | None = None) -> pd.DataFrame: """Run all checks on a daily rainfall series (mm). Parameters ---------- prcp : pd.Series with a DatetimeIndex (daily). cfg : overrides for DEFAULTS. Returns a DataFrame with the values, one boolean flag column per check, 'n_flags' and a 'status' column (ok / suspect / fail / missing). """ cfg = {**DEFAULTS, **(cfg or {})} s = prcp.astype(float).copy() s.index = pd.to_datetime(s.index) s = s.sort_index() if s.index.has_duplicates: # keep first occurrence s = s[~s.index.duplicated(keep="first")] if s.empty: raise ValueError("empty series") # reindex to a complete daily calendar: gaps in the record become NaN -> 'missing' s = s.reindex(pd.date_range(s.index.min(), s.index.max(), freq="D")) flags = pd.DataFrame(index=s.index) flags["missing"] = s.isna() | s.isin(list(cfg["missing_codes"])) flags["negative"] = s < 0 flags["above_limit"] = s > cfg["upper_limit"] flags["repeated"] = _check_repeated(s, cfg["max_repeated"]) flags["outlier"] = _check_outlier(s, cfg["iqr_k"], cfg["min_wet_days"]) flags["wet_spell"] = _check_spell(s > 0, cfg["max_wet_spell"]) flags["dry_spell"] = _check_spell(s == 0, cfg["max_dry_spell"]) prev = s.shift(1) flags["jump"] = s.notna() & prev.notna() & \ ((s - prev).abs() > cfg["jump_threshold"]) flags["monthly_total"] = _check_monthly(s, cfg["max_monthly_total"]) flags["n_flags"] = flags[list(FLAGS)].sum(axis=1) flags["status"] = np.select( [flags["missing"], flags[list(HARD)].any(axis=1), flags[list(SOFT)].any(axis=1)], ["missing", "fail", "suspect"], default="ok", ) return pd.concat([s.rename("prcp"), flags], axis=1) def qc_report(qc: pd.DataFrame, n_examples: int = 6) -> str: """Human-readable summary of the QC results.""" lines = [f"Daily rainfall QC | {len(qc)} days " f"({qc.index.min():%Y-%m-%d} to {qc.index.max():%Y-%m-%d})", "-" * 74] for col in FLAGS: hits = qc.index[qc[col]] if hits.size: ex = ", ".join(hits[:n_examples].strftime("%Y-%m-%d")) dots = " ..." if hits.size > n_examples else "" lines.append(f"{col:<13}{len(hits):>6} e.g. {ex}{dots}") counts = qc["status"].value_counts().reindex( ["ok", "suspect", "fail", "missing"], fill_value=0) lines.append("-" * 74) lines.append("status: " + " ".join(f"{k}={v}" for k, v in counts.items())) return "\n".join(lines) # --------------------------------------------------------------------------- # # Demo with synthetic data + injected errors # # --------------------------------------------------------------------------- # if __name__ == "__main__": rng = np.random.default_rng(42) idx = pd.date_range("2020-01-01", "2022-12-31", freq="D") n = idx.size rain = np.zeros(n) wet = rng.random(n) < 0.35 rain[wet] = rng.gamma(0.6, 8.0, int(wet.sum())) rain = np.round(rain, 1) # --- inject typical real-world problems -------------------------------- rain[100] = -2.0 # negative value rain[200] = 5000.0 # impossible magnitude rain[300:305] = 2.0 # stuck sensor (5 identical days) rain[400:402] = [0.0, 200.0] # day-to-day jump rain[600:637] = rng.gamma(0.7, 3.0, 37) # 37-day wet spell rain[900:1100] = 0.0 # long dry spell s = pd.Series(rain, index=idx, name="rain_mm") s = s.drop(s.index[[500, 501, 800]]) # gaps in the record qc = run_qc(s) print(qc_report(qc)) # Real data: # df = pd.read_csv("station.csv", parse_dates=["date"], index_col="date") # qc = run_qc(df["rain_mm"]) ``` ## Example output (illustrative) ``` Daily rainfall QC | 1096 days (2020-01-01 to 2022-12-31) -------------------------------------------------------------------------- missing 3 e.g. 2021-05-15, 2021-05-16, 2022-03-11 negative 1 e.g. 2020-04-10 above_limit 1 e.g. 2020-07-19 repeated 5 e.g. 2020-10-27, 2020-10-28, 2020-10-29 ... outlier 9 e.g. 2020-09-05, 2021-06-11 ... <- varies with data wet_spell 37 e.g. 2021-08-23, 2021-08-24 ... dry_spell 196 e.g. 2022-06-19, 2022-06-20 ... jump 4 e.g. 2020-07-19, 2020-07-20, 2021-02-05 ... -------------------------------------------------------------------------- status: ok=841 suspect=245 fail=7 missing=3 ``` ## Practical notes - **Use with your data**: `qc = run_qc(df["rain_mm"], cfg={"upper_limit": 300, "jump_threshold": 100})` — thresholds are climate-specific (arid vs. monsoon regions differ a lot). - **Flag, don't delete (yet)**: standard practice is to flag values, inspect them (plots, station metadata, neighboring stations), and only then set bad values to NaN, e.g. `qc.loc[qc["status"] == "fail", "prcp"] = np.nan`. - **Iterate**: after removing confirmed errors, re-run the QC since the outlier fence depends on the data. - **Visual inspection** is essential in QC: ```python import matplotlib.pyplot as plt ax = qc["prcp"].plot(style=".", ms=3, label="ok") qc.loc[qc["status"] == "suspect", "prcp"].plot(style="o", color="orange", ls="none", ax=ax, label="suspect") qc.loc[qc["status"] == "fail", "prcp"].plot(style="X", color="red", ls="none", ax=ax, label="fail") ax.set_ylabel("rain (mm/day)"); ax.legend() ``` - **Extensions** if you need them: spatial consistency (compare with distance-weighted values from nearby stations), double-mass curves to detect gauge drift/accumulation errors, trace-value handling, and homogeneity testing. Everything is vectorized with pandas/numpy, so it runs quickly even on multi-decade daily series.

Hydrometeorological Quality Control