农业数据挖掘中的特征工程：从原始气象数据到可建模特征

import pandas as pd

df = pd.read_csv("data/weather_raw.csv", parse_dates=["date"])
missing = df.isna().mean().sort_values(ascending=False)
print(missing[missing > 0])

# 生长季（4-9 月）的缺失率远高于非生长季
df["is_growing"] = df["date"].dt.month.between(4, 9)
print(df.groupby("is_growing")["precipitation"].apply(lambda s: s.isna().mean()))

def impute_by_season(df: pd.DataFrame, cols: list[str]) -> pd.DataFrame:
    """
    对气象列按「生长季 / 非生长季」分组，用组内中位数填充。
    相比全局均值，这样能保留季节性差异。
    """
    df = df.copy()
    df["is_growing"] = df["date"].dt.month.between(4, 9)

    for col in cols:
        # 先按站点 + 生长季分组填充，最贴近该站点当季的典型值
        grouped = df.groupby(["station_id", "is_growing"])[col]
        df[col] = df[col].fillna(grouped.transform("median"))
        # 站点内依然全空的列，再退回全局中位数兜底
        df[col] = df[col].fillna(df[col].median())

    return df


weather_cols = ["temperature", "precipitation", "sunshine_hours", "humidity"]
df = impute_by_season(df, weather_cols)

import numpy as np

def add_agronomic_features(df: pd.DataFrame) -> pd.DataFrame:
    df = df.copy().sort_values(["station_id", "date"])

    # 有效积温：日平均温度高于基温（10℃）的部分累加
    BASE_TEMP = 10
    df["effective_temp"] = (df["temperature"] - BASE_TEMP).clip(lower=0)

    group = df.groupby(["station_id", "year"])
    df["gdd"] = group["effective_temp"].cumsum()        # 生长度日
    df["cum_precip"] = group["precipitation"].cumsum()  # 累计降水

    # 极端天气计数：高温日、暴雨日
    df["heat_days"] = group["temperature"].transform(lambda s: (s > 35).cumsum())
    df["heavy_rain_days"] = group["precipitation"].transform(
        lambda s: (s > 50).cumsum()
    )

    # 水分平衡：降水与蒸散的差值，反映水分是否充足
    df["water_balance"] = df["cum_precip"] - df["potential_et"]

    return df

STAGES = {
    "sowing": (4, 5),       # 播种期
    "tillering": (6, 6),    # 分蘖期
    "jointing": (7, 7),     # 拔节期
    "filling": (8, 9),      # 灌浆期
}

def aggregate_by_stage(df: pd.DataFrame) -> pd.DataFrame:
    df = df.copy()
    rows = []

    for (station, year), sub in df.groupby(["station_id", "year"]):
        record = {"station_id": station, "year": year}

        for stage, (start, end) in STAGES.items():
            stage_df = sub[sub["date"].dt.month.between(start, end)]
            if stage_df.empty:
                continue
            # 每个生育期取均值与极值，既反映水平也反映波动
            record[f"{stage}_gdd"] = stage_df["gdd"].max()
            record[f"{stage}_precip"] = stage_df["precipitation"].sum()
            record[f"{stage}_temp_mean"] = stage_df["temperature"].mean()
            record[f"{stage}_temp_max"] = stage_df["temperature"].max()
            record[f"{stage}_sunshine"] = stage_df["sunshine_hours"].mean()

        rows.append(record)

    return pd.DataFrame(rows)

from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_val_score

numeric_features = X.select_dtypes(include=np.number).columns.tolist()

preprocess = ColumnTransformer([
    ("num", Pipeline([
        ("impute", SimpleImputer(strategy="median")),
        ("scale", StandardScaler()),
    ]), numeric_features),
])

model = Pipeline([
    ("preprocess", preprocess),
    ("regressor", GradientBoostingRegressor(
        n_estimators=400,
        learning_rate=0.05,
        max_depth=3,
        random_state=42,
    )),
])

# 交叉验证：标准化只在训练折上拟合，测试折不参与，杜绝数据泄漏
scores = cross_val_score(model, X, y, cv=5, scoring="r2")
print(f"R² = {scores.mean():.3f} ± {scores.std():.3f}")