From bd3d73d140f9a8889d1ce8f2051a533a44e12306 Mon Sep 17 00:00:00 2001
From: zhenghu <1831829219@qq.com>
Date: Tue, 14 Apr 2026 15:49:37 +0800
Subject: [PATCH] =?UTF-8?q?feat:=20=E5=B0=86=E6=A8=A1=E6=8B=9F=E5=B9=B3?=
=?UTF-8?q?=E5=8F=B0=E4=BB=8E=E8=A5=BF=E7=8F=AD=E7=89=99=20Demo=20?=
=?UTF-8?q?=E6=95=B0=E6=8D=AE=E8=BF=81=E7=A7=BB=E8=87=B3=E4=B8=AD=E5=9B=BD?=
=?UTF-8?q?=E5=A4=9A=E7=9C=81=E4=BB=BD=E5=9C=BA=E6=99=AF?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
提交正文：
- simulator.py: 新增 SyntheticWeatherDataProvider，基于河南、黑龙江、湖北、
新疆、四川五省气候模板生成合成日气象数据
- simulator.py: 引入国内农事日历 CROP_CALENDAR_CN 与省份作物配置 PROVINCE_CROPS
- simulator.py: 移除对 GridWeatherDataProvider / AgroManagementDataProvider /
fetch_sitedata 的西班牙 Demo 数据依赖
- app.py: 侧边栏支持省份选择，年份范围扩展为 2019-2023
- app.py: 全面移除自定义 CSS，改用 Streamlit 原生组件（st.metric / st.info /
st.success / st.divider 等）简化界面
- app.py: 图表回归 Plotly 原生 add_vline，移除 hex_to_rgba / add_milestone_line
辅助函数
---
app.py | 508 ++++++++++++++++-----------------------------------
simulator.py | 392 ++++++++++++++++++++++++++++++---------
2 files changed, 460 insertions(+), 440 deletions(-)
diff --git a/app.py b/app.py
index 4633487..bdf28cf 100644
--- a/app.py
+++ b/app.py
@@ -4,39 +4,15 @@
import streamlit as st
import plotly.graph_objects as go
-from plotly.subplots import make_subplots
-from simulator import run_wofost, run_multi_crop, list_available_crops, CROP_META
+from simulator import (
+ run_wofost,
+ run_multi_crop,
+ list_available_provinces,
+ list_available_crops,
+ CROP_META,
+)
-
-def hex_to_rgba(hex_color: str, alpha: float = 0.1) -> str:
- hex_color = hex_color.lstrip("#")
- r = int(hex_color[0:2], 16)
- g = int(hex_color[2:4], 16)
- b = int(hex_color[4:6], 16)
- return f"rgba({r}, {g}, {b}, {alpha})"
-
-
-def add_milestone_line(fig, x, color: str, text: str, xref: str = "x"):
- """用 shape + annotation 手动添加竖线，绕过 Plotly add_vline 对 datetime.date 的 bug。"""
- fig.add_shape(
- type="line",
- x0=x, x1=x, y0=0, y1=1,
- xref=xref, yref="paper",
- line=dict(color=color, width=1.5, dash="dot"),
- )
- fig.add_annotation(
- x=x, y=1.02,
- xref=xref, yref="paper",
- text=text,
- showarrow=False,
- font=dict(color=color, size=10),
- bgcolor="rgba(255,255,255,0.8)",
- borderpad=2,
- )
-
-
-# ─── Page Config ────────────────────────────────────────────────────────────
st.set_page_config(
page_title="麦麦智农",
page_icon="🌾",
@@ -44,279 +20,108 @@ st.set_page_config(
initial_sidebar_state="expanded",
)
-# ─── Custom CSS ──────────────────────────────────────────────────────────────
-st.markdown("""
-
-""", unsafe_allow_html=True)
-
-
-# ─── Sidebar Inputs ──────────────────────────────────────────────────────────
+# ─── Sidebar ────────────────────────────────────────────────────────────────
with st.sidebar:
- st.markdown('', unsafe_allow_html=True)
- st.markdown('', unsafe_allow_html=True)
- st.markdown("
", unsafe_allow_html=True)
+ st.header("🌾 麦麦智农")
+ st.caption("基于 PCSE/WOFOST 的真实作物模拟")
+ st.divider()
- st.markdown('', unsafe_allow_html=True)
- st.markdown("""
-
- 西班牙 · 埃斯特雷马杜拉
- Grid 31031 (37.64°N, 6.09°W)
-
- """, unsafe_allow_html=True)
+ province = st.selectbox("省份", list_available_provinces())
+ year = st.selectbox("年份", [2019, 2020, 2021, 2022, 2023], index=3)
- st.markdown('', unsafe_allow_html=True)
- year = st.selectbox("年份", [2000], index=0, help="Demo 数据库当前仅提供 2000 年数据")
-
- crops = list_available_crops(grid_no=31031, year=year)
+ crops = list_available_crops(province)
crop_options = {c["crop_no"]: f"{c['emoji']} {c['name']}" for c in crops}
- selected_crop_no = st.selectbox("作物", list(crop_options.keys()), format_func=lambda k: crop_options[k])
+ selected_crop_no = st.selectbox(
+ "作物", list(crop_options.keys()), format_func=lambda k: crop_options[k]
+ )
- mode = st.radio("生产模式", ["pp", "wlp"], format_func=lambda m: {"pp": "潜在生产 (PP)", "wlp": "水分限制生产 (WLP)"}[m])
+ mode = st.radio(
+ "生产模式",
+ ["pp", "wlp"],
+ format_func=lambda m: {"pp": "潜在生产 (PP)", "wlp": "水分限制生产 (WLP)"}[m],
+ )
- st.markdown('', unsafe_allow_html=True)
- area = st.number_input("种植面积 (公顷)", min_value=0.1, max_value=10000.0, value=100.0, step=10.0)
+ area = st.number_input(
+ "种植面积 (公顷)", min_value=0.1, max_value=10000.0, value=100.0, step=10.0
+ )
-# ─── Run Simulation ──────────────────────────────────────────────────────────
+# ─── Run Simulation ─────────────────────────────────────────────────────────
@st.cache_data(show_spinner=False)
-def cached_run(crop_no, year, mode):
- return run_wofost(crop_no=crop_no, year=year, mode=mode)
+def cached_run(province, crop_no, year, mode):
+ return run_wofost(province=province, crop_no=crop_no, year=year, mode=mode)
@st.cache_data(show_spinner=False)
-def cached_multi_crop(year, mode):
- return run_multi_crop(year=year, mode=mode)
+def cached_multi_crop(province, year, mode):
+ return run_multi_crop(province=province, year=year, mode=mode)
with st.spinner("正在运行 WOFOST 作物模拟，请稍候..."):
- result = cached_run(selected_crop_no, year, mode)
- multi_crop_result = cached_multi_crop(year, mode)
+ result = cached_run(province, selected_crop_no, year, mode)
+ multi_crop_result = cached_multi_crop(province, year, mode)
meta = result["meta"]
summary = result["summary"]
df = result["df"]
crop_info = CROP_META.get(selected_crop_no, {"name": "未知作物", "emoji": "🌱", "color": "#888"})
-# 单位换算
twso = summary["twso"] if summary["twso"] is not None else 0.0
tagp = summary["tagp"] if summary["tagp"] is not None else 0.0
max_lai = summary["max_lai"] if summary["max_lai"] is not None else 0.0
total_twso_tons = twso * area / 1000.0
-# ─── Main Layout ─────────────────────────────────────────────────────────────
-st.markdown(f"""
-
-
麦麦智农 · {crop_info['emoji']} {crop_info['name']} 生长模拟
-
- {meta['year']} 年 · {meta['mode_label']} · 基于 WOFOST 7.2 真实作物模型
-
-""", unsafe_allow_html=True)
+# ─── Header ─────────────────────────────────────────────────────────────────
+st.title(f"{crop_info['emoji']} {crop_info['name']} 生长模拟")
+st.caption(f"{province} · {year} 年 · {meta['mode_label']} · 基于 WOFOST 7.2 真实作物模型")
-st.markdown("
", unsafe_allow_html=True)
-
-# KPI row
+# ─── KPIs ───────────────────────────────────────────────────────────────────
k1, k2, k3, k4 = st.columns(4)
-with k1:
- st.markdown(f"""
-
-
{twso:,.0f}
-
kg / 公顷
-
🌾 经济产量 (TWSO)
-
-
{tagp:,.0f}
-
kg / 公顷
-
🌿 总生物量 (TAGP)
-
-
{max_lai:.2f}
-
最大叶面积指数
-
☘️ LAI max
-
-
{total_twso_tons:,.1f}
-
吨 / 总产量
-
📦 {area:.0f} 公顷总产
-
", unsafe_allow_html=True)
+st.divider()
-# ─── Charts Row 1 ────────────────────────────────────────────────────────────
+# ─── Charts Row 1 ───────────────────────────────────────────────────────────
chart_left, chart_right = st.columns(2)
with chart_left:
- st.markdown('', unsafe_allow_html=True)
+ st.subheader("📈 叶面积指数 (LAI) 动态")
if not df.empty and "LAI" in df.columns:
fig_lai = go.Figure()
- fig_lai.add_trace(go.Scatter(
- x=df.index, y=df["LAI"],
- mode="lines",
- line=dict(color=crop_info["color"], width=2.5),
- fill="tozeroy",
- fillcolor=hex_to_rgba(crop_info["color"], 0.1),
- name="LAI",
- ))
- # 添加生育期标记
+ fig_lai.add_trace(
+ go.Scatter(
+ x=df.index,
+ y=df["LAI"],
+ mode="lines",
+ line=dict(color=crop_info["color"], width=2.5),
+ fill="tozeroy",
+ name="LAI",
+ )
+ )
ms = summary.get("milestones", {})
if "flowering" in ms:
- add_milestone_line(fig_lai, x=ms["flowering"], color="#d4a574", text="开花")
+ fig_lai.add_vline(
+ x=ms["flowering"],
+ line=dict(color="#d4a574", width=1.5, dash="dot"),
+ annotation_text="开花",
+ annotation_position="top left",
+ )
if "maturity" in ms:
- add_milestone_line(fig_lai, x=ms["maturity"], color="#7c5e42", text="成熟")
+ fig_lai.add_vline(
+ x=ms["maturity"],
+ line=dict(color="#7c5e42", width=1.5, dash="dash"),
+ annotation_text="成熟",
+ annotation_position="top right",
+ )
fig_lai.update_layout(
- xaxis=dict(title="日期", color="#5a5a5a", gridcolor="rgba(0,0,0,0.06)"),
- yaxis=dict(title="LAI", color="#5a5a5a", gridcolor="rgba(0,0,0,0.06)"),
- paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)",
- font=dict(color="#2c2c2c", size=10),
- margin=dict(t=20, b=36, l=50, r=20), height=280,
+ xaxis_title="日期",
+ yaxis_title="LAI",
+ margin=dict(t=20, b=36, l=50, r=20),
+ height=300,
showlegend=False,
)
st.plotly_chart(fig_lai, use_container_width=True)
@@ -324,60 +129,62 @@ with chart_left:
st.info("暂无 LAI 数据")
with chart_right:
- st.markdown('', unsafe_allow_html=True)
+ st.subheader("📊 生物量与产量积累")
if not df.empty and "TAGP" in df.columns:
fig_bio = go.Figure()
- fig_bio.add_trace(go.Scatter(
- x=df.index, y=df["TAGP"],
- mode="lines",
- line=dict(color="#7a9e7e", width=2.5),
- name="TAGP",
- ))
+ fig_bio.add_trace(
+ go.Scatter(
+ x=df.index, y=df["TAGP"], mode="lines", name="TAGP", line=dict(color="#7a9e7e", width=2.5)
+ )
+ )
if "TWSO" in df.columns:
- fig_bio.add_trace(go.Scatter(
- x=df.index, y=df["TWSO"],
- mode="lines",
- line=dict(color=crop_info["color"], width=2.5),
- name="TWSO",
- ))
+ fig_bio.add_trace(
+ go.Scatter(
+ x=df.index,
+ y=df["TWSO"],
+ mode="lines",
+ name="TWSO",
+ line=dict(color=crop_info["color"], width=2.5),
+ )
+ )
ms = summary.get("milestones", {})
if "flowering" in ms:
- add_milestone_line(fig_bio, x=ms["flowering"], color="#d4a574", text="开花")
+ fig_bio.add_vline(x=ms["flowering"], line=dict(color="#d4a574", width=1.5, dash="dot"))
if "maturity" in ms:
- add_milestone_line(fig_bio, x=ms["maturity"], color="#7c5e42", text="成熟")
+ fig_bio.add_vline(x=ms["maturity"], line=dict(color="#7c5e42", width=1.5, dash="dash"))
fig_bio.update_layout(
- xaxis=dict(title="日期", color="#5a5a5a", gridcolor="rgba(0,0,0,0.06)"),
- yaxis=dict(title="干物质 (kg/ha)", color="#5a5a5a", gridcolor="rgba(0,0,0,0.06)"),
- paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)",
- font=dict(color="#2c2c2c", size=10),
- margin=dict(t=20, b=36, l=50, r=20), height=280,
- legend=dict(orientation="h", y=-0.18, font=dict(size=10)),
+ xaxis_title="日期",
+ yaxis_title="干物质 (kg/ha)",
+ margin=dict(t=20, b=36, l=50, r=20),
+ height=300,
+ legend=dict(orientation="h", y=-0.18),
)
st.plotly_chart(fig_bio, use_container_width=True)
else:
st.info("暂无生物量数据")
-# ─── Charts Row 2 ────────────────────────────────────────────────────────────
+# ─── Charts Row 2 ───────────────────────────────────────────────────────────
chart2_left, chart2_right = st.columns(2)
with chart2_left:
- st.markdown('', unsafe_allow_html=True)
+ st.subheader("💧 土壤水分动态")
if not df.empty and "SM" in df.columns:
fig_sm = go.Figure()
- fig_sm.add_trace(go.Scatter(
- x=df.index, y=df["SM"],
- mode="lines",
- line=dict(color="#5a8f9e", width=2.5),
- fill="tozeroy",
- fillcolor="rgba(90, 143, 158, 0.1)",
- name="SM",
- ))
+ fig_sm.add_trace(
+ go.Scatter(
+ x=df.index,
+ y=df["SM"],
+ mode="lines",
+ line=dict(color="#5a8f9e", width=2.5),
+ fill="tozeroy",
+ name="SM",
+ )
+ )
fig_sm.update_layout(
- xaxis=dict(title="日期", color="#5a5a5a", gridcolor="rgba(0,0,0,0.06)"),
- yaxis=dict(title="土壤含水量 (cm³/cm³)", color="#5a5a5a", gridcolor="rgba(0,0,0,0.06)"),
- paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)",
- font=dict(color="#2c2c2c", size=10),
- margin=dict(t=20, b=36, l=50, r=20), height=260,
+ xaxis_title="日期",
+ yaxis_title="土壤含水量 (cm³/cm³)",
+ margin=dict(t=20, b=36, l=50, r=20),
+ height=280,
showlegend=False,
)
st.plotly_chart(fig_sm, use_container_width=True)
@@ -385,84 +192,87 @@ with chart2_left:
st.info("暂无土壤水分数据")
with chart2_right:
- st.markdown('', unsafe_allow_html=True)
+ st.subheader("🏅 同一年份作物产量对比")
if multi_crop_result:
- names = [f"{CROP_META.get(r['crop_no'], {}).get('emoji', '🌱')} {r['name']}" for r in multi_crop_result]
+ names = [
+ f"{CROP_META.get(r['crop_no'], {}).get('emoji', '🌱')} {r['name']}"
+ for r in multi_crop_result
+ ]
twsos = [r["twso"] if r["twso"] is not None else 0.0 for r in multi_crop_result]
colors = [CROP_META.get(r["crop_no"], {}).get("color", "#888") for r in multi_crop_result]
fig_bar = go.Figure()
- fig_bar.add_trace(go.Bar(
- x=names, y=twsos,
- marker=dict(color=colors, opacity=0.85,
- line=dict(color="rgba(0,0,0,0.08)", width=1)),
- text=[f"{v:,.0f}" for v in twsos],
- textposition="outside",
- textfont=dict(color="#5a5a5a", size=10),
- ))
+ fig_bar.add_trace(
+ go.Bar(
+ x=names,
+ y=twsos,
+ marker=dict(color=colors, opacity=0.85),
+ text=[f"{v:,.0f}" for v in twsos],
+ textposition="outside",
+ )
+ )
fig_bar.update_layout(
- xaxis=dict(color="#5a5a5a", gridcolor="rgba(0,0,0,0.04)"),
- yaxis=dict(title="经济产量 (kg/ha)", color="#5a5a5a", gridcolor="rgba(0,0,0,0.06)"),
- paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)",
- font=dict(color="#2c2c2c", size=11),
- margin=dict(t=20, b=30, l=60, r=20), height=260,
+ yaxis_title="经济产量 (kg/ha)",
+ margin=dict(t=20, b=30, l=60, r=20),
+ height=280,
showlegend=False,
)
st.plotly_chart(fig_bar, use_container_width=True)
else:
st.info("暂无对比数据")
-# ─── Advisory Panel ───────────────────────────────────────────────────────────
-st.markdown('', unsafe_allow_html=True)
+# ─── Advisory Panel ─────────────────────────────────────────────────────────
+st.divider()
+st.subheader("💡 农艺建议与生育期")
+
adv1, adv2 = st.columns(2)
with adv1:
- advisories = []
if twso > 6000:
- advisories.append(("good", f"{crop_info['name']} 模拟产量达到 {twso:,.0f} kg/ha，属于高产水平，气候与土壤条件匹配良好。"))
+ st.success(
+ f"{crop_info['name']} 模拟产量达到 {twso:,.0f} kg/ha，属于高产水平，气候与土壤条件匹配良好。"
+ )
elif twso > 3000:
- advisories.append(("good", f"{crop_info['name']} 模拟产量为 {twso:,.0f} kg/ha，处于中等水平，可通过优化水肥管理进一步提升。"))
+ st.info(
+ f"{crop_info['name']} 模拟产量为 {twso:,.0f} kg/ha，处于中等水平，可通过优化水肥管理进一步提升。"
+ )
else:
- advisories.append(("warn", f"{crop_info['name']} 模拟产量仅 {twso:,.0f} kg/ha，建议检查品种适宜性或水分胁迫情况。"))
+ st.warning(
+ f"{crop_info['name']} 模拟产量仅 {twso:,.0f} kg/ha，建议检查品种适宜性或水分胁迫情况。"
+ )
if mode == "wlp":
if not df.empty and "SM" in df.columns:
sm_min = df["SM"].min()
if sm_min < 0.15:
- advisories.append(("warn", f"模拟期间土壤水分最低降至 {sm_min:.3f}，出现明显水分胁迫，建议评估灌溉方案。"))
+ st.warning(
+ f"模拟期间土壤水分最低降至 {sm_min:.3f}，出现明显水分胁迫，建议评估灌溉方案。"
+ )
else:
- advisories.append(("good", "模拟期间土壤水分状况总体良好，未出现极端干旱胁迫。"))
+ st.success("模拟期间土壤水分状况总体良好，未出现极端干旱胁迫。")
else:
- advisories.append(("good", "当前为潜在生产模式，结果反映理想水肥条件下的产量上限。"))
+ st.info("当前为潜在生产模式，结果反映理想水肥条件下的产量上限。")
if max_lai > 5:
- advisories.append(("good", f"最大 LAI 达到 {max_lai:.2f}，冠层覆盖充分，光能截获效率高。"))
+ st.success(f"最大 LAI 达到 {max_lai:.2f}，冠层覆盖充分，光能截获效率高。")
elif max_lai < 2:
- advisories.append(("warn", f"最大 LAI 仅 {max_lai:.2f}，冠层发育不足，可能存在播期或品种问题。"))
-
- for typ, msg in advisories:
- css_class = "alert-good" if typ == "good" else "alert-warn"
- st.markdown(f'{msg}
', unsafe_allow_html=True)
+ st.warning(f"最大 LAI 仅 {max_lai:.2f}，冠层发育不足，可能存在播期或品种问题。")
with adv2:
ms = summary.get("milestones", {})
if ms:
- st.markdown(f"""
-
- 生育期里程碑
- 播种/出苗：{ms.get('start', '—')}
- 开花期：{ms.get('flowering', '—')}
- 成熟期：{ms.get('maturity', '—')}
- 收获/结束：{ms.get('end', '—')}
- 生育期天数：{summary.get('duration', '—')} 天
-
- """, unsafe_allow_html=True)
- else:
- st.markdown('暂无生育期数据
', unsafe_allow_html=True)
+ st.info(
+ f"""
+**生育期里程碑**
-# ─── Footer ───────────────────────────────────────────────────────────────────
-st.markdown("
", unsafe_allow_html=True)
-st.markdown("""
-
- 麦麦智农 · 基于 PCSE/WOFOST 真实作物生长模型 · 结果仅供参考
-
-""", unsafe_allow_html=True)
+- 播种/出苗：{ms.get('start', '—')}
+- 开花期：{ms.get('flowering', '—')}
+- 成熟期：{ms.get('maturity', '—')}
+- 收获/结束：{ms.get('end', '—')}
+- 生育期天数：{summary.get('duration', '—')} 天
+ """
+ )
+ else:
+ st.info("暂无生育期数据")
+
+st.divider()
+st.caption("麦麦智农 · 基于 PCSE/WOFOST 真实作物生长模型 · 结果仅供参考")
diff --git a/simulator.py b/simulator.py
index 0e2a1c9..1904d00 100644
--- a/simulator.py
+++ b/simulator.py
@@ -1,37 +1,27 @@
"""
-麦麦智农 - PCSE 模拟引擎封装
+麦麦智农 - PCSE 模拟引擎封装（国内场景版）
基于 WOFOST 7.2 潜在生产（PP）与水分限制生产（WLP）模型
"""
+import datetime
+import math
import sqlite3
from collections import namedtuple
from typing import Literal
+import numpy as np
import pandas as pd
-import pcse
from pcse.models import Wofost72_PP, Wofost72_WLP_CWB
from pcse.settings import settings
-from pcse.base import ParameterProvider
-from pcse.tests.db_input import (
- GridWeatherDataProvider,
- fetch_soildata,
- fetch_sitedata,
- fetch_cropdata,
- AgroManagementDataProvider,
-)
+from pcse.base import ParameterProvider, WeatherDataContainer, WeatherDataProvider
+from pcse.tests.db_input import fetch_soildata, fetch_cropdata
+from pcse.util import reference_ET
-def namedtuple_factory(cursor, row):
- fields = [column[0] for column in cursor.description]
- cls = namedtuple("Row", fields)
- return cls._make(row)
-
-
-# Demo 数据库配置
+# ─── 常量配置 ────────────────────────────────────────────────────────────────
DB_PATH = settings.PCSE_USER_HOME + "/pcse.db"
-# 作物映射（crop_no -> 中文名/元数据）
CROP_META = {
1: {"name": "冬小麦", "emoji": "🌾", "color": "#c69c5d", "en": "WINTER WHEAT"},
2: {"name": "玉米", "emoji": "🌽", "color": "#e8a93f", "en": "GRAIN MAIZE"},
@@ -46,6 +36,129 @@ MODE_LABELS = {
"wlp": "水分限制生产",
}
+# 省份支持作物（黑龙江太冷，去掉越冬作物）
+PROVINCE_CROPS = {
+ "河南": [1, 2, 3, 7, 10, 11],
+ "黑龙江": [2, 3, 7, 11],
+ "湖北": [1, 2, 3, 7, 10, 11],
+ "新疆": [1, 2, 3, 7, 10, 11],
+ "四川": [1, 2, 3, 7, 10, 11],
+}
+
+# 国内农事日历（播种月日 -> 收获月日）
+CROP_CALENDAR_CN = {
+ "冬小麦": {
+ "河南": {"start": (10, 15), "end": (5, 25)},
+ "湖北": {"start": (10, 25), "end": (5, 15)},
+ "新疆": {"start": (9, 20), "end": (6, 25)},
+ "四川": {"start": (10, 20), "end": (5, 10)},
+ },
+ "玉米": {
+ "河南": {"start": (6, 5), "end": (9, 25)},
+ "黑龙江": {"start": (5, 1), "end": (9, 20)},
+ "湖北": {"start": (3, 25), "end": (7, 20)},
+ "新疆": {"start": (4, 15), "end": (9, 10)},
+ "四川": {"start": (3, 20), "end": (7, 15)},
+ },
+ "春大麦": {
+ "河南": {"start": (3, 1), "end": (5, 25)},
+ "黑龙江": {"start": (4, 10), "end": (7, 15)},
+ "湖北": {"start": (2, 15), "end": (5, 20)},
+ "新疆": {"start": (3, 25), "end": (7, 10)},
+ "四川": {"start": (2, 20), "end": (5, 15)},
+ },
+ "马铃薯": {
+ "河南": {"start": (2, 25), "end": (5, 30)},
+ "黑龙江": {"start": (4, 20), "end": (8, 15)},
+ "湖北": {"start": (1, 15), "end": (5, 10)},
+ "新疆": {"start": (4, 1), "end": (7, 25)},
+ "四川": {"start": (1, 20), "end": (5, 5)},
+ },
+ "冬油菜": {
+ "河南": {"start": (9, 15), "end": (5, 15)},
+ "湖北": {"start": (9, 20), "end": (5, 5)},
+ "新疆": {"start": (8, 25), "end": (6, 15)},
+ "四川": {"start": (9, 10), "end": (4, 25)},
+ },
+ "向日葵": {
+ "河南": {"start": (5, 25), "end": (8, 30)},
+ "黑龙江": {"start": (5, 15), "end": (9, 10)},
+ "湖北": {"start": (4, 5), "end": (8, 5)},
+ "新疆": {"start": (4, 20), "end": (8, 20)},
+ "四川": {"start": (3, 25), "end": (7, 25)},
+ },
+}
+
+# 省份气候模板（每月 [TMAX, TMIN, RAIN(mm/day), IRRAD(kJ/m2/day)]）
+PROVINCES = {
+ "河南": {
+ "lat": 34.76, "lon": 113.65, "elev": 110,
+ "rh": 0.65, "wind_base": 2.0,
+ "climate": [
+ [5, -4, 0.5, 9000], [9, -1, 0.6, 12000], [15, 4, 1.0, 16000],
+ [22, 10, 1.5, 19000], [28, 16, 2.0, 22000], [33, 21, 2.5, 24000],
+ [32, 23, 4.0, 21000], [31, 22, 3.5, 19000], [27, 17, 2.0, 16000],
+ [22, 11, 1.2, 13000], [14, 3, 0.8, 9000], [7, -2, 0.4, 8000],
+ ],
+ },
+ "黑龙江": {
+ "lat": 45.80, "lon": 126.53, "elev": 150,
+ "rh": 0.60, "wind_base": 3.0,
+ "climate": [
+ [-13, -24, 0.2, 6000], [-8, -20, 0.2, 9000], [2, -9, 0.4, 14000],
+ [13, 1, 1.0, 18000], [21, 8, 1.8, 22000], [26, 15, 3.0, 24000],
+ [28, 18, 4.5, 22000], [26, 16, 4.0, 19000], [20, 9, 1.5, 15000],
+ [11, 1, 0.8, 11000], [-1, -10, 0.4, 7000], [-11, -21, 0.3, 5000],
+ ],
+ },
+ "湖北": {
+ "lat": 30.59, "lon": 114.31, "elev": 30,
+ "rh": 0.75, "wind_base": 1.8,
+ "climate": [
+ [8, 0, 1.5, 8000], [11, 3, 1.8, 10000], [16, 7, 2.5, 13000],
+ [22, 13, 3.0, 16000], [27, 18, 3.5, 19000], [30, 22, 4.5, 20000],
+ [33, 25, 5.5, 22000], [33, 25, 4.0, 21000], [29, 21, 3.0, 17000],
+ [24, 15, 2.0, 14000], [17, 8, 1.5, 10000], [11, 2, 1.0, 8000],
+ ],
+ },
+ "新疆": {
+ "lat": 44.30, "lon": 87.60, "elev": 450,
+ "rh": 0.45, "wind_base": 2.5,
+ "climate": [
+ [-8, -16, 0.1, 7000], [-4, -12, 0.2, 11000], [6, -2, 0.3, 16000],
+ [18, 7, 0.5, 20000], [25, 13, 0.8, 24000], [30, 18, 1.0, 26000],
+ [32, 20, 1.2, 27000], [30, 18, 0.8, 25000], [24, 12, 0.5, 20000],
+ [15, 4, 0.3, 15000], [3, -4, 0.2, 9000], [-5, -13, 0.1, 6000],
+ ],
+ },
+ "四川": {
+ "lat": 30.67, "lon": 104.07, "elev": 500,
+ "rh": 0.80, "wind_base": 1.2,
+ "climate": [
+ [10, 3, 0.6, 6000], [13, 5, 0.8, 8000], [17, 9, 1.2, 11000],
+ [23, 14, 2.0, 13000], [27, 18, 3.0, 15000], [29, 21, 4.0, 16000],
+ [31, 23, 5.0, 17000], [31, 23, 4.5, 16000], [27, 19, 3.0, 12000],
+ [22, 15, 1.5, 9000], [16, 10, 1.0, 7000], [11, 5, 0.5, 6000],
+ ],
+ },
+}
+
+SITEDATA = {
+ "IFUNRN": 0.0,
+ "SSMAX": 0.0,
+ "NOTINF": 0.0,
+ "SSI": 0.0,
+ "WAV": 100.0,
+ "SMLIM": 0.0,
+}
+
+
+# ─── 数据库辅助 ───────────────────────────────────────────────────────────────
+def namedtuple_factory(cursor, row):
+ fields = [column[0] for column in cursor.description]
+ cls = namedtuple("Row", fields)
+ return cls._make(row)
+
def get_db_conn():
conn = sqlite3.connect(DB_PATH)
@@ -53,68 +166,172 @@ def get_db_conn():
return conn
-def list_available_years(grid_no: int = 31031) -> list[int]:
- """返回 Demo 数据库中指定 grid 有气象数据的年份列表。"""
- conn = get_db_conn()
- c = conn.cursor()
- rows = c.execute(
- "SELECT DISTINCT strftime('%Y', day) as year FROM grid_weather WHERE grid_no=? ORDER BY year",
- (grid_no,),
- ).fetchall()
- conn.close()
- return [int(r.year) for r in rows]
+# ─── 气象数据生成器 ───────────────────────────────────────────────────────────
+class SyntheticWeatherDataProvider(WeatherDataProvider):
+ """基于省份气候模板和年份生成合成日气象数据。"""
+
+ angstA = 0.29
+ angstB = 0.49
+
+ def __init__(self, province: str, start_date: datetime.date, end_date: datetime.date, seed: int):
+ WeatherDataProvider.__init__(self)
+ self.province = PROVINCES[province]
+ self.latitude = self.province["lat"]
+ self.longitude = self.province["lon"]
+ self.elevation = self.province["elev"]
+ self.description = f"Synthetic weather for {province}"
+
+ rng = np.random.default_rng(seed)
+ days = (end_date - start_date).days + 1
+ days_in_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
+ if (start_date.year % 4 == 0 and start_date.year % 100 != 0) or (start_date.year % 400 == 0):
+ days_in_month[1] = 29
+
+ for i in range(days):
+ day = start_date + datetime.timedelta(days=i)
+ month = day.month - 1
+ next_month = (month + 1) % 12
+ d = day.day - 1
+ frac = (d + 0.5) / days_in_month[month]
+ w2 = 0.5 - 0.5 * math.cos(math.pi * frac)
+ w1 = 1 - w2
+
+ c = self.province["climate"][month]
+ c_next = self.province["climate"][next_month]
+
+ tmax_base = c[0] * w1 + c_next[0] * w2
+ tmin_base = c[1] * w1 + c_next[1] * w2
+ rain_mean = c[2] * w1 + c_next[2] * w2
+ irrad_base = (c[3] * w1 + c_next[3] * w2) * 1000.0 # kJ -> J
+
+ tmax = tmax_base + rng.normal(0, 1.5)
+ tmin = tmin_base + rng.normal(0, 1.5)
+ tmin = min(tmin, tmax - 1.0)
+ irrad = max(0, irrad_base * (1 + rng.normal(0, 0.1)))
+
+ p_rain = min(0.8, max(0.05, rain_mean / 5.0))
+ if rng.random() < p_rain:
+ scale = rain_mean / (2 * p_rain)
+ rain = rng.gamma(2, scale) if scale > 0 else 0.0
+ else:
+ rain = 0.0
+ rain = max(0.0, rain)
+
+ tmean = (tmax + tmin) / 2.0
+ es = 0.6108 * math.exp(17.27 * tmean / (tmean + 237.3))
+ rh = self.province["rh"] + rng.normal(0, 0.05)
+ rh = max(0.2, min(0.95, rh))
+ vap = es * rh * 10.0
+
+ wind = max(0.5, self.province["wind_base"] + rng.normal(0, 0.8))
+
+ e0, es0, et0 = reference_ET(
+ DAY=day, LAT=self.latitude, ELEV=self.elevation,
+ TMIN=float(tmin), TMAX=float(tmax), IRRAD=float(irrad),
+ VAP=float(vap), WIND=float(wind),
+ ANGSTA=self.angstA, ANGSTB=self.angstB,
+ )
+
+ t = {
+ "DAY": day,
+ "LAT": self.latitude,
+ "LON": self.longitude,
+ "ELEV": self.elevation,
+ "TMAX": float(tmax),
+ "TMIN": float(tmin),
+ "VAP": float(vap),
+ "WIND": float(wind),
+ "RAIN": float(rain) / 10.0,
+ "IRRAD": float(irrad),
+ "E0": float(e0) / 10.0,
+ "ES0": float(es0) / 10.0,
+ "ET0": float(et0) / 10.0,
+ "SNOWDEPTH": 0.0,
+ }
+ wdc = WeatherDataContainer(**t)
+ self._store_WeatherDataContainer(wdc, day)
-def list_available_crops(grid_no: int = 31031, year: int = 2000) -> list[dict]:
- """返回指定 grid 与年份有作物日历的作物列表。"""
- conn = get_db_conn()
- c = conn.cursor()
- rows = c.execute(
- "SELECT crop_no FROM crop_calendar WHERE grid_no=? AND year=? ORDER BY crop_no",
- (grid_no, year),
- ).fetchall()
- conn.close()
- result = []
- for r in rows:
- meta = CROP_META.get(r.crop_no)
- if meta:
- result.append({"crop_no": r.crop_no, **meta})
- return result
+# ─── 农事管理构造 ───────────────────────────────────────────────────────────
+def build_agromanagement(crop_no: int, province: str, year: int):
+ meta = CROP_META[crop_no]
+ crop_name = meta["en"]
+ cal = CROP_CALENDAR_CN[meta["name"]][province]
+ sm, sd = cal["start"]
+ em, ed = cal["end"]
+
+ start_date = datetime.date(year, sm, sd)
+ end_date = datetime.date(year, em, ed)
+ if end_date < start_date:
+ end_date = datetime.date(year + 1, em, ed)
+ campaign_start = datetime.date(year, 1, 1)
+ else:
+ campaign_start = datetime.date(year, 1, 1)
+
+ return [{
+ campaign_start: {
+ "CropCalendar": {
+ "crop_name": crop_name,
+ "variety_name": f"{crop_name}_{province}_{year}",
+ "crop_start_date": start_date,
+ "crop_start_type": "emergence",
+ "crop_end_date": end_date,
+ "crop_end_type": "earliest",
+ "max_duration": 365,
+ },
+ "TimedEvents": None,
+ "StateEvents": None,
+ }
+ }]
+
+
+# ─── 公开接口 ───────────────────────────────────────────────────────────────
+def list_available_provinces() -> list[str]:
+ return list(PROVINCES.keys())
+
+
+def list_available_crops(province: str) -> list[dict]:
+ return [{"crop_no": c, **CROP_META[c]} for c in PROVINCE_CROPS.get(province, [])]
def run_wofost(
- grid_no: int = 31031,
- crop_no: int = 1,
- year: int = 2000,
+ province: str,
+ crop_no: int,
+ year: int,
mode: Literal["pp", "wlp"] = "wlp",
) -> dict:
- """
- 运行 WOFOST 模拟，返回包含输出时间序列与关键指标的字典。
- """
+ agromanagement = build_agromanagement(crop_no, province, year)
+
+ # 确定气象起止日期（覆盖整个生育期 + 前后缓冲）
+ campaign = agromanagement[0]
+ campaign_date = list(campaign.keys())[0]
+ cc = campaign[campaign_date]["CropCalendar"]
+ start = campaign_date - datetime.timedelta(days=30)
+ end = cc["crop_end_date"] + datetime.timedelta(days=30)
+
+ seed = year * 1000 + hash(province) % 100000
+ wdp = SyntheticWeatherDataProvider(province, start, end, seed)
+
conn = get_db_conn()
try:
- agromanagement = AgroManagementDataProvider(conn, grid_no, crop_no, year)
- sited = fetch_sitedata(conn, grid_no, year)
- cropd = fetch_cropdata(conn, grid_no, year, crop_no)
- soild = fetch_soildata(conn, grid_no)
- parvalues = ParameterProvider(sitedata=sited, soildata=soild, cropdata=cropd)
- wdp = GridWeatherDataProvider(conn, grid_no=grid_no)
-
- if mode == "pp":
- wofsim = Wofost72_PP(parvalues, wdp, agromanagement)
- else:
- wofsim = Wofost72_WLP_CWB(parvalues, wdp, agromanagement)
-
- wofsim.run_till_terminate()
- output = wofsim.get_output()
+ # 作物/土壤参数复用 PCSE Demo DB（生理参数通用）
+ cropd = fetch_cropdata(conn, 31031, 2000, crop_no)
+ soild = fetch_soildata(conn, 31031)
+ parvalues = ParameterProvider(sitedata=SITEDATA, soildata=soild, cropdata=cropd)
finally:
conn.close()
+ if mode == "pp":
+ wofsim = Wofost72_PP(parvalues, wdp, agromanagement)
+ else:
+ wofsim = Wofost72_WLP_CWB(parvalues, wdp, agromanagement)
+
+ wofsim.run_till_terminate()
+ output = wofsim.get_output()
df = pd.DataFrame(output)
if not df.empty and "day" in df.columns:
df = df.set_index("day")
- # 关键指标
if not df.empty:
last = df.iloc[-1]
max_lai = df["LAI"].max() if "LAI" in df.columns else None
@@ -126,28 +343,27 @@ def run_wofost(
else:
max_lai = tagp = twso = sm_end = dvs_end = duration = None
- # 生育期里程碑（基于 DVS）
milestones = {}
+ crop_start = cc["crop_start_date"]
+ crop_end = cc["crop_end_date"]
if not df.empty and "DVS" in df.columns:
dvs_series = df["DVS"]
- # 出苗/开始生长 DVS ~0（实际在 calendar 开始日）
- # 开花 DVS ~1.0
- # 成熟 DVS ~2.0
flowering = dvs_series[dvs_series >= 1.0]
if not flowering.empty:
milestones["flowering"] = flowering.index[0]
maturity = dvs_series[dvs_series >= 2.0]
if not maturity.empty:
milestones["maturity"] = maturity.index[0]
- # 开始日
- milestones["start"] = df.index[0]
- milestones["end"] = df.index[-1]
+ milestones["start"] = crop_start
+ milestones["end"] = df.index[-1] if df.index[-1] <= crop_end else crop_end
+
+ actual_duration = (milestones.get("end", crop_end) - crop_start).days + 1 if milestones else None
meta = CROP_META.get(crop_no, {"name": "未知作物", "emoji": "🌱", "color": "#888"})
return {
"meta": {
- "grid_no": grid_no,
+ "province": province,
"crop_no": crop_no,
"year": year,
"mode": mode,
@@ -156,30 +372,27 @@ def run_wofost(
},
"df": df,
"summary": {
- "tagp": tagp, # 总地上生物量 kg/ha
- "twso": twso, # 经济产量 kg/ha
+ "tagp": tagp,
+ "twso": twso,
"max_lai": max_lai,
"sm_end": sm_end,
"dvs_end": dvs_end,
- "duration": duration,
+ "duration": actual_duration,
"milestones": milestones,
},
}
def run_multi_crop(
- grid_no: int = 31031,
- year: int = 2000,
+ province: str,
+ year: int,
mode: Literal["pp", "wlp"] = "wlp",
) -> list[dict]:
- """
- 对指定年份所有可用作物运行模拟，返回各作物汇总结果列表（按 TWSO 排序）。
- """
- crops = list_available_crops(grid_no, year)
+ crops = list_available_crops(province)
results = []
for c in crops:
try:
- res = run_wofost(grid_no=grid_no, crop_no=c["crop_no"], year=year, mode=mode)
+ res = run_wofost(province=province, crop_no=c["crop_no"], year=year, mode=mode)
results.append({
"crop_no": c["crop_no"],
"name": c["name"],
@@ -191,26 +404,23 @@ def run_multi_crop(
"duration": res["summary"]["duration"],
})
except Exception:
- # 某些作物可能因数据缺失导致模拟失败，静默跳过
continue
results.sort(key=lambda x: (x["twso"] if x["twso"] is not None else -1), reverse=True)
return results
def run_multi_year(
- grid_no: int = 31031,
- crop_no: int = 1,
+ province: str,
+ crop_no: int,
mode: Literal["pp", "wlp"] = "wlp",
+ years: list[int] | None = None,
) -> list[dict]:
- """
- 对指定作物在所有可用年份运行模拟，返回各年份汇总结果列表。
- """
- years = list_available_years(grid_no)
+ if years is None:
+ years = list(range(2019, 2024))
results = []
- meta = CROP_META.get(crop_no, {"name": "未知作物", "emoji": "🌱", "color": "#888"})
for year in years:
try:
- res = run_wofost(grid_no=grid_no, crop_no=crop_no, year=year, mode=mode)
+ res = run_wofost(province=province, crop_no=crop_no, year=year, mode=mode)
results.append({
"year": year,
"tagp": res["summary"]["tagp"],