安徽农户持续种粮调研：公开数据分析与结果复核

from pathlib import Path
import hashlib
import json
import platform
import re

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from IPython.display import display

candidates = [Path.cwd(), Path.cwd().parent, Path.cwd() / "public-release"]
ROOT = next(
    (path.resolve() for path in candidates if (path / "data").is_dir()),
    None,
)
if ROOT is None:
    raise FileNotFoundError("未找到公开仓库根目录及 data 目录。")

DATA_ROOT = ROOT / "data"
RESULT_DIR = DATA_ROOT / "results"
FIGURE_DATA_DIR = DATA_ROOT / "charts"

import sys
from matplotlib_inline.backend_inline import set_matplotlib_formats
sys.path.insert(0, str(ROOT / "scripts"))
from plot_style import configure, finish_figure
from table_style import show_table
font_audit = configure()
set_matplotlib_formats("svg", "png")
print("绘图字体：" + " / ".join(font_audit["families"]))
np.random.seed(20260914)

print("公开仓库根目录：已定位")
print(f"Python：{platform.python_version()}")
print(f"pandas：{pd.__version__}；NumPy：{np.__version__}")

data_files = sorted(
    path for path in DATA_ROOT.rglob("*")
    if path.is_file() and path.suffix.lower() in {".csv", ".json"}
)

inventory = []
for path in data_files:
    relative = path.relative_to(ROOT).as_posix()
    if path.suffix.lower() == ".csv":
        frame = pd.read_csv(path)
        inventory.append(
            {"文件": relative, "格式": "CSV", "行数": len(frame), "字段数": len(frame.columns)}
        )
    else:
        value = json.loads(path.read_text(encoding="utf-8-sig"))
        inventory.append(
            {
                "文件": relative,
                "格式": "JSON",
                "行数": len(value) if hasattr(value, "__len__") else 1,
                "字段数": len(value) if isinstance(value, dict) else np.nan,
            }
        )

inventory_df = pd.DataFrame(inventory)
show_table(inventory_df, "表1　公开数据文件清单")
print(f"公开研究数据文件：{len(data_files)} 个")

risky_field_parts = {
    "姓名", "手机号", "联系电话", "身份证", "邮箱", "电子邮件",
    "详细地址", "家庭住址", "答卷id", "受访者编号",
    "respondent_id", "phone", "mobile", "email", "id_card",
    "identity_number", "ip_address", "openid", "unionid",
}
value_patterns = {
    "中国大陆手机号": re.compile(r"(?<![\d./])1[3-9]\d{9}(?![\d./])"),
    "电子邮箱": re.compile(r"\b[A-Za-z0-9._%+-]+@[A-Za-z0-9.-]+\.[A-Za-z]{2,}\b"),
    "IPv4": re.compile(
        r"(?<![\d.])(?:25[0-5]|2[0-4]\d|1?\d?\d)"
        r"(?:\.(?:25[0-5]|2[0-4]\d|1?\d?\d)){3}(?![\d.])"
    ),
    "微信ID": re.compile(r"(?i)\bwxid_[a-z0-9]+\b"),
    "Windows个人路径": re.compile(r"(?i)\b[A-Z]:\\"),
    "Linux个人路径": re.compile(r"/home/[A-Za-z0-9._-]+/"),
}

risky_fields = []
direct_value_hits = []
for path in data_files:
    relative = path.relative_to(ROOT).as_posix()
    text = path.read_text(encoding="utf-8-sig")
    if path.suffix.lower() == ".csv":
        fields = list(pd.read_csv(path, nrows=0).columns)
    else:
        value = json.loads(text)
        fields = []
        stack = [value]
        while stack:
            item = stack.pop()
            if isinstance(item, dict):
                fields.extend(map(str, item.keys()))
                stack.extend(item.values())
            elif isinstance(item, list):
                stack.extend(item)
    for field in fields:
        normalized = re.sub(r"[\s-]+", "_", field.casefold())
        if any(part in normalized for part in risky_field_parts):
            risky_fields.append((relative, field))
    for label, pattern in value_patterns.items():
        if pattern.search(text):
            direct_value_hits.append((relative, label))

assert not risky_fields, f"发现高风险字段：{risky_fields}"
assert not direct_value_hits, f"发现直接标识模式：{direct_value_hits}"
print("直接标识检查：通过")
print("高风险字段：0；手机号/邮箱/IP/微信ID/个人路径模式：0")

land_status = pd.read_csv(RESULT_DIR / "land-status-2025.csv")
labor_groups = pd.read_csv(RESULT_DIR / "fallow-groups-2025.csv")
show_table(land_status, "表2　2025年土地经营状态", {"percent": "{:.2f}%"})
show_table(labor_groups, "表3　不同务农劳动力组的撂荒比例", {"fallow_percent": "{:.2f}%"})

fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))
sns.barplot(data=land_status, x="percent", y="state", color="#555555", ax=axes[0])
axes[0].set(title="2025年土地经营状态", xlabel="比例（%）", ylabel="")
sns.barplot(data=labor_groups, x="group", y="fallow_percent", color="#555555", ax=axes[1])
axes[1].set(title="不同家庭务农劳动力组的撂荒比例", xlabel="家庭务农劳动力", ylabel="撂荒比例（%）")
finish_figure(fig)

fallow_model = pd.read_csv(RESULT_DIR / "fallow-model-2025.csv")
model_view = fallow_model.loc[fallow_model["term"] != "截距"].copy()
model_view["置信区间是否跨1"] = ~(
    (model_view["ci_low_or"] > 1) | (model_view["ci_high_or"] < 1)
)
show_table(
    model_view[
        ["term", "odds_ratio", "ci_low_or", "ci_high_or", "bootstrap_sign_p", "n", "置信区间是否跨1"]
    ], "表4　2025年撂荒模型优势比与区间",
        {
            "odds_ratio": "{:.3f}", "ci_low_or": "{:.3f}",
            "ci_high_or": "{:.3f}", "bootstrap_sign_p": "{:.4f}",
        }
)

ordered = model_view.sort_values("odds_ratio")
y = np.arange(len(ordered))
fig, ax = plt.subplots(figsize=(9, 5))
ax.errorbar(
    ordered["odds_ratio"], y,
    xerr=[ordered["odds_ratio"] - ordered["ci_low_or"], ordered["ci_high_or"] - ordered["odds_ratio"]],
    fmt="o", markerfacecolor="white", color="#555555", ecolor="#555555", capsize=4,
)
ax.axvline(1, color="#999999", linestyle="--", linewidth=1)
ax.set(yticks=y, yticklabels=ordered["term"], xlabel="优势比 OR（95%置信区间）", title="2025年撂荒模型森林图")
finish_figure(fig)

prediction = pd.read_csv(RESULT_DIR / "predictions-2025.csv")
fig, ax = plt.subplots(figsize=(9, 5))
sns.lineplot(
    data=prediction,
    x="operated_area_mu", y="predicted_fallow_percent",
    hue="disaster", style="disaster", linewidth=1.3, ax=ax,
)
ax.set(xlabel="经营面积（亩）", ylabel="预测撂荒比例（%）", title="经营面积与受灾情景下的预测撂荒比例")
ax.legend(title="受灾情景")
finish_figure(fig)

sensitivity = pd.read_csv(RESULT_DIR / "sensitivity-2025.csv")
show_table(sensitivity, "表5　2025年敏感性检验", {"fallow_percent": "{:.2f}%", "labor_or": "{:.3f}", "area_or": "{:.3f}", "disaster_or": "{:.3f}"})

sensitivity_long = sensitivity.melt(
    id_vars=["scenario", "n"],
    value_vars=["labor_or", "area_or", "disaster_or"],
    var_name="指标", value_name="OR",
)
label_map = {"labor_or": "务农劳动力", "area_or": "经营面积", "disaster_or": "受灾"}
sensitivity_long["指标"] = sensitivity_long["指标"].map(label_map)
fig, ax = plt.subplots(figsize=(11, 5))
sns.pointplot(data=sensitivity_long, x="OR", y="scenario", hue="指标", dodge=0.35,
              markers=["o", "s", "^"], linestyles=["-", "--", ":"], ax=ax)
ax.axvline(1, color="#999999", linestyle="--", linewidth=1)
ax.set(xlabel="优势比 OR", ylabel="", title="不同敏感性口径的关键优势比")
finish_figure(fig)

profiles = pd.read_csv(RESULT_DIR / "profiles-2025.csv")
profile_meta = json.loads((RESULT_DIR / "cluster-metadata-2025.json").read_text(encoding="utf-8"))
profile_columns = [
    "profile", "n", "share_percent", "labor", "operated_area", "grain_area",
    "fallow_percent", "managed_percent", "transfer_out_percent", "transfer_in_percent",
]
show_table(profiles[profile_columns], "表6　2025年经营画像规模与特征", {column: "{:.2f}%" if column == "fallow_percent" else "{:.2f}" for column in profile_columns if column not in {"profile", "n"}})
show_table(pd.Series(profile_meta, name="值").rename_axis("参数").reset_index(), "表7　画像模型参数与稳定性")

fig, axes = plt.subplots(1, 2, figsize=(13, 4.8))
sns.barplot(data=profiles, x="share_percent", y="profile", color="#555555", ax=axes[0])
axes[0].set(title="四类经营画像的样本占比", xlabel="占比（%）", ylabel="")
profile_rates = profiles.melt(
    id_vars="profile", value_vars=["fallow_percent", "managed_percent", "transfer_out_percent", "transfer_in_percent"],
    var_name="指标", value_name="比例",
)
rate_names = {
    "fallow_percent": "撂荒", "managed_percent": "托管",
    "transfer_out_percent": "转出", "transfer_in_percent": "转入",
}
profile_rates["指标"] = profile_rates["指标"].map(rate_names)
sns.barplot(data=profile_rates, x="比例", y="profile", hue="指标", ax=axes[1])
axes[1].set(title="不同画像的土地经营特征", xlabel="比例（%）", ylabel="")
finish_figure(fig)

household = pd.read_csv(RESULT_DIR / "households-2026.csv")
household["主口径复算"] = (household["success"] / household["n"] * 100).round(2)
household["敏感性复算"] = (household["sensitivity_success"] / household["sensitivity_n"] * 100).round(2)
show_table(household[["label", "n", "success", "percent", "主口径复算", "denominator_rule"]],
           "表8　2026年农户端指标与主口径复算", {"percent": "{:.2f}%", "主口径复算": "{:.2f}%"})
show_table(household[["label", "sensitivity_n", "sensitivity_success", "sensitivity_percent", "敏感性复算"]],
           "表9　2026年农户端排除短时答卷后的复算", {"sensitivity_percent": "{:.2f}%", "敏感性复算": "{:.2f}%"})

fig, ax = plt.subplots(figsize=(10, 7))
plot_household = household.sort_values("percent")
sns.barplot(data=plot_household, x="percent", y="label", color="#555555", ax=ax)
ax.set(title="2026年农户端公开指标", xlabel="比例（%）", ylabel="")
finish_figure(fig)

village = pd.read_csv(RESULT_DIR / "villages-2026.csv")
village["主口径复算"] = (village["success"] / village["n"] * 100).round(2)
village["敏感性复算"] = (village["sensitivity_success"] / village["sensitivity_n"] * 100).round(2)
show_table(village[["label", "n", "success", "percent", "主口径复算", "denominator_rule"]],
           "表10　2026年村级端指标与主口径复算", {"percent": "{:.2f}%", "主口径复算": "{:.2f}%"})
show_table(village[["label", "sensitivity_n", "sensitivity_success", "sensitivity_percent", "敏感性复算"]],
           "表11　2026年村级端排除短时答卷后的复算", {"sensitivity_percent": "{:.2f}%", "敏感性复算": "{:.2f}%"})

fig, ax = plt.subplots(figsize=(10, 6))
plot_village = village.sort_values("percent")
sns.barplot(data=plot_village, x="percent", y="label", color="#555555", ax=ax)
ax.set(title="2026年村级公开指标", xlabel="比例（%）", ylabel="")
finish_figure(fig)

parallel = pd.read_csv(RESULT_DIR / "household-village-2026.csv")
show_table(parallel, "表12　2026年农户与村级并列指标", {"household_percent": "{:.2f}%", "village_percent": "{:.2f}%"})

parallel_long = pd.concat(
    [
        parallel[["theme", "household_percent"]].rename(columns={"household_percent": "比例"}).assign(层级="农户端"),
        parallel[["theme", "village_percent"]].rename(columns={"village_percent": "比例"}).assign(层级="村级端"),
    ],
    ignore_index=True,
)
fig, ax = plt.subplots(figsize=(9, 4.8))
sns.barplot(data=parallel_long, x="theme", y="比例", hue="层级", ax=ax)
ax.set(title="2026年农户端与村级端并列指标", xlabel="主题", ylabel="比例（%）")
finish_figure(fig)

figure_prediction = pd.read_csv(FIGURE_DATA_DIR / "fig4-4_area_disaster_prediction.csv").rename(
    columns={"category": "disaster", "value": "predicted_fallow_percent"}
)
prediction_sorted = prediction.sort_values(["disaster", "operated_area_mu"]).reset_index(drop=True)
figure_prediction_sorted = figure_prediction.sort_values(["disaster", "operated_area_mu"]).reset_index(drop=True)

model_figure = pd.read_csv(FIGURE_DATA_DIR / "fig4-3_logit_or.csv")
model_expected = model_view.rename(columns={"term": "category"})[
    ["category", "odds_ratio", "ci_low_or", "ci_high_or", "n"]
].sort_values("category").reset_index(drop=True)
model_figure_sorted = model_figure.sort_values("category").reset_index(drop=True)

profile_figure = pd.read_csv(FIGURE_DATA_DIR / "fig5-3_profile_metrics.csv")
profile_expected = profiles[profile_figure.columns].sort_values("cluster").reset_index(drop=True)
profile_figure_sorted = profile_figure.sort_values("cluster").reset_index(drop=True)

checks = {
    "公开数据文件数为24": len(data_files) == 24,
    "2025土地状态最大有效分母为2079": int(land_status["n"].max()) == 2079,
    "2025撂荒模型有效样本为2055": fallow_model["n"].nunique() == 1 and int(fallow_model["n"].iloc[0]) == 2055,
    "敏感性检验状态全部PASS": sensitivity["status"].eq("PASS").all(),
    "四类画像占比合计为100%": np.isclose(profiles["share_percent"].sum(), 100.0),
    "预测情景制图数据一致": np.allclose(
        prediction_sorted[["operated_area_mu", "predicted_fallow_percent"]],
        figure_prediction_sorted[["operated_area_mu", "predicted_fallow_percent"]],
    ) and prediction_sorted["disaster"].equals(figure_prediction_sorted["disaster"]),
    "模型森林图数据一致": model_expected["category"].equals(model_figure_sorted["category"])
    and np.allclose(model_expected.drop(columns="category"), model_figure_sorted.drop(columns="category")),
    "画像制图数据一致": profile_expected["profile"].equals(profile_figure_sorted["profile"])
    and np.allclose(
        profile_expected.drop(columns="profile"),
        profile_figure_sorted.drop(columns="profile"),
    ),
    "2026农户端比例复算一致": np.allclose(household["percent"], household["主口径复算"], atol=0.011)
    and np.allclose(household["sensitivity_percent"], household["敏感性复算"], atol=0.011),
    "2026村级端比例复算一致": np.allclose(village["percent"], village["主口径复算"], atol=0.011)
    and np.allclose(village["sensitivity_percent"], village["敏感性复算"], atol=0.011),
    "直接标识检查通过": not risky_fields and not direct_value_hits,
}
checks_df = pd.Series(checks, name="通过").rename_axis("检查项").to_frame()
show_table(checks_df.reset_index(), "表13　数据一致性检查结果")
assert all(checks.values()), checks
print(f"一致性检查：{sum(checks.values())}/{len(checks)} 项通过")