2. 数据获取与来源整合

作物产量预测模型的精度，很大程度上取决于输入数据的完整性与一致性。一个可用的产量预测数据集通常需要同时覆盖四类信息：气象条件（温度、降水、辐射）、土壤属性（质地、有机质、pH、持水量）、遥感观测（植被指数、地表温度、物候）以及田间管理记录（播期、品种、施肥、灌溉、收获）。这四类数据在时间分辨率、空间尺度、格式和获取方式上差异极大，因此“获取”与“整合”必须作为一个整体来设计，而不是先随便下载再临时拼接。

一、四类数据源的采集渠道与特征

气象数据最常用的公开渠道是中国气象数据网提供的站点观测数据，以及 NASA POWER、ERA5-Land 等再分析产品。站点数据精度高但空间稀疏，再分析数据覆盖连续但存在偏差，实践中常以站点数据校准再分析数据。土壤数据可来自第二次全国土壤普查的县级剖面资料、HWSD 世界土壤数据库，或通过土壤类型图与实测采样点结合生成栅格属性。遥感数据以 Sentinel-2、Landsat 8/9 和 MODIS 为主，分别对应 10 m、30 m 和 250 m～1 km 的空间分辨率，可通过 Copernicus Data Space、USGS EarthExplorer 或 Google Earth Engine 获取。田间管理数据一般没有公开 API，需要从农场记录本、合作社台账或农机作业日志中整理，通常以 Excel 或 CSV 形式存在。

这四类数据在接入方式上可归纳为三种：API 在线调用（气象再分析、部分遥感目录服务）、文件批量读取（站点 CSV、土壤属性表、管理记录表）和栅格影像处理（遥感 GeoTIFF）。下面的图示展示了从四类来源到统一数据表的整合路径，以及每类数据在接入时需要处理的关键字段。

二、API 调用：以气象再分析数据为例

NASA POWER API 是教学与原型开发中较易上手的气象数据接口，支持按经纬度和日期范围返回 JSON 或 CSV。以下示例请求某地块（经度 116.4，纬度 39.9）2023 年 4 月至 9 月的日最高温、最低温和降水，参数标识为 T2M_MAX、T2M_MIN、PRECTOTCORR：

import requests
import pandas as pd

params = {
    "parameters": "T2M_MAX,T2M_MIN,PRECTOTCORR",
    "community": "AG",
    "longitude": 116.4,
    "latitude": 39.9,
    "start": "20230401",
    "end": "20230930",
    "format": "JSON"
}
url = "https://power.larc.nasa.gov/api/temporal/daily/point"
resp = requests.get(url, params=params, timeout=30)
resp.raise_for_status()
data = resp.json()["properties"]["parameter"]

df_weather = pd.DataFrame(data)
df_weather.index = pd.to_datetime(df_weather.index, format="%Y%m%d")
df_weather = df_weather.replace(-999.0, pd.NA)   # 缺失值哨兵
print(df_weather.head())

调用时需要注意三点：一是接口返回的缺失值常用 -999 表示，必须显式替换，否则会污染后续统计；二是再分析数据存在数天延迟，批量拉取时应按年份分段并加入重试；三是同一地块的经纬度应固定记录在配置表中，避免不同批次数据落在不同网格。对于需要多地块的场景，可把请求封装为函数，按地块循环调用并追加 site_id 列。

三、文件读取：站点气象与管理台账

站点气象数据通常以 CSV 提供，常见问题是编码（GBK 与 UTF-8 混用）、日期格式不统一、列名含中文与单位。读取时应显式指定编码与日期解析列：

df_station = pd.read_csv(
    "weather_station_2023.csv",
    encoding="gbk",
    parse_dates=["观测日期"],
    na_values=["", "999999", "—"]
)
df_station = df_station.rename(columns={
    "观测日期": "date", "平均气温": "tmean",
    "降水量": "precip", "站点号": "station_id"
})

田间管理台账多为人工填写，字段命名随意，例如“播种时间”“播期”“播种日期”可能同时存在。建议先做字段映射表，再统一为 sowing_date、variety、n_rate、irrigation_mm 等标准列。对于 Excel 多工作表文件，可用 pd.read_excel(path, sheet_name=None) 一次性读入全部工作表，再按表名筛选。管理数据的关键是地块编号必须与气象、遥感数据使用同一套编码，否则后续无法按主键合并。

四、遥感栅格：从影像到地块指数

遥感数据以 GeoTIFF 为主，读取可使用 rasterio。典型流程是：读取影像 → 应用云掩膜 → 按地块边界做分区统计 → 得到每个地块每个日期的 NDVI 均值。下面示例读取单景 NDVI 并按矢量地块提取均值：

import rasterio
import geopandas as gpd
from rasterstats import zonal_stats

with rasterio.open("NDVI_20230715.tif") as src:
    ndvi = src.read(1)
    transform = src.transform
    nodata = src.nodata

plots = gpd.read_file("plots.geojson")
stats = zonal_stats(
    plots, ndvi, affine=transform,
    nodata=nodata, stats=["mean", "count"]
)
plots["ndvi_mean"] = [s["mean"] for s in stats]

实际处理中还要解决两个问题：一是时间合成，单景影像常受云污染，需按旬或半月做最大值合成；二是坐标系统一，遥感影像多为 UTM 投影，而地块边界可能是 WGS84 经纬度，分区统计前必须重投影到同一坐标系。若地块很小而影像分辨率较粗（如 MODIS 250 m），一个像元可能覆盖多个地块，此时应记录像元覆盖比例，避免把混合像元值当作地块真实值。

五、整合阶段的关键对齐规则

四类数据接入后，整合的核心是“时空对齐”。时间上，气象是日尺度、遥感是旬尺度、管理是事件尺度，建议统一聚合到生育期阶段（播种—出苗、出苗—拔节、拔节—抽雄等），每个阶段计算气象累计值或均值、遥感指数均值，管理措施按发生日期归入对应阶段。空间上，以地块编号为主键，气象按最近站点或网格插值、土壤按地块质心提取、遥感按分区统计。整合完成后应输出一份数据字典，记录每个字段的来源、单位、时间范围和缺失率，这份字典在后续建模和排障时比代码本身更有价值。