California Wildfire 5-Day Warning Modeling with Random Forest and XGBoost¶

# Install only packages missing from the active kernel.
from importlib.util import find_spec
import subprocess
import sys

REQUIRED_PACKAGES = {
    "numpy": "numpy>=1.26",
    "pandas": "pandas>=2.1",
    "pyarrow": "pyarrow>=14",
    "joblib": "joblib>=1.3",
    "xgboost": "xgboost>=3.0",
    "sklearn": "scikit-learn>=1.5",
}

missing = [
    package
    for import_name, package in REQUIRED_PACKAGES.items()
    if find_spec(import_name) is None
]

if missing:
    print("Installing missing packages:", ", ".join(missing))
    subprocess.check_call(
        [sys.executable, "-m", "pip", "install", "--quiet", *missing]
    )
else:
    print("All required packages are already available.")

# Download, verify, and safely extract Dataset A.
from pathlib import Path
from urllib.request import urlretrieve
import hashlib
import zipfile

DATASET_URL = (
    "https://storage.i-guide.io/datasets/"
    "764ff862-9f65-48df-b325-5dc616ab826f/"
    "iguide_california_wildfire_ml_2015_2020.zip"
)
EXPECTED_ARCHIVE_SHA256 = (
    "256d371290bcc920f88d1d878a3e68f4ab66e33ddd56f68d488646aaf2e017be"
)

WORK_DIR = Path.cwd() / "iguide_wildfire_model_workspace"
DOWNLOAD_DIR = WORK_DIR / "downloads"
EXTRACT_DIR = WORK_DIR / "dataset"
ZIP_PATH = DOWNLOAD_DIR / "iguide_california_wildfire_ml_2015_2020.zip"

DOWNLOAD_DIR.mkdir(parents=True, exist_ok=True)
EXTRACT_DIR.mkdir(parents=True, exist_ok=True)

def sha256_file(path: Path, chunk_size: int = 1024 * 1024) -> str:
    digest = hashlib.sha256()
    with path.open("rb") as stream:
        for chunk in iter(lambda: stream.read(chunk_size), b""):
            digest.update(chunk)
    return digest.hexdigest()

def safe_extract_zip(zip_path: Path, destination: Path) -> None:
    destination = destination.resolve()
    with zipfile.ZipFile(zip_path) as archive:
        for member in archive.infolist():
            target = (destination / member.filename).resolve()
            if destination != target and destination not in target.parents:
                raise ValueError(f"Unsafe archive member: {member.filename}")
        archive.extractall(destination)

parquet_matches = list(EXTRACT_DIR.rglob("wildfire_features.parquet"))

if not parquet_matches:
    if not ZIP_PATH.exists():
        print("Downloading Dataset A from I-GUIDE...")
        urlretrieve(DATASET_URL, ZIP_PATH)

    actual_hash = sha256_file(ZIP_PATH)
    if actual_hash != EXPECTED_ARCHIVE_SHA256:
        ZIP_PATH.unlink(missing_ok=True)
        raise ValueError(
            "Dataset A checksum mismatch. "
            f"Expected {EXPECTED_ARCHIVE_SHA256}, received {actual_hash}."
        )

    print("Archive checksum verified. Extracting Dataset A...")
    safe_extract_zip(ZIP_PATH, EXTRACT_DIR)
    parquet_matches = list(EXTRACT_DIR.rglob("wildfire_features.parquet"))

if len(parquet_matches) != 1:
    raise FileNotFoundError(
        "Expected exactly one wildfire_features.parquet after extraction; "
        f"found {len(parquet_matches)}."
    )

PARQUET_PATH = parquet_matches[0]
print(f"Dataset ready: {PARQUET_PATH}")

# Load and validate the published data contract.
import json
import numpy as np
import pandas as pd
import joblib
import xgboost as xgb
if not hasattr(xgb.XGBClassifier, "_estimator_type"):
    xgb.XGBClassifier._estimator_type = "classifier"
from IPython.display import display
from sklearn.ensemble import RandomForestClassifier
from sklearn.impute import SimpleImputer

df = pd.read_parquet(PARQUET_PATH)

EXPECTED_METADATA_COLUMNS = {"sample_id", "date", "year", "split", "label"}
missing_metadata = EXPECTED_METADATA_COLUMNS.difference(df.columns)
if missing_metadata:
    raise ValueError(f"Missing required columns: {sorted(missing_metadata)}")

feature_cols = [
    column for column in df.columns
    if column not in EXPECTED_METADATA_COLUMNS
]

assert df.shape == (18_000, 92), f"Unexpected table shape: {df.shape}"
assert len(feature_cols) == 87, f"Expected 87 features, found {len(feature_cols)}"
assert int((df["label"] == 0).sum()) == 15_000
assert int((df["label"] == 1).sum()) == 3_000
assert int((df["split"] == "train").sum()) == 14_118
assert int((df["split"] == "test").sum()) == 3_882
assert set(df["split"].dropna().unique()) == {"train", "test"}

df["date"] = pd.to_datetime(df["date"], errors="raise")
assert df.loc[df["split"] == "train", "year"].between(2015, 2019).all()
assert (df.loc[df["split"] == "test", "year"] == 2020).all()

validation_summary = pd.DataFrame(
    {
        "item": [
            "Rows",
            "Total columns",
            "Model features",
            "Positive samples",
            "Negative samples",
            "Training samples",
            "Testing samples",
        ],
        "value": [
            len(df),
            len(df.columns),
            len(feature_cols),
            int((df["label"] == 1).sum()),
            int((df["label"] == 0).sum()),
            int((df["split"] == "train").sum()),
            int((df["split"] == "test").sum()),
        ],
    }
)
display(validation_summary)

# Prepare temporal train/test arrays using training-derived preprocessing only.
X_all = df[feature_cols].to_numpy(dtype=np.float32)
y_all = df["label"].to_numpy(dtype=np.int8)

train_mask = df["split"].eq("train").to_numpy()
test_mask = df["split"].eq("test").to_numpy()

X_train, y_train = X_all[train_mask], y_all[train_mask]
X_test, y_test = X_all[test_mask], y_all[test_mask]

imputer = SimpleImputer(strategy="median")
X_train_imputed = imputer.fit_transform(X_train)
X_test_imputed = imputer.transform(X_test)

feature_min = X_train_imputed.min(axis=0)
feature_max = X_train_imputed.max(axis=0)
scale_denominator = np.maximum(feature_max - feature_min, 1e-8)

X_train_scaled = np.clip(
    (X_train_imputed - feature_min) / scale_denominator,
    0.0,
    1.0,
)
X_test_scaled = np.clip(
    (X_test_imputed - feature_min) / scale_denominator,
    0.0,
    1.0,
)

print("Training matrix:", X_train_scaled.shape)
print("Testing matrix:", X_test_scaled.shape)

# Train and evaluate Random Forest and XGBoost.
from sklearn.metrics import (
    accuracy_score,
    precision_score,
    recall_score,
    f1_score,
    roc_auc_score,
    confusion_matrix,
)

def evaluate_binary_classifier(name, y_true, probabilities, threshold=0.5):
    predictions = (probabilities >= threshold).astype(np.int8)
    tn, fp, fn, tp = confusion_matrix(
        y_true, predictions, labels=[0, 1]
    ).ravel()
    return {
        "model": name,
        "accuracy": accuracy_score(y_true, predictions),
        "precision": precision_score(y_true, predictions, zero_division=0),
        "recall": recall_score(y_true, predictions, zero_division=0),
        "f1": f1_score(y_true, predictions, zero_division=0),
        "roc_auc": roc_auc_score(y_true, probabilities),
        "true_negative": int(tn),
        "false_positive": int(fp),
        "false_negative": int(fn),
        "true_positive": int(tp),
        "threshold": threshold,
    }

print("Training Random Forest...")
rf_model = RandomForestClassifier(
    n_estimators=223,
    max_depth=11,
    min_samples_leaf=1,
    min_samples_split=3,
    random_state=42,
    n_jobs=-1,
)
rf_model.fit(X_train_scaled, y_train)
rf_probabilities = rf_model.predict_proba(X_test_scaled)[:, 1]

print("Training XGBoost...")
xgb_model = xgb.XGBClassifier(
    n_estimators=100,
    max_depth=6,
    learning_rate=0.1,
    random_state=42,
    eval_metric="logloss",
    n_jobs=-1,
)
xgb_model.fit(X_train_scaled, y_train)
xgb_probabilities = xgb_model.predict_proba(X_test_scaled)[:, 1]

metrics = [
    evaluate_binary_classifier("Random Forest", y_test, rf_probabilities),
    evaluate_binary_classifier("XGBoost", y_test, xgb_probabilities),
]
results_df = pd.DataFrame(metrics)
display(
    results_df[
        ["model", "accuracy", "precision", "recall", "f1", "roc_auc"]
    ].style.format(precision=3)
)

# Export a complete, reusable model bundle into the notebook workspace.
OUTPUT_DIR = WORK_DIR / "outputs" / "model_bundle"
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

rf_path = OUTPUT_DIR / "random_forest_model.joblib"
xgb_path = OUTPUT_DIR / "xgboost_model.json"
preprocessor_path = OUTPUT_DIR / "preprocessor_state.json"
metadata_path = OUTPUT_DIR / "model_metadata.json"

joblib.dump(rf_model, rf_path)
xgb_model.save_model(xgb_path)

preprocessor_state = {
    "feature_names": feature_cols,
    "feature_count": len(feature_cols),
    "imputer_strategy": "median",
    "imputer_medians": imputer.statistics_.astype(float).tolist(),
    "scaler_min": feature_min.astype(float).tolist(),
    "scaler_max": feature_max.astype(float).tolist(),
    "classification_threshold": 0.5,
}
preprocessor_path.write_text(
    json.dumps(preprocessor_state, indent=2),
    encoding="utf-8",
)

model_metadata = {
    "dataset_url": DATASET_URL,
    "dataset_archive_sha256": EXPECTED_ARCHIVE_SHA256,
    "training_period": "2015-2019",
    "testing_period": "2020",
    "sample_count": int(len(df)),
    "feature_count": int(len(feature_cols)),
    "train_count": int(train_mask.sum()),
    "test_count": int(test_mask.sum()),
    "metrics": metrics,
    "random_forest_parameters": rf_model.get_params(),
    "xgboost_parameters": xgb_model.get_params(),
}
metadata_path.write_text(
    json.dumps(model_metadata, indent=2, default=str),
    encoding="utf-8",
)

exported = pd.DataFrame(
    {
        "artifact": [
            "Random Forest",
            "XGBoost",
            "Preprocessor state",
            "Model metadata",
        ],
        "path": [
            str(rf_path),
            str(xgb_path),
            str(preprocessor_path),
            str(metadata_path),
        ],
    }
)
display(exported)