"""Trend-detective analysis: did fire activity in each Ukrainian region change after Feb 2022, beyond what weather and the pre-war trend explain? Usage: python3 -I analyze.py Per region and sensor (MODIS, VIIRS kept separate, never summed): model log1p(monthly detections) ~ month-of-year + linear year trend + weather anomalies (MERRA-2 via NASA POWER: T2M_MAX, PRECTOTCORR, RH2M, GWETTOP) fit on pre-war months only. Baseline: 2001-2021 (MODIS) / 2013-2021 (VIIRS); Donetsk, Luhansk, Crimea: 2001-2013 (conflict/occupation from 2014), MODIS only. pre-war residuals are leave-one-year-out (out of sample), like the post-war ones. effect mean annual residual Mar 2022 onward -> ratio observed / expected test Welch t-test post vs pre annual residuals, Benjamini-Hochberg across regions also Mann-Kendall and Pettitt on raw annual totals, to show what a pure trend test says. Writes results.csv, annual.csv (observed vs expected per region-year), national.csv. """ import os import sys import numpy as np import pandas as pd import pymannkendall as mk from scipy import stats WX = ["T2M_MAX", "PRECTOTCORR", "RH2M", "GWETTOP"] EARLY = {"Donetsk", "Luhansk", "Crimea"} POST_START = (2022, 3) def pettitt(x): x = np.asarray(x, float) n = len(x) s = np.sign(x[None, :] - x[:, None]) u = np.array([s[: t + 1, t + 1:].sum() for t in range(n - 1)]) k = np.argmax(np.abs(u)) K = np.abs(u[k]) p = 2 * np.exp(-6 * K ** 2 / (n ** 3 + n ** 2)) return k, min(p, 1.0) def bh(p): p = np.asarray(p, float) ok = ~np.isnan(p) q = np.full_like(p, np.nan) pv = p[ok] order = np.argsort(pv) ranked = pv[order] * len(pv) / (np.arange(len(pv)) + 1) ranked = np.minimum.accumulate(ranked[::-1])[::-1] out = np.empty_like(pv) out[order] = np.minimum(ranked, 1) q[ok] = out return q USE_TREND = os.environ.get("NO_TREND") != "1" METRIC = os.environ.get("METRIC", "n") # "n" = detections, "frp_sum" = summed fire radiative power (MW) # Regions with major ground fighting or occupation in 2022, from public records (not from the fire data). FRONTLINE_2022 = {"Kyiv", "Chernihiv", "Sumy", "Kharkiv", "Donetsk", "Luhansk", "Zaporizhzhia", "Kherson", "Mykolaiv"} def design(df, base_mask): X = [np.ones(len(df))] for m in range(2, 13): X.append((df["month"] == m).astype(float).values) if USE_TREND: X.append((df["year"] - 2010).values / 10.0) for v in WX: X.append(df[v + "_z"].values) return np.column_stack(X) def fit_region(d, base_years): d = d.sort_values(["year", "month"]).reset_index(drop=True) post = (d["year"] > POST_START[0]) | ((d["year"] == POST_START[0]) & (d["month"] >= POST_START[1])) base = d["year"].isin(base_years) & ~post # weather anomalies: z-score per calendar month against the baseline for v in WX: g = d[base].groupby("month")[v] mu, sd = d["month"].map(g.mean()), d["month"].map(g.std()).replace(0, np.nan) d[v + "_z"] = ((d[v] - mu) / sd).fillna(0.0) y = np.log1p(d[METRIC].values) X = design(d, base) pred = np.full(len(d), np.nan) # leave-one-year-out predictions for baseline years for yr in base_years: tr = base & (d["year"] != yr) te = base & (d["year"] == yr) if te.sum() == 0: continue b, *_ = np.linalg.lstsq(X[tr], y[tr], rcond=None) pred[te] = X[te] @ b b, *_ = np.linalg.lstsq(X[base], y[base], rcond=None) pred[post] = X[post] @ b d["pred"] = pred d["resid"] = y - pred d["exp_n"] = np.expm1(np.clip(pred, 0, None)) d["obs_m"] = d[METRIC] d["post"] = post return d[base | post] def main(work): fm = pd.read_csv(os.path.join(work, "fires_monthly.csv")) wx = pd.read_csv(os.path.join(work, "weather_monthly.csv")) df = fm.merge(wx, on=["region", "year", "month"], how="left") rows, annual = [], [] for (region, sensor), d in df.groupby(["region", "sensor"]): yrs = sorted(d["year"].unique()) if region in EARLY: if sensor == "VIIRS": continue base_years = [y for y in yrs if y <= 2013] else: base_years = [y for y in yrs if y <= 2021] if len(base_years) < 8 or d.loc[d.year.isin(base_years), "n"].sum() < 50: continue r = fit_region(d.copy(), base_years) ay = r.groupby(["year", "post"]).agg(resid=("resid", "mean"), obs=("obs_m", "sum"), exp=("exp_n", "sum")).reset_index() pre_res = ay.loc[~ay.post, "resid"].values post_res = ay.loc[ay.post, "resid"].values if len(post_res) == 0: continue t, p = stats.ttest_ind(post_res, pre_res, equal_var=False) tot = d.groupby("year")[METRIC].sum() mkr = mk.original_test(tot.values) k, pp = pettitt(tot.values) rows.append(dict(region=region, sensor=sensor, base=f"{base_years[0]}-{base_years[-1]}", post_years=",".join(str(int(v)) for v in ay.loc[ay.post, "year"]), ratio=float(np.exp(post_res.mean())), p=float(p), pre_sd=float(pre_res.std(ddof=1)), obs_post=float(ay.loc[ay.post, "obs"].sum()), exp_post=float(ay.loc[ay.post, "exp"].sum()), mk_trend=mkr.trend, mk_p=float(mkr.p), sen_slope=float(mkr.slope), pettitt_year=int(tot.index[k + 1]), pettitt_p=float(pp))) ay["region"], ay["sensor"] = region, sensor annual.append(ay) res = pd.DataFrame(rows) for s in res.sensor.unique(): m = res.sensor == s res.loc[m, "q"] = bh(res.loc[m, "p"].values) res["group"] = np.where(res.region.isin(FRONTLINE_2022), "frontline_2022", np.where(res.region == "Crimea", "occupied_2014", "rear")) res = res.sort_values(["sensor", "ratio"], ascending=[True, False]) tag = ("" if USE_TREND else "_notrend") + ("" if METRIC == "n" else "_" + METRIC) res.to_csv(os.path.join(work, f"results{tag}.csv"), index=False) pd.concat(annual).to_csv(os.path.join(work, f"annual{tag}.csv"), index=False) # group test: frontline vs rear log-ratios (regions as units), per sensor print("metric:", METRIC, "| model:", "month + trend + weather" if USE_TREND else "month + weather (no trend)") for s in res.sensor.unique(): a = np.log(res[(res.sensor == s) & (res.group == "frontline_2022")].ratio) b = np.log(res[(res.sensor == s) & (res.group == "rear")].ratio) u = stats.mannwhitneyu(a, b, alternative="two-sided") print(f"{s}: frontline median ratio {np.exp(a.median()):.2f} (n={len(a)}) | " f"rear {np.exp(b.median()):.2f} (n={len(b)}) | contrast x{np.exp(a.median()-b.median()):.2f} " f"| Mann-Whitney p={u.pvalue:.4f}") nat = df.groupby(["sensor", "year"])["n"].sum().unstack(0) nat.to_csv(os.path.join(work, "national.csv")) pd.set_option("display.width", 200) print(res[["region", "sensor", "base", "post_years", "ratio", "p", "q", "obs_post", "exp_post", "mk_trend", "pettitt_year", "pettitt_p"]].round(3).to_string(index=False)) if __name__ == "__main__": main(sys.argv[1])