#!/usr/bin/env python3 """ ============================================================ SIX-CRITERION MAGNITUDE VARIANCE CONSENSUS DETECTION FRAMEWORK Eight-Step Operational Protocol — Complete Implementation ============================================================ Author : Ramakrishna Pasupuleti ORCID : 0009-0008-8418-1430 Email : workisfun415@gmail.com Version : 2.0 (June 2026) GitHub : https://github.com/workisfun415/six-criterion-earthquake-zones Zenodo : https://doi.org/10.5281/zenodo.20613231 Preprint: https://doi.org/10.21203/rs.3.rs-10000129/v1 BSSA : BSSA-D-26-00168 (under review) EIGHT STEPS: Step 1 : Catalog loading and preparation Step 2 : Grid definition and Mc estimation Step 3 : Baseline computation (training period 2000-2010) Step 4 : Six-signal extraction (C, Q, S, b, n90, Cv) Step 5 : Consensus voting (6 scoring functions F1-F6) Step 6 : Zone declaration and forecast generation Step 7 : Weekly traffic-light monitoring Step 8 : Molchan diagram + N-test validation USAGE: python six_criterion_framework.py \ --catalog your_catalog.csv \ --region India \ --lat_min 8 --lat_max 37 \ --lon_min 68 --lon_max 98 \ --radius 200 \ --min_mag 3.0 \ --assess 2026-06-12 \ --output results/ LICENSE: CC BY 4.0 ============================================================ """ import os, sys, argparse, warnings import numpy as np import pandas as pd from scipy.ndimage import uniform_filter1d from scipy.stats import mannwhitneyu, skew as scipy_skew from datetime import datetime, timedelta import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import matplotlib.patches as mpatches import matplotlib.dates as mdates import matplotlib.gridspec as gridspec warnings.filterwarnings('ignore') # ───────────────────────────────────────────────────────────── # CONSTANTS # ───────────────────────────────────────────────────────────── VERSION = '2.0' TRAIN_START = pd.Timestamp('2000-01-01') TRAIN_END = pd.Timestamp('2010-01-01') WINDOW_W = 20 # rolling window size (events) WINDOW_W0 = 40 # minimum events for baseline SMOOTH_K = 12 # smoothing kernel size Q_THRESHOLD = 0.90 # quiescence threshold TOP_N = 30 # top-N zones per scoring function MIN_VOTES = 3 # minimum votes for critical zone # Step 7 thresholds C_GREEN = 0.30 C_ORANGE = 0.50 C_DROP_RED = 0.15 # C collapse without earthquake = RED C_DROP_REL = 0.05 # C after drop = RELEASED EQ_CHECK_DAYS = 60 # days to check for recent earthquake EQ_CHECK_MAG = 4.5 # minimum magnitude for release check # ───────────────────────────────────────────────────────────── # UTILITY FUNCTIONS # ───────────────────────────────────────────────────────────── def haversine(lat1, lon1, lat2, lon2): """Haversine distance in km between two points.""" R = 6371.0 dlat = np.radians(lat2 - lat1) dlon = np.radians(lon2 - lon1) a = (np.sin(dlat/2)**2 + np.cos(np.radians(lat1)) * np.cos(np.radians(lat2)) * np.sin(dlon/2)**2) return R * 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a)) def haversine_vec(zlat, zlon, lats, lons): """Vectorized haversine for array inputs.""" R = 6371.0 dlat = np.radians(lats - zlat) dlon = np.radians(lons - zlon) a = (np.sin(dlat/2)**2 + np.cos(np.radians(zlat)) * np.cos(np.radians(lats)) * np.sin(dlon/2)**2) return R * 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a)) def compute_b_value(mags, Mc=3.0): """Gutenberg-Richter b-value (Aki 1965 MLE).""" m = mags[mags >= Mc] if len(m) < 15: return np.nan mm = float(np.mean(m)) if mm <= Mc: return np.nan return round(float(np.log10(np.e) / (mm - Mc)), 3) def estimate_Mc(mags): """Simple Mc estimation using maximum curvature method.""" if len(mags) < 20: return float(np.min(mags)) counts, bins = np.histogram(mags, bins=np.arange( np.floor(mags.min()*2)/2, np.ceil(mags.max()*2)/2 + 0.5, 0.1)) if len(counts) == 0: return float(np.min(mags)) peak_idx = np.argmax(counts) return float(bins[peak_idx]) def log(msg, level='INFO'): """Simple logger.""" ts = datetime.now().strftime('%H:%M:%S') print(f'[{ts}] {level:5s} | {msg}') # ───────────────────────────────────────────────────────────── # STEP 1: CATALOG LOADING AND PREPARATION # ───────────────────────────────────────────────────────────── def step1_load_catalog(catalog_path, min_mag=3.0, lat_min=None, lat_max=None, lon_min=None, lon_max=None): """ Step 1: Load and prepare earthquake catalog. Accepts CSV with columns (flexible naming): time/date, latitude/lat, longitude/lon, magnitude/mag Returns cleaned DataFrame with standardized columns: time_n (datetime), latitude, longitude, magnitude """ log('STEP 1: Loading catalog...') df = pd.read_csv(catalog_path, low_memory=False) df.columns = [c.lower().strip() for c in df.columns] # Standardize column names col_map = {} for c in df.columns: if c in ('time','date','datetime','origin_time','origintime'): col_map[c] = 'time_raw' elif c in ('latitude','lat','y'): col_map[c] = 'latitude' elif c in ('longitude','lon','long','x'): col_map[c] = 'longitude' elif c in ('magnitude','mag','ml','mw','ms','mb'): col_map[c] = 'magnitude' elif c in ('depth','dep'): col_map[c] = 'depth' df = df.rename(columns=col_map) # Parse time if 'time_raw' in df.columns: df['time_n'] = pd.to_datetime( df['time_raw'], utc=True, errors='coerce' ).dt.tz_localize(None) else: log('No time column found!', 'ERROR') sys.exit(1) # Filter df = df.dropna(subset=['time_n','latitude','longitude','magnitude']) df = df[df['magnitude'] >= min_mag] if lat_min is not None: df = df[(df.latitude >= lat_min) & (df.latitude <= lat_max) & (df.longitude >= lon_min) & (df.longitude <= lon_max)] df = df.sort_values('time_n').reset_index(drop=True) log(f' Loaded {len(df):,} events | ' f'M{df.magnitude.min():.1f}–{df.magnitude.max():.1f} | ' f'{df.time_n.min().date()} to {df.time_n.max().date()}') # Simplified declustering (Gardner-Knopoff style) log(' Applying simplified declustering...') is_mainshock = np.ones(len(df), dtype=bool) big = df[df.magnitude >= 5.0].copy() for _, row in big.iterrows(): t0 = row.time_n # Window: 5 days + 150 km window = df[ (df.time_n > t0) & (df.time_n <= t0 + pd.Timedelta(days=5)) ] for i2, r2 in window.iterrows(): d = haversine(row.latitude, row.longitude, r2.latitude, r2.longitude) if d <= 150 and r2.magnitude < row.magnitude: is_mainshock[i2] = False df_clean = df[is_mainshock].reset_index(drop=True) log(f' After declustering: {len(df_clean):,} events ' f'({len(df)-len(df_clean):,} removed)') return df_clean # ───────────────────────────────────────────────────────────── # STEP 2: GRID DEFINITION # ───────────────────────────────────────────────────────────── def step2_define_grid(lat_min, lat_max, lon_min, lon_max, grid_step=1.0): """ Step 2: Define 1x1 degree grid over study region. Returns list of (lat, lon) grid cell centers. """ log('STEP 2: Defining grid...') lats = np.arange(lat_min, lat_max + grid_step, grid_step) lons = np.arange(lon_min, lon_max + grid_step, grid_step) grid = [(round(la, 1), round(lo, 1)) for la in lats for lo in lons] log(f' Grid: {len(lats)} lat x {len(lons)} lon = ' f'{len(grid)} cells') return grid # ───────────────────────────────────────────────────────────── # STEP 3: BASELINE COMPUTATION # ───────────────────────────────────────────────────────────── def step3_compute_baseline(sub_df, Mc=3.0): """ Step 3: Compute baseline statistics from training period (TRAIN_START to TRAIN_END). Returns: sigma0 : baseline magnitude std rate0 : baseline event rate (events/year) mu_s : mean skewness in training period std_s : std of skewness in training period """ train = sub_df[ (sub_df.time_n >= TRAIN_START) & (sub_df.time_n < TRAIN_END) ] if len(train) < WINDOW_W0: # Fall back to first third of data n = max(WINDOW_W0, len(sub_df)//3) train = sub_df.iloc[:n] mags = train.magnitude.values sigma0 = max(float(np.std(mags)), 0.05) days = max((train.time_n.max() - train.time_n.min()).days, 1) rate0 = max(len(train) / (days / 365.25), 0.1) # Skewness baseline skews = [] for k in range(WINDOW_W, len(train)): w = train.magnitude.values[k-WINDOW_W:k] skews.append(float(scipy_skew(w))) mu_s = float(np.mean(skews)) if skews else 0.0 std_s = float(np.std(skews)) if skews else 1.0 std_s = max(std_s, 0.01) return sigma0, rate0, mu_s, std_s # ───────────────────────────────────────────────────────────── # STEP 4: SIX-SIGNAL EXTRACTION # ───────────────────────────────────────────────────────────── def step4_extract_signals(sub_df, sigma0, rate0, mu_s, std_s, assess_date): """ Step 4: Compute six signals for a single grid cell. Signals: C = magnitude compression = max(0, 1 - sigma_W/sigma0) Q = seismicity quiescence = max(0, 1 - rate_W/rate0) S = skewness anomaly = max(0, (skew-mu_s)/(3*std_s)) Psi= C x Q x (1 + S) [combined signal] b = Gutenberg-Richter b-value n90= event count in last 90 days Cv = std(IET)/mean(IET) [inter-event clustering] Returns dict of current signal values and full C waveform (for Step 7 monitoring). """ mags = sub_df.magnitude.values times = sub_df.time_n.values n = len(sub_df) if n < WINDOW_W0 + WINDOW_W: return None # ── Full C waveform ────────────────────────────────────── C_raw = [] Q_raw = [] S_raw = [] T_raw = [] for k in range(WINDOW_W, n): # C — magnitude compression w_mags = mags[k-WINDOW_W:k] sigma_w = float(np.std(w_mags)) C = max(0.0, 1.0 - sigma_w / sigma0) C_raw.append(C) # Q — quiescence t_now = sub_df.time_n.iloc[k] t_prev = sub_df.time_n.iloc[k-WINDOW_W] dt_days = max((t_now - t_prev).days, 1) rate_w = WINDOW_W / (dt_days / 365.25) Q = max(0.0, 1.0 - rate_w / rate0) Q_raw.append(Q) # S — skewness anomaly sk = float(scipy_skew(w_mags)) S = max(0.0, (sk - mu_s) / (3 * std_s)) S_raw.append(S) T_raw.append(sub_df.time_n.iloc[k]) C_arr = np.array(C_raw) Q_arr = np.array(Q_raw) S_arr = np.array(S_raw) T_arr = pd.DatetimeIndex(T_raw) # Smooth C_sm = uniform_filter1d(C_arr, size=SMOOTH_K) Q_sm = uniform_filter1d(Q_arr, size=SMOOTH_K) S_sm = uniform_filter1d(S_arr, size=SMOOTH_K) # ── Current values at assessment date ──────────────────── def get_at(arr, t): d = np.abs(T_arr - pd.Timestamp(t)) if len(d) == 0: return 0.0 i = int(np.argmin(d)) return float(arr[i]) if i < len(arr) else 0.0 C_now = get_at(C_sm, assess_date) Q_now = get_at(Q_sm, assess_date) S_now = get_at(S_sm, assess_date) # Psi — combined signal Psi = C_now * Q_now * (1 + S_now) # Trend (90-day change) C_90d = get_at(C_sm, pd.Timestamp(assess_date) - pd.Timedelta(days=90)) trend = round(C_now - C_90d, 3) # Peak C in last 12 months mask12 = T_arr >= pd.Timestamp(assess_date) - pd.Timedelta(days=365) if mask12.sum() > 0: pk_C = float(np.max(C_sm[mask12])) pk_T = T_arr[mask12][np.argmax(C_sm[mask12])] else: pk_C = C_now pk_T = pd.Timestamp(assess_date) c_drop = round(pk_C - C_now, 3) # b-value (recent 1 year) cut1yr = pd.Timestamp(assess_date) - pd.Timedelta(days=365) mags1yr = sub_df[sub_df.time_n >= cut1yr].magnitude.values b = compute_b_value(mags1yr) if np.isnan(b): b = compute_b_value(mags) # n90 — event count last 90 days cut90 = pd.Timestamp(assess_date) - pd.Timedelta(days=90) n90 = int(len(sub_df[sub_df.time_n >= cut90])) # Cv — inter-event time coefficient of variation iet = np.diff([t.timestamp() for t in sub_df.time_n.values[-WINDOW_W:]]) iet = iet[iet > 0] / 86400.0 # convert to days if len(iet) > 2: Cv = float(np.std(iet) / max(np.mean(iet), 0.001)) else: Cv = 1.0 return { 'C': round(C_now, 3), 'Q': round(Q_now, 3), 'S': round(S_now, 3), 'Psi': round(Psi, 3), 'b': round(b, 3) if not np.isnan(b) else 1.0, 'n90': n90, 'Cv': round(Cv, 3), 'trend': trend, 'pk_C': round(pk_C, 3), 'pk_T': str(pk_T.date()) if pk_T else '---', 'c_drop':c_drop, 'C_sm': C_sm, 'Q_sm': Q_sm, 'T_arr': T_arr, 'n_events': n, } # ───────────────────────────────────────────────────────────── # STEP 5: CONSENSUS VOTING # ───────────────────────────────────────────────────────────── def step5_consensus_voting(zone_signals): """ Step 5: Apply six scoring functions and consensus vote. Scoring functions: F1 = Psi + n90/100 F2 = n90/100 + Cv F3 = |dC| + |dS| F4 = 1/b F5 = Psi + n90/100 + 1/b + Cv [MEGA] F6 = n90/100 + 1/b Zone declared critical if top-30 in >= MIN_VOTES functions. Returns zone_signals dict with added 'votes' and 'critical'. """ log('STEP 5: Consensus voting...') keys = list(zone_signals.keys()) sigs = [zone_signals[k] for k in keys] def safe_b(s): b = s.get('b', 1.0) return 1.0/max(b, 0.3) if not np.isnan(b) else 1.0 # Compute scores F1 = [s['Psi'] + s['n90']/100 for s in sigs] F2 = [s['n90']/100 + s['Cv'] for s in sigs] F3 = [abs(s['c_drop']) + abs(s['S']) for s in sigs] F4 = [safe_b(s) for s in sigs] F5 = [s['Psi'] + s['n90']/100 + safe_b(s) + s['Cv'] for s in sigs] F6 = [s['n90']/100 + safe_b(s) for s in sigs] all_F = [F1, F2, F3, F4, F5, F6] names = ['F1','F2','F3','F4','F5','F6'] # Rank each function — get top-N indices votes = {k: 0 for k in keys} for Fi, fname in zip(all_F, names): ranked = np.argsort(Fi)[::-1] for idx in ranked[:TOP_N]: votes[keys[idx]] += 1 # Declare critical zones n_critical = 0 for k in keys: zone_signals[k]['votes'] = votes[k] zone_signals[k]['critical'] = votes[k] >= MIN_VOTES if votes[k] >= MIN_VOTES: n_critical += 1 log(f' Critical zones: {n_critical} of {len(keys)} ' f'({n_critical/max(len(keys),1)*100:.1f}% of grid)') return zone_signals # ───────────────────────────────────────────────────────────── # STEP 6: ZONE DECLARATION AND FORECAST # ───────────────────────────────────────────────────────────── def step6_declare_zones(zone_signals, assess_date, region_name, radius, output_dir): """ Step 6: Declare critical zones and generate forecast list. Computes: - Expected earthquake window (based on mean C lead time) - Magnitude estimate from b-value - Saves CSV forecast file Returns sorted list of critical zones. """ log('STEP 6: Declaring critical zones...') assess = pd.Timestamp(assess_date) # Mean lead time from historical analysis MEAN_LEAD = {'Japan': 90, 'Turkey': 85, 'Italy': 110, 'India': 115} lead = MEAN_LEAD.get(region_name, 100) critical = [] for (zlat, zlon), sig in zone_signals.items(): if not sig.get('critical', False): continue # Expected window pk_t = pd.Timestamp(sig['pk_T']) \ if sig['pk_T'] != '---' \ else assess eq_start = pk_t + pd.Timedelta(days=lead//2) eq_end = pk_t + pd.Timedelta(days=lead*2) window = f'{eq_start.date()} to {eq_end.date()}' # Magnitude estimate from b-value b = sig['b'] if b < 0.65: mag_est = 'M5.5-7.0' elif b < 0.80: mag_est = 'M5.0-6.5' elif b < 1.00: mag_est = 'M4.5-6.0' else: mag_est = 'M4.0-5.5' critical.append({ 'lat': zlat, 'lon': zlon, 'b': sig['b'], 'C': sig['C'], 'Q': sig['Q'], 'S': sig['S'], 'Psi': sig['Psi'], 'n90': sig['n90'], 'Cv': sig['Cv'], 'trend': sig['trend'], 'pk_C': sig['pk_C'], 'pk_T': sig['pk_T'], 'c_drop': sig['c_drop'], 'votes': sig['votes'], 'window': window, 'mag_est': mag_est, 'region': region_name, 'radius': radius, 'assess': str(assess_date), }) # Sort by C descending critical.sort(key=lambda x: -x['C']) # Save CSV os.makedirs(output_dir, exist_ok=True) csv_path = os.path.join( output_dir, f'critical_zones_{region_name}_{assess_date}.csv' ) pd.DataFrame(critical).to_csv(csv_path, index=False) log(f' Saved {len(critical)} critical zones → {csv_path}') # Print summary print() print('='*70) print(f'CRITICAL ZONES — {region_name} — {assess_date}') print('='*70) print(f"{'#':>3} {'Lat':>6} {'Lon':>6} {'b':>6} " f"{'C':>6} {'Trend':>7} {'Votes':>6} " f"{'Window':<35} {'Mag'}") print('-'*70) for i, z in enumerate(critical[:30]): arr = '↑' if z['trend'] > 0.03 \ else '↓' if z['trend'] < -0.03 else '→' print(f"{i+1:>3} {z['lat']:>5.1f}N {z['lon']:>5.1f}E " f"{z['b']:>6.3f} {z['C']:>6.3f} " f"{arr}{abs(z['trend']):>5.3f} {z['votes']:>6} " f"{z['window']:<35} {z['mag_est']}") print() return critical # ───────────────────────────────────────────────────────────── # STEP 7: WEEKLY TRAFFIC-LIGHT MONITORING # ───────────────────────────────────────────────────────────── def step7_traffic_light(zone_signals, critical_zones, catalog_df, assess_date, output_dir): """ Step 7: Weekly traffic-light monitoring. For each critical zone assigns: GREEN : C < C_GREEN ORANGE : C = C_GREEN to C_ORANGE, rising RED : C > C_ORANGE OR C dropped >0.15 without earthquake RELEASED : C dropped >0.15 and M>=4.5 confirmed in 60d CRITICAL: Always checks for recent earthquakes before classifying C-collapse as RED. Returns zones with 'flag' and 'timing' added. """ log('STEP 7: Traffic-light monitoring...') assess = pd.Timestamp(assess_date) cut60 = assess - pd.Timedelta(days=EQ_CHECK_DAYS) flagged = [] counts = {'RED':0,'ORANGE':0,'GREEN':0, 'RELEASED':0,'STRESS':0} for z in critical_zones: zlat = z['lat'] zlon = z['lon'] key = (zlat, zlon) sig = zone_signals.get(key, {}) C = z['C'] trend = z['trend'] c_drop= z['c_drop'] b = z['b'] # ── Check for recent earthquake ────────────────────── dists = haversine_vec( zlat, zlon, catalog_df.latitude.values, catalog_df.longitude.values ) zone_cat = catalog_df[dists <= z['radius']].copy() recent_big = zone_cat[ (zone_cat.time_n >= cut60) & (zone_cat.magnitude >= EQ_CHECK_MAG) ] has_recent_eq = len(recent_big) > 0 max_recent_mag = (float(recent_big.magnitude.max()) if has_recent_eq else 0.0) recent_eq_date = (str(recent_big.sort_values( 'magnitude', ascending=False )['time_n'].iloc[0].date()) if has_recent_eq else '---') # ── C drop caused by recent earthquake? ───────────── eq_caused_drop = (has_recent_eq and c_drop > C_DROP_RED and C < C_DROP_REL * 3) # ── Assign flag ────────────────────────────────────── if eq_caused_drop: flag = 'RELEASED 🔵' timing = (f'M{max_recent_mag:.1f} on ' f'{recent_eq_date} released stress') urgency = 0 elif (C > C_ORANGE or (c_drop > C_DROP_RED and C < 0.10 and not has_recent_eq)): flag = 'RED 🔴' if trend > 0.02: d2p = max(0, int((0.55-C)/(trend/90))) else: d2p = 0 peak_est = assess + pd.Timedelta(days=d2p) w1 = peak_est + pd.Timedelta(days=20) w2 = peak_est + pd.Timedelta(days=150) timing = f'{w1.date()} to {w2.date()}' urgency = C * 10 + (1/max(b,0.3)) * 2 elif (C >= C_GREEN and trend > 0.03 and not eq_caused_drop): flag = 'ORANGE 🟠' d2p = max(0, int((0.40-C)/(trend/90))) \ if trend > 0.02 else 90 peak_est = assess + pd.Timedelta(days=d2p) w1 = peak_est + pd.Timedelta(days=60) w2 = peak_est + pd.Timedelta(days=200) timing = f'{w1.date()} to {w2.date()}' urgency = C * 5 + trend * 3 elif C >= C_GREEN and not eq_caused_drop: flag = 'ORANGE 🟠' w1 = assess + pd.Timedelta(days=90) w2 = assess + pd.Timedelta(days=270) timing = f'{w1.date()} to {w2.date()}' urgency = C * 4 elif b < 0.70 and C < 0.10: flag = 'STRESS ⚠️' timing = 'Extreme stress — escalate when C rises' urgency = 1/max(b,0.3) else: flag = 'GREEN 🟢' timing = 'Normal monitoring' urgency = 0 zz = dict(z) zz['flag'] = flag zz['timing'] = timing zz['urgency'] = round(urgency, 2) zz['has_recent_eq'] = has_recent_eq zz['max_recent_mag'] = max_recent_mag zz['recent_eq_date'] = recent_eq_date zz['eq_caused_drop'] = eq_caused_drop flagged.append(zz) # Count for k in counts: if k in flag: counts[k] += 1 break # Sort by urgency flagged.sort(key=lambda x: -x['urgency']) # Print Step 7 status print('='*70) print(f'STEP 7 MONITORING STATUS — {assess_date}') print('='*70) print(f" 🔴 RED : {counts['RED']}") print(f" 🟠 ORANGE : {counts['ORANGE']}") print(f" 🟢 GREEN : {counts['GREEN']}") print(f" 🔵 RELEASED : {counts['RELEASED']}") print(f" ⚠️ STRESS : {counts['STRESS']}") print() # Detailed RED zones reds = [z for z in flagged if 'RED' in z['flag']] if reds: print('RED ZONE DETAIL:') for z in reds: print(f" {z['lat']:.1f}N {z['lon']:.1f}E " f"b={z['b']:.3f} C={z['C']:.3f} " f"trend={z['trend']:+.3f}") print(f" Window: {z['timing']}") print(f" Mag est: {z['mag_est']}") if z['has_recent_eq']: print(f" ⚠️ Recent M{z['max_recent_mag']:.1f} " f"on {z['recent_eq_date']} — " f"{'CAUSED DROP' if z['eq_caused_drop'] else 'check'}") print() # Save Step 7 CSV csv7 = os.path.join( output_dir, f'step7_monitoring_{assess_date}.csv' ) pd.DataFrame(flagged).to_csv(csv7, index=False) log(f' Step 7 results saved → {csv7}') return flagged # ───────────────────────────────────────────────────────────── # STEP 8: VALIDATION — MOLCHAN + N-TEST # ───────────────────────────────────────────────────────────── def step8_validation(catalog_df, zone_signals, assess_date, radius, target_mag=5.0, forecast_days=365, output_dir='.'): """ Step 8: Molchan diagram + N-test validation. Tests whether zone votes predict future earthquakes better than random allocation. Returns validation metrics dict. """ log('STEP 8: Validation (Molchan + N-test)...') assess = pd.Timestamp(assess_date) cutoff = assess + pd.Timedelta(days=forecast_days) # Future earthquakes future = catalog_df[ (catalog_df.time_n > assess) & (catalog_df.time_n <= cutoff) & (catalog_df.magnitude >= target_mag) ].copy() n_future = len(future) if n_future == 0: log(' No future earthquakes found for validation.') return {} log(f' Future M>={target_mag} events: {n_future}') # All zones and their votes all_zones = list(zone_signals.keys()) n_total = len(all_zones) # For each threshold (top-5, 10, 15, 20, 25, 30) results = [] for top_n in [5, 10, 15, 20, 25, TOP_N]: # Select top-N zones by votes sorted_zones = sorted( all_zones, key=lambda k: -zone_signals[k].get('votes', 0) ) top_zones = sorted_zones[:top_n] # Area fraction tau = top_n / max(n_total, 1) # Detection: did future eq fall in any top zone? detected = 0 for _, eq in future.iterrows(): for (zlat, zlon) in top_zones: d = haversine(zlat, zlon, eq.latitude, eq.longitude) if d <= radius: detected += 1 break nu = 1 - (detected / max(n_future, 1)) # miss rate skill = 1 - nu / max(1 - tau, 0.001) det_pct= round(detected/max(n_future,1)*100, 1) results.append({ 'top_n': top_n, 'tau': round(tau, 3), 'nu': round(nu, 3), 'skill': round(skill, 3), 'detected': detected, 'n_future': n_future, 'det_pct': det_pct, }) log(f' Top-{top_n:2d}: tau={tau:.2f} ' f'nu={nu:.2f} skill={skill:.3f} ' f'detected={detected}/{n_future} ({det_pct}%)') # Mann-Whitney test: C before earthquakes vs background # C in 90d before future M>=target_mag earthquakes pre_C = [] back_C = [] for _, eq in future.iterrows(): eq_t = eq.time_n cut90 = eq_t - pd.Timedelta(days=90) # Find nearest zone best_d = 9999 best_z = None for (zlat,zlon) in all_zones: d = haversine(zlat,zlon,eq.latitude,eq.longitude) if d < best_d: best_d = d; best_z = (zlat,zlon) if best_z and best_d <= radius*1.5: sig = zone_signals[best_z] T_arr = sig.get('T_arr', pd.DatetimeIndex([])) C_sm = sig.get('C_sm', np.array([])) if len(T_arr) > 0 and len(C_sm) > 0: mask = (T_arr >= cut90) & (T_arr < eq_t) if mask.sum() > 0: pre_C.append(float(np.mean(C_sm[mask]))) mask_b = (T_arr >= eq_t - pd.Timedelta(days=365)) & \ (T_arr < cut90) if mask_b.sum() > 0: back_C.append(float(np.mean(C_sm[mask_b]))) mw_p = np.nan if len(pre_C) >= 5 and len(back_C) >= 5: _, mw_p = mannwhitneyu(pre_C, back_C, alternative='greater') log(f' Mann-Whitney p = {mw_p:.4f} ' f'(C elevated before earthquakes)') # Save validation val_df = pd.DataFrame(results) csv8 = os.path.join(output_dir, f'step8_validation_{assess_date}.csv') val_df.to_csv(csv8, index=False) log(f' Validation saved → {csv8}') best = max(results, key=lambda x: x['skill']) return { 'molchan_skill': best['skill'], 'detection_rate': best['det_pct'], 'mann_whitney_p': round(mw_p, 4) if not np.isnan(mw_p) else np.nan, 'n_future': n_future, 'results': results, } # ───────────────────────────────────────────────────────────── # FIGURE GENERATION # ───────────────────────────────────────────────────────────── def generate_figures(zone_signals, critical_zones, flagged_zones, val_metrics, catalog_df, assess_date, region_name, output_dir): """Generate all publication-quality figures.""" log('Generating figures...') os.makedirs(output_dir, exist_ok=True) # ── Figure 1: Critical zones map ───────────────────────── fig, ax = plt.subplots(figsize=(14, 10)) fig.patch.set_facecolor('#0A0A1A') ax.set_facecolor('#0A0A1A') # Background seismicity recent = catalog_df[ catalog_df.time_n >= pd.Timestamp('2020-01-01') ] ax.scatter(recent.longitude, recent.latitude, c='#222244', s=recent.magnitude**2*1.5, alpha=0.3, zorder=1) # Color by flag col_map = { 'RED': '#CC0000', 'ORANGE': '#FF8C00', 'STRESS': '#8B0080', 'RELEASED': '#0066CC', 'GREEN': '#228B22', } sz_map = {'RED':500,'ORANGE':350,'STRESS':280, 'RELEASED':200,'GREEN':150} for z in flagged_zones: flag = z.get('flag','GREEN') col = next((c for k,c in col_map.items() if k in flag), '#228B22') sz = next((s for k,s in sz_map.items() if k in flag), 150) ax.scatter([z['lon']], [z['lat']], s=sz, c=col, marker='s', edgecolors='white', lw=0.8, alpha=0.9, zorder=4) ax.text(z['lon']+0.15, z['lat']+0.12, f"C={z['C']:.2f}", fontsize=7, color='white', zorder=5) # Major earthquakes big = catalog_df[ (catalog_df.magnitude >= 5.5) & (catalog_df.time_n >= pd.Timestamp('2020-01-01')) ] ax.scatter(big.longitude, big.latitude, s=big.magnitude**3*2, c='yellow', marker='*', alpha=0.7, zorder=6, edgecolors='white', lw=0.5) ax.set_xlabel('Longitude (°E)', fontsize=11, color='white') ax.set_ylabel('Latitude (°N)', fontsize=11, color='white') ax.tick_params(colors='#AAAAAA') for sp in ax.spines.values(): sp.set_color('#333355') ax.grid(True, alpha=0.15, color='#334') ax.set_title( f'{region_name} Critical Zones | {assess_date}\n' f'Six-Criterion Framework | Ramakrishna Pasupuleti', fontsize=12, fontweight='bold', color='white') patches = [ mpatches.Patch(color='#CC0000', label='RED 🔴'), mpatches.Patch(color='#FF8C00', label='ORANGE 🟠'), mpatches.Patch(color='#8B0080', label='STRESS ⚠️'), mpatches.Patch(color='#0066CC', label='RELEASED 🔵'), mpatches.Patch(color='#228B22', label='GREEN 🟢'), plt.Line2D([0],[0], marker='*', color='w', markerfacecolor='yellow', markersize=12, label='M≥5.5 (2020-)'), ] ax.legend(handles=patches, loc='lower left', fontsize=9, facecolor='#111122', edgecolor='#333355', labelcolor='white') fig_path = os.path.join( output_dir, f'fig_critical_zones_{region_name}.png') plt.savefig(fig_path, dpi=300, bbox_inches='tight', facecolor='#0A0A1A') plt.close() log(f' Zones map saved → {fig_path}') # ── Figure 2: Step 7 C waveforms ──────────────────────── top_zones = [z for z in flagged_zones if 'RED' in z.get('flag','') or 'ORANGE' in z.get('flag','')][:4] if top_zones: fig2, axes = plt.subplots( len(top_zones), 1, figsize=(16, 5*len(top_zones))) if len(top_zones) == 1: axes = [axes] fig2.patch.set_facecolor('white') fig2.suptitle( f'Step 7 Monitoring — C Signal Waveforms | ' f'{region_name} | {assess_date}', fontsize=13, fontweight='bold') for ax2, z in zip(axes, top_zones): key = (z['lat'], z['lon']) sig = zone_signals.get(key, {}) C_sm = sig.get('C_sm', np.array([])) T_arr = sig.get('T_arr', pd.DatetimeIndex([])) if len(C_sm) == 0: continue flag = z.get('flag','GREEN') col = ('#CC0000' if 'RED' in flag else '#FF8C00' if 'ORANGE' in flag else '#228B22') ax2.fill_between(T_arr, 0, C_sm, where=C_sm >= C_GREEN, alpha=0.35, color=col) ax2.fill_between(T_arr, 0, C_sm, where=C_sm < C_GREEN, alpha=0.10, color='gray') ax2.plot(T_arr, C_sm, color=col, lw=2.5) ax2.axhline(C_ORANGE, color='red', lw=1.5, ls='--', alpha=0.8, label='RED threshold') ax2.axhline(C_GREEN, color='orange', lw=1.2, ls=':', alpha=0.8, label='ORANGE threshold') ax2.scatter([T_arr[-1]], [C_sm[-1]], s=200, c='gold', marker='D', zorder=7, edgecolors='k', lw=1.0, label=f'Now C={z["C"]:.3f}') ax2.set_title( f'{z["lat"]:.1f}N {z["lon"]:.1f}E | ' f'b={z["b"]:.3f} | {flag} | ' f'Expected: {z["timing"]}', fontsize=10, fontweight='bold', color=col, loc='left') ax2.set_ylabel('Compression C', fontsize=10) ax2.set_ylim(-0.03, 1.05) ax2.grid(True, alpha=0.2) ax2.xaxis.set_major_formatter( mdates.DateFormatter('%Y')) ax2.legend(fontsize=8, loc='upper right') plt.tight_layout() fig2_path = os.path.join( output_dir, f'fig_step7_waveforms_{region_name}.png') plt.savefig(fig2_path, dpi=300, bbox_inches='tight') plt.close() log(f' Step 7 waveforms saved → {fig2_path}') # ── Figure 3: Molchan diagram ───────────────────────── if val_metrics and 'results' in val_metrics: fig3, ax3 = plt.subplots(figsize=(8, 8)) res = val_metrics['results'] tau = [r['tau'] for r in res] nu = [r['nu'] for r in res] x = np.linspace(0, 1, 100) ax3.plot(x, 1-x, 'k--', lw=1.5, alpha=0.6, label='Random diagonal') ax3.fill_between(x, 1-x, 1, alpha=0.05, color='gray') ax3.plot(tau, nu, 'o-', color='#1F3864', lw=2.5, ms=8, label='Six-Criterion Framework') best = min(res, key=lambda r: r['nu']) ax3.scatter([best['tau']], [best['nu']], s=300, marker='*', c='gold', zorder=5, edgecolors='k', lw=1.0, label=f'Best: skill={val_metrics["molchan_skill"]}') ax3.set_xlabel('Area fraction τ', fontsize=12) ax3.set_ylabel('Miss rate ν', fontsize=12) ax3.set_title( f'Molchan Diagram | {region_name}\n' f'Molchan skill = ' f'{val_metrics["molchan_skill"]:.3f} | ' f'Detection = {val_metrics["detection_rate"]}%', fontsize=12, fontweight='bold') ax3.legend(fontsize=10) ax3.set_xlim(-0.01, 0.6) ax3.set_ylim(-0.02, 1.02) ax3.grid(True, alpha=0.25) fig3_path = os.path.join( output_dir, f'fig_step8_molchan_{region_name}.png') plt.savefig(fig3_path, dpi=300, bbox_inches='tight') plt.close() log(f' Molchan diagram saved → {fig3_path}') log('All figures generated.') # ───────────────────────────────────────────────────────────── # MAIN ORCHESTRATOR — ALL 8 STEPS # ───────────────────────────────────────────────────────────── def run_eight_steps(catalog_path, region_name='Region', lat_min=8, lat_max=37, lon_min=68, lon_max=98, radius=200, min_mag=3.0, assess_date='2026-06-12', grid_step=1.0, output_dir='results', validate=True): """ Run all eight steps of the Six-Criterion Framework. Parameters: ----------- catalog_path : str — path to earthquake catalog CSV region_name : str — name of study region lat_min/max : float — latitude bounds lon_min/max : float — longitude bounds radius : float — detection radius in km min_mag : float — minimum magnitude assess_date : str — assessment date YYYY-MM-DD grid_step : float — grid cell size in degrees (default 1.0) output_dir : str — output directory validate : bool — run Step 8 validation Returns: -------- flagged_zones : list of dicts with all signal values and flags val_metrics : dict of validation metrics (if validate=True) """ print() print('='*70) print(' SIX-CRITERION FRAMEWORK — EIGHT-STEP PROTOCOL') print(f' Region: {region_name} | Assessment: {assess_date}') print(f' Version: {VERSION} | Ramakrishna Pasupuleti') print('='*70) print() # ── STEP 1 ─────────────────────────────────────────────── catalog_df = step1_load_catalog( catalog_path, min_mag, lat_min, lat_max, lon_min, lon_max ) # ── STEP 2 ─────────────────────────────────────────────── grid = step2_define_grid( lat_min, lat_max, lon_min, lon_max, grid_step ) # ── STEPS 3-4: Per-cell signal extraction ──────────────── log('STEPS 3-4: Baseline + Signal extraction...') cat_lats = catalog_df.latitude.values cat_lons = catalog_df.longitude.values zone_signals = {} n_total = len(grid) n_processed = 0 for i, (zlat, zlon) in enumerate(grid): # Get events in zone dists = haversine_vec(zlat, zlon, cat_lats, cat_lons) idx = dists <= radius sub = catalog_df[idx].reset_index(drop=True) if len(sub) < WINDOW_W0 + WINDOW_W: continue # insufficient data # Step 3: Baseline sigma0, rate0, mu_s, std_s = step3_compute_baseline(sub) # Step 4: Signals sig = step4_extract_signals( sub, sigma0, rate0, mu_s, std_s, assess_date ) if sig is None: continue zone_signals[(zlat, zlon)] = sig n_processed += 1 if (i+1) % 50 == 0 or i+1 == n_total: log(f' Processed {i+1}/{n_total} cells ' f'({n_processed} with data)') log(f' Zones with sufficient data: {n_processed}') # ── STEP 5 ─────────────────────────────────────────────── zone_signals = step5_consensus_voting(zone_signals) # ── STEP 6 ─────────────────────────────────────────────── critical_zones = step6_declare_zones( zone_signals, assess_date, region_name, radius, output_dir ) if not critical_zones: log('No critical zones found.', 'WARN') return [], {} # ── STEP 7 ─────────────────────────────────────────────── flagged_zones = step7_traffic_light( zone_signals, critical_zones, catalog_df, assess_date, output_dir ) # ── STEP 8 ─────────────────────────────────────────────── val_metrics = {} if validate: val_metrics = step8_validation( catalog_df, zone_signals, assess_date, radius, target_mag=5.0, forecast_days=365, output_dir=output_dir ) # ── FIGURES ────────────────────────────────────────────── generate_figures( zone_signals, critical_zones, flagged_zones, val_metrics, catalog_df, assess_date, region_name, output_dir ) # ── FINAL SUMMARY ───────────────────────────────────────── print() print('='*70) print('FINAL SUMMARY') print('='*70) print(f' Region : {region_name}') print(f' Assessment : {assess_date}') print(f' Total cells : {n_total}') print(f' With data : {n_processed}') print(f' Critical : {len(critical_zones)}') red = sum(1 for z in flagged_zones if 'RED' in z.get('flag','')) org = sum(1 for z in flagged_zones if 'ORANGE' in z.get('flag','')) rel = sum(1 for z in flagged_zones if 'RELEASED' in z.get('flag','')) st = sum(1 for z in flagged_zones if 'STRESS' in z.get('flag','')) grn = sum(1 for z in flagged_zones if 'GREEN' in z.get('flag','')) print(f' 🔴 RED : {red}') print(f' 🟠 ORANGE : {org}') print(f' 🟢 GREEN : {grn}') print(f' 🔵 RELEASED : {rel}') print(f' ⚠️ STRESS : {st}') if val_metrics: print(f' Skill score : {val_metrics.get("molchan_skill","N/A")}') print(f' Detection : {val_metrics.get("detection_rate","N/A")}%') print(f' Mann-Whitney : p={val_metrics.get("mann_whitney_p","N/A")}') print(f' Output dir : {output_dir}') print('='*70) print() return flagged_zones, val_metrics # ───────────────────────────────────────────────────────────── # COMMAND-LINE INTERFACE # ───────────────────────────────────────────────────────────── def main(): parser = argparse.ArgumentParser( description='Six-Criterion Framework — Eight-Step Protocol', formatter_class=argparse.RawDescriptionHelpFormatter, epilog=""" EXAMPLES: # India: python six_criterion_framework.py \\ --catalog india_ISC_M3plus.csv \\ --region India \\ --lat_min 8 --lat_max 37 \\ --lon_min 68 --lon_max 98 \\ --radius 200 --min_mag 3.0 \\ --assess 2026-06-12 \\ --output results/india/ # Japan: python six_criterion_framework.py \\ --catalog jma_M3plus_2000_2023.csv \\ --region Japan \\ --lat_min 24 --lat_max 48 \\ --lon_min 122 --lon_max 150 \\ --radius 200 --min_mag 3.0 \\ --assess 2023-12-31 \\ --output results/japan/ # Turkey: python six_criterion_framework.py \\ --catalog turkey_AFAD_M3plus.csv \\ --region Turkey \\ --lat_min 36 --lat_max 42.5 \\ --lon_min 26 --lon_max 45 \\ --radius 100 --min_mag 3.0 \\ --assess 2026-05-28 \\ --output results/turkey/ # Italy: python six_criterion_framework.py \\ --catalog italy_ISC_M2plus.csv \\ --region Italy \\ --lat_min 36 --lat_max 48 \\ --lon_min 6 --lon_max 20 \\ --radius 150 --min_mag 2.0 \\ --assess 2026-06-12 \\ --output results/italy/ CATALOG FORMAT (CSV): Required columns: time, latitude, longitude, magnitude Time format: ISO 8601 (YYYY-MM-DD HH:MM:SS or similar) Accepts: JMA, ISC, AFAD, USGS formats automatically OUTPUT FILES: critical_zones_{region}_{date}.csv — zone list step7_monitoring_{date}.csv — traffic-light status step8_validation_{date}.csv — Molchan results fig_critical_zones_{region}.png — zones map fig_step7_waveforms_{region}.png — C waveforms fig_step8_molchan_{region}.png — Molchan diagram """ ) parser.add_argument('--catalog', required=True, help='Path to earthquake catalog CSV') parser.add_argument('--region', default='Region', help='Region name (Japan/Turkey/Italy/India)') parser.add_argument('--lat_min', type=float, required=True) parser.add_argument('--lat_max', type=float, required=True) parser.add_argument('--lon_min', type=float, required=True) parser.add_argument('--lon_max', type=float, required=True) parser.add_argument('--radius', type=float, default=200, help='Detection radius in km (default 200)') parser.add_argument('--min_mag', type=float, default=3.0, help='Minimum magnitude (default 3.0)') parser.add_argument('--assess', default='2026-06-12', help='Assessment date YYYY-MM-DD') parser.add_argument('--grid_step',type=float, default=1.0, help='Grid cell size in degrees (default 1.0)') parser.add_argument('--output', default='results', help='Output directory (default: results)') parser.add_argument('--no_validate', action='store_true', help='Skip Step 8 validation') args = parser.parse_args() flagged, metrics = run_eight_steps( catalog_path = args.catalog, region_name = args.region, lat_min = args.lat_min, lat_max = args.lat_max, lon_min = args.lon_min, lon_max = args.lon_max, radius = args.radius, min_mag = args.min_mag, assess_date = args.assess, grid_step = args.grid_step, output_dir = args.output, validate = not args.no_validate, ) return 0 # ───────────────────────────────────────────────────────────── # QUICK-START: Run directly with existing catalogs # ───────────────────────────────────────────────────────────── def run_india_quick(): """Quick-start for India using existing catalog.""" return run_eight_steps( catalog_path='/mnt/user-data/outputs/india_ISC_M3plus_2000_2026.csv', region_name ='India', lat_min=8, lat_max=37, lon_min=68, lon_max=98, radius=200, min_mag=3.0, assess_date='2026-06-12', output_dir ='/mnt/user-data/outputs/results_india/', validate=False, ) def run_italy_quick(): """Quick-start for Italy using existing catalog.""" return run_eight_steps( catalog_path='/mnt/user-data/uploads/italy_M2plus_2000_2026.csv', region_name ='Italy', lat_min=36, lat_max=48, lon_min=6, lon_max=20, radius=150, min_mag=2.0, assess_date='2026-06-12', output_dir ='/mnt/user-data/outputs/results_italy/', validate=False, ) if __name__ == '__main__': if len(sys.argv) > 1: sys.exit(main()) else: # Demo mode — run India print('No arguments provided — running India demo...') print('Usage: python six_criterion_framework.py --help') print() run_india_quick()