1. 先搞清楚这份数据集是什么、能干什么、不能干什么
聊“国家级地面气象站基本气象要素日值数据集”之前，我习惯先说清楚它的物理形态：它是一份以“逐站、逐日”为粒度的长表数据，全国两千四百多个国家级地面观测站，从建站或1951年前后起，每天一条记录，覆盖气温、气压、相对湿度、降水、蒸发、风向风速、日照时数、地温这几类最基础的气象要素。数据以纯文本分发，按要素拆成多个文件，每个文件里一年一个分片，字段用逗号分隔，数值全部做过毫米、摄氏度、百帕这类物理量的“放大十倍”的整数化处理，同时附一列质量控制码。它的价值不在于“新”，而在于“长”和“全”——几十年的连续观测、统一口径、统一质控规范，这是任何商业气象接口都给不了的东西。这篇内容我打算把它拆到能直接上手用的程度：结构怎么读、坑在哪、代码怎么写、不同目标下参数怎么选，适合刚接触气象数据的数据分析人员、做农业和能源评估的工程人员，以及需要给模型找驱动数据的研究者参考。
1.1 三句话定位：它是一份“逐站逐日”的长表
第一句，它是地面观测，不是再分析、不是卫星反演、不是格点产品。这意味着每一条记录背后对应一个具体的观测场，有经纬度、有海拔、有建站和迁站历史，空间上是离散的点，不做插值时你拿到的是“点真值”，做插值时你必须自己承担插值误差。
第二句，它是日值，不是小时值，也不是月值。日值意味着观测时次被聚合过了，比如日平均气温是几个定时观测的平均，日降水量是20时至次日20时的累计（部分要素还有08—20时、20—08时的分段量），日最高最低气温带出现时间。聚合方式决定了你不能拿日值去反推小时尺度的极端强度，比如“一小时下了多少毫米”这种问题，日值回答不了。
第三句，它是经过质控但保留质控标记的。这一点极其关键：数据发布方没有把可疑值直接删掉，而是原样保留、另起一列打标记。这个设计对老手是福音，对新手是陷阱——你要是忽略那一列质控码，直接对全表求平均，那平均值的可信度就只能靠运气了。
注意：不同批次、不同年份的分发版本，在字段顺序、分隔符、缺测码取值上都存在细微差异，凡是本文给出的字段名和码值，都以“常见形态”描述，正式使用前必须拿随数据附带的说明文档核对一遍。
1.2 谁在真正用它：四类典型用户
第一类是气候与生态环境方向的分析者。他们关心的是趋势：某区域近六十年平均气温的线性倾向率是多少、无霜期是不是在延长、日降水≥50毫米的暴雨日数有没有增加、连续无降水日数（气象干旱的一个粗指标）如何变化。这类工作对数据的连续性要求极高，站点一旦迁站或中断观测，序列就断了，必须做均一性检验。
第二类是农业与生态模型的使用者。作物生长模型、蒸散模型、病虫害发生期预测模型，几乎都需要日尺度的最高气温、最低气温、降水、相对湿度、日照时数、风速作为驱动。他们要的是“没有缺测的连续序列”，所以插补工作量往往比分析工作量还大。
第三类是能源与工程设计方向的人。光伏选址要看日照时数和晴空比例，风电初筛要看平均风速和极大风速，建筑暖通设计要取近三十年最冷月平均温度、极端最低温度这类统计量。这类应用对“极值”的敏感度远高于对“均值”的敏感度，一个错误的极端值就能让设备选型偏好几个档。
第四类是做数据建模和教学的。这份数据字段干净、量级适中、真实世界的脏数据特征齐全（缺测、异常、格式不一致），是非常好的数据清洗与特征工程练手材料。
1.3 三条能力边界，越早知道越好
边界一：空间代表性有限。国家级站大多设在县城附近的开阔地，对城市建成区、山区、荒漠腹地的代表性并不好。你要研究城市热岛，用城郊一个站的日值去代表整片城区，结论会偏得离谱。
边界二：日值不能替代过程。日平均气温相同的两天，可能一天是昼夜温差两度，另一天是温差二十度。涉及作物冷害、材料老化、人体舒适度这类对温度过程敏感的问题，日值只能做粗筛。
边界三：质控码不是万能的。标记为“正确”的值不代表它在物理上一定合理，标记为“可疑”的值也不代表一定错。极端天气事件本身就容易触发质控规则的报警。我见过把一次真实寒潮的日最低气温当成错误值删掉的案例，删完之后序列平滑了，可结论也废了。
2. 拆开文件看结构：字段、单位与质量控制码
2.1 文件命名与目录组织规律
这类数据集的文件名通常遵循“数据集前缀 + 要素代码 + 版本或区划号 + 年份”的组合方式，比如某要素的年度文件形如SURF_CLI_CHN_MUL_DAY-TEM-001-2010.TXT。前缀表明这是地面日值产品，中间的三字母代码是要素标识（TEM 气温、PRS 气压、RHU 相对湿度、PRE 降水、EVP 蒸发、WIN 风、SSD 日照、GST 地温），后面的数字是版本或区划编号，最后是年份。
理解这个命名规则的实际意义在于：你可以用文件名直接做批处理，而不必把几十个 G 的文本全部读进内存。一个典型的目录会是这样：
raw/ TEM/ SURF_CLI_CHN_MUL_DAY-TEM-001-1951.TXT SURF_CLI_CHN_MUL_DAY-TEM-001-1952.TXT ... PRE/ PRS/我一般会在raw同级建一个parquet目录，把清洗后的结果按“要素 + 年份”或直接按“要素全量”存成 parquet。文本读取慢、占空间、每次都重新解析，转换一次之后后续所有分析都快得多，这一步的投入产出比非常高。
2.2 八个核心要素的字段清单与单位换算表
字段的通用头部是一致的：区站号、纬度、经度、观测场海拔、年、月、日。后面才是各要素的具体观测项，末尾跟对应的质控码。区站号是五位数的主键，跨表关联全靠它。
单位是这个数据集最典型的操作门槛。原始文件里的数值几乎全部是整数化的物理量，你要自己除回去。常见换算关系我整理成了下面这张表，这张表建议直接存进你的项目笔记里：
| 要素 | 典型字段 | 原始存储单位 | 换算方式 | 合理值域参考（换算前） |
|---|---|---|---|---|
| 气温 | 平均/最高/最低气温 | 0.1 摄氏度 | 除以 10 | -600 至 500 |
| 气压 | 平均/最高/最低气压 | 0.1 百帕 | 除以 10 | 5000 至 11000 |
| 相对湿度 | 平均/最小相对湿度 | 1 百分比 | 保持原值 | 0 至 100 |
| 降水 | 20—20时降水量等 | 0.1 毫米 | 除以 10 | 0 至 3000 |
| 蒸发 | 小型/大型蒸发量 | 0.1 毫米 | 除以 10 | 0 至 1000 |
| 风 | 平均风速/最大风速/极大风速 | 0.1 米每秒 | 除以 10 | 0 至 1000 |
| 日照 | 日照时数 | 0.1 小时 | 除以 10 | 0 至 240 |
| 地温 | 0厘米地温等 | 0.1 摄氏度 | 除以 10 | -600 至 800 |
这张表最大的用处不是换算，而是做值域扫描。你读完一年数据之后，对每个数值列求一次最大值和最小值，凡是超出上表区间范围的，全部打出来看。我做过统计，一个中等规模的分析项目里，八成的数据问题都能被这一步提前抓出来——剩下的两成才是真正需要人工判断的疑难杂症。
还有两个容易被忽略的字段。一是经纬度，部分分发版本里它是以“度分”形式的整数存储的，整数部分代表度，剩余部分代表分，需要除以相应系数才能换算成十进制度。别死记文档，读完之后把换算结果和你已知的某个站点位置对一次，误差在两三公里以内就说明换算方式对了。二是观测场海拔，这个字段在做气压订正和温度垂直订正时是必需项，很多人把它当成无关列直接扔了，等到要做海拔相关分析时又得重新读一遍原始文件。
2.3 质量控制码：最容易被跳过、也最容易毁掉结论的一列
质控码的常见取值形态是：0 表示正确，1 表示可疑，2 表示错误，8 表示缺测，9 表示未做质量控制。不同版本、不同要素的码表可能有细微差别，所以第一件事永远是拿说明文档核对，第二件事是在代码里把码值分布打印出来看一眼。
这里有一个非常典型的判断题：可疑值到底留不留？
我的做法是分场景走两条路。做气候平均态、趋势分析时，我倾向于剔除质控码为 1 和 2 的记录，因为这些工作对单个值的误差不敏感，但对系统性偏差敏感，留着可疑值会让趋势线抖得不干净。做极端事件统计分析时，我做反向操作：只剔除码值 2 和 8，保留 1，同时在论文或报告的脚注里写清楚“可疑值予以保留，以避免系统性低估极端事件强度”。原因很实在——极端值最容易触发质控规则，一刀切剔掉，你的极端事件频率会被显著压低。
提示：不论选哪种策略，都建议在清洗脚本里保留一个
qc_flag原始列，而不是直接覆盖。数据分析最怕的就是“不可回溯”，半年后你自己都不记得当时为什么删掉那批值。
还有一类陷阱是质控码与缺测码的混淆。缺测在观测值列里通常用一个特殊整数表示（日值产品里常见的是 32766，也有版本用 999xxx 系列），而质控码列用 8 表示。这两件事要分开处理：先把特殊整数替换成 NaN，再根据质控码决定是否进一步剔除。顺序搞反了，你的均值会被 32766 拉到大气层外。
3. 从原始TXT到可分析宽表的完整实操
3.1 环境准备与目录规划
工具链不用复杂，pandas + numpy + geopandas（可选）+ matplotlib 就够了。版本上建议 pandas 1.5 以上，主要是为了read_csv的dtype和na_values参数行为更稳定。
目录我会这样规划，目的是让“原始数据永不被修改”成为硬约束：
project/ raw/ # 原始TXT，只读 interim/ # 解压、转码后的中间文件 parquet/ # 清洗后的列存 meta/ # 台站元数据、说明文档 src/ # 脚本 out/ # 图表、统计结果这个习惯看着啰嗦，实际上是救命的。我遇到过最惨的一次事故，是写了个循环直接在原始目录里做to_csv覆盖，结果某年的数据被截断了一半，只能重新申请下载。从那以后raw目录的权限我直接设成只读。
3.2 先试读三行，别急着写循环
新手最常见的动作是：一看文件有几百个，立刻写一个 for 循环批量读。然后报错，然后加 try/except 把错误吞掉，最后得到一份莫名其妙少了三分之一站点的数据集。
正确的第一步是打开一个文件，读三行，肉眼看。
path = "raw/TEM/SURF_CLI_CHN_MUL_DAY-TEM-001-2010.TXT" # 编码是最常见的第一个坑，先试 gbk，不行再试 utf-8 with open(path, "r", encoding="gbk", errors="replace") as f: for _ in range(3): print(repr(f.readline()))打印出来的东西要确认四件事：分隔符是逗号、制表符还是空格；有没有表头行；字段数量是多少；行尾有没有多余的空白或回车符。这四件事确认完，你后面写的读取参数才是有依据的，而不是从某篇博客里抄来的。
我踩过一次很隐蔽的坑：同一批数据里，早年份的文件用空格分隔，晚年份的文件用逗号分隔。如果不做逐文件探测，直接用统一分隔符读，会得到一堆只有一两列的数据框，而且不报错——因为 pandas 会把整行当成一个字段。所以我现在会写一个极简的探测函数：
def sniff_sep(path, enc="gbk"): with open(path, "r", encoding=enc, errors="replace") as f: line = f.readline() for sep in [",", "\t", " "]: if line.count(sep) >= 5: return sep return ","粗糙，但足够用。
3.3 批量读取与合并的可用代码
确认格式之后，就可以写正式的读取逻辑了。下面这段是我在多个项目里反复用过的骨架，改改字段名就能跑：
import glob import os import pandas as pd import numpy as np TEM_COLS = [ "STA", "LAT", "LON", "ALT", "YEAR", "MON", "DAY", "TEM_AVG", "TEM_MAX", "TEM_MIN", "QC_AVG", "QC_MAX", "QC_MIN", ] def read_temp_year(path): df = pd.read_csv( path, header=None, names=TEM_COLS, sep=sniff_sep(path), engine="python", encoding="gbk", na_values=[32766, 999998, 999999], dtype={"STA": str}, # 区站号必须当字符串，前导零会丢 ) df["SRC_FILE"] = os.path.basename(path) return df files = sorted(glob.glob("raw/TEM/*.TXT")) frames = [] for p in files: try: frames.append(read_temp_year(p)) except Exception as e: print("FAILED:", p, e) tem = pd.concat(frames, ignore_index=True)这段代码里有三个细节值得单独说。
dtype={"STA": str}：区站号是五位数字，但只要有一批数据里出现了以 0 开头的编号，用默认的整型读入就会丢掉前导零，之后和其他表关联时全部对不上。这个错误不会报错，只会让你在 merge 之后发现行数为 0，然后怀疑人生。
na_values=[32766, ...]：把特殊缺测码在读取阶段就转成 NaN，比读完之后再替换安全得多，因为你不会忘记。
SRC_FILE这一列：我几乎在所有批处理场景里都会加这么一列来源标识。出问题的时候，你能立刻定位到是哪一年的哪个文件，而不是对着一百万行数据发呆。
合并完之后马上做一次体检，别急着往下走：
print(tem.shape) print(tem["YEAR"].value_counts().sort_index().head(20)) print(tem[["TEM_AVG", "TEM_MAX", "TEM_MIN"]].describe()) print(tem[["QC_AVG", "QC_MAX", "QC_MIN"]].apply(pd.Series.value_counts))第一行看总量，第二行看年份分布有没有断档，第三行看值域是否离谱（比如最大值 32000 就说明缺测码没替换干净），第四行看质控码分布是否异常（如果某个年份 90% 都是 8，那年的数据基本可以判定为整体缺失）。
3.4 三步清洗：缺测、质控码、单位
清洗顺序很重要，我固定按“先物理、后逻辑、再业务”的三步走。
第一步，物理层清洗：处理缺测和不可能值。缺测在读取阶段已经处理，这一步主要抓“物理上不可能”的值，比如相对湿度大于 100、日照时数大于 24 小时、风速为负。这类值数量通常很少，但必须处理，因为它们几乎必然是格式错位或字段错位造成的。
第二步，逻辑层清洗：处理质控码和内部一致性。这里要建立一套规则，比如“最高气温必须大于等于最低气温”、“20—20时降水量应大于等于08—20时与20—08时之和的较小者（允许小量观测误差）”。违反一致性规则的行，我一般标记为新的一列flag_consistency，而不是直接删，因为这类行往往对应真实的大风、强对流天气过程，删了会损失信息。
第三步，业务层清洗：按你的分析目标做取舍。做趋势分析的，剔除质控码 1、2；做极端统计的，只剔 2 和 8。这一步的规则要写死在配置里，别散落在代码各处。
单位换算我会用一个字典统一管理，避免“这列除了 10、那列忘了除”的低级错误：
SCALE = { "TEM_AVG": 0.1, "TEM_MAX": 0.1, "TEM_MIN": 0.1, "PRE_20_20": 0.1, "WIN_AVG": 0.1, "SSD": 0.1, } for col, k in SCALE.items(): if col in tem.columns: tem[col] = tem[col] * k乘 0.1 比除 10 在浮点运算上略稳妥一点点，虽然差异微乎其微，但这是个好习惯——尤其是当你要处理几十个要素、上亿行数据的时候，任何一点系统性的数值处理差异都会被放大。
3.5 站点元数据对齐与空间核查
日值文件和台站元数据是两张表，靠区站号关联。元数据里通常包含站名、经纬度、观测场海拔、建站时间、撤站时间、台站级别、迁站记录。这张表的价值远不止“加个站名好看”，它决定了你的分析是否成立。
第一个核查点是站点数量随时间的变化。把每年有有效观测的站点数画出来，你会看到一条不平滑的曲线——早年站点少，中间逐步增加，某些年份有一批站点集体消失（通常是迁站或观测任务调整）。如果你的分析区间跨越了这些断点，就必须说明“站点集合在时间上不一致”，否则算出来的区域平均值会混入“站点变化”这个伪信号。
第二个核查点是迁站记录。一个站从海拔 50 米搬到海拔 300 米的山坡上，气温序列会在迁站那年出现一个台阶式的跳变，这个跳变是真实的，但它不代表气候变化。做长序列趋势分析时的标准做法是做均一性检验（比如 SNHT 方法），把这类非气候突变识别出来并做订正，或者干脆把序列截断在迁站之后。
第三个核查点是经纬度与海拔的合理性。部分老站的经纬度记录精度很低，甚至有过录入错误。我做空间图的时候会先把所有站点打在地图上，凡是落在国境线外、海洋里、或者和相邻站点位置完全重合的，全部挑出来人工核对。这一步看着笨，但极其有效——我抓到过三个坐标录反了经纬度的站点，其中一个还被我拿去算过区域平均，纠正前后那片区域的结果差了 0.4 摄氏度。
4. 不同目标下的用法差异与参数选择
4.1 气候趋势与极端事件统计
趋势分析最忌讳的是“直接对原始序列做线性回归然后报一个斜率”。正确做法至少包含三步：均一性处理、趋势估计、显著性检验。
趋势估计我一般同时给出两种结果：普通最小二乘的线性倾向率，和 Sen 斜率。前者大家熟悉，后者对异常值更稳健，两者差得多就说明序列里有离群点需要单独看。显著性检验用 Mann-Kendall，不做检验直接报斜率在那个年代已经过不去了。
import pymannkendall as mk s = tem.groupby("YEAR")["TEM_AVG"].mean().dropna() res = mk.original_test(s.values, alpha=0.05) print(res.trend, res.p, res.slope)极端事件指标我常用这一组，每一项都有明确的农业或水文含义：
| 指标 | 定义 | 典型用途 |
|---|---|---|
| 霜冻日数 | 日最低气温 < 0 摄氏度的日数 | 农业生长期、无霜期长度 |
| 暴雨日数 | 日降水量 ≥ 50 毫米的日数 | 城市内涝风险、水土流失 |
| 连续无降水日数 | 连续日降水量 < 0.1 毫米的最长天数 | 气象干旱粗指标 |
| 暖夜日数 | 日最低气温 > 25 摄氏度的日数 | 人体健康、作物呼吸消耗 |
| 极端高温日数 | 日最高气温 ≥ 35 摄氏度的日数 | 电力负荷、人群热暴露 |
计算这些指标时有一个必须提醒的细节：阈值判定要用换算后的物理量。我见过有人拿原始整数直接和 35 比，结果 350（即 35.0 摄氏度）被判成极端高温，而 34.9 摄氏度也被判成 350，全表都是极端高温日。这个错误在代码里看起来毫不违和，因为df["TEM_MAX"] >= 35这句本身没有任何语法问题。
4.2 农业气象与积温计算
积温是最常用的派生量，公式简单：GDD = Σ max(Tavg - Tbase, 0)，Tbase 按作物取，喜温作物取 10 摄氏度，耐寒作物取 5 或 0 摄氏度。看起来一行代码就能算完，但实际落地的坑全在细节上。
第一个细节：用(Tmax + Tmin) / 2算日平均，还是直接用观测的日平均气温？两者在多数站点差别不大，但在昼夜温差大的干旱区能差出 0.5 摄氏度以上。农业模型里我倾向用观测的日平均气温，因为它更接近气象站的正式统计口径。
第二个细节：上限截断。标准做法是把超过 30 摄氏度（部分作物用 35）的部分截掉，因为高温对作物发育的贡献不是线性的。不截断会让盛夏的积温被高估。
def gdd(tavg, tmax, tmin, base=10.0, cap=30.0): if pd.isna(tavg): return np.nan t = max(tavg, base) t = min(t, cap) return t - base tem["GDD10"] = [gdd(a, x, n) for a, x, n in zip(tem["TEM_AVG"], tem["TEM_MAX"], tem["TEM_MIN"])] tem["GDD10_CUM"] = tem.groupby(["STA", "YEAR"])["GDD10"].cumsum()第三个细节，也是最容易忽略的：跨年作物的积温窗口。冬小麦的生育期是从上一年秋季到当年夏季，如果你按自然年 groupby，积温会被切断在 12 月 31 日。正确做法是先给每个日期分配一个“农业年”标签（比如 8 月 1 日之后算下一年），再按农业年聚合。
4.3 太阳能与风能资源初筛
日照时数到太阳辐射的转换，业内常用 Ångström–Prescott 公式：
H = H0 * (a + b * n / N)其中 H 是日总辐射，H0 是天文辐射（按纬度和日序计算），n 是实际日照时数，N 是天文可照时数，a、b 是经验系数。国内不同气候区的推荐值差异不小，常见取值在 a≈0.20~0.25、b≈0.50~0.60 之间。这里要强调一点：这些系数是区域经验值，直接套用外区参数是太阳能资源评估里最常见的错误来源之一，有条件的话用邻近辐射站的实测值做一次本地标定。
天文辐射和可照时数可以自己算：
def astro_radiation(lat_deg, doy): phi = np.radians(lat_deg) Gsc = 0.0820 # 太阳常数, MJ m-2 min-1 dr = 1 + 0.033 * np.cos(2 * np.pi * doy / 365) delta = 0.409 * np.sin(2 * np.pi * doy / 365 - 1.39) ws = np.arccos(np.clip(-np.tan(phi) * np.tan(delta), -1, 1)) Ra = (24 * 60 / np.pi) * Gsc * dr * ( ws * np.sin(phi) * np.sin(delta) + np.cos(phi) * np.cos(delta) * np.sin(ws) ) N = 24 / np.pi * ws return Ra, N # MJ m-2 d-1, 小时风电初筛则绕不开风速高度订正。气象站的风速传感器一般在 10 米高度，而风机轮毂在 70 到 100 米，必须换算。幂律廓线是最常用的粗算方法：v(z) = v10 * (z / 10) ** alpha。alpha 的取值是全部学问所在——开阔平坦地区取 0.12 到 0.16，植被和建筑密集区取 0.25 到 0.30，城市下垫面甚至更高。这个参数没法从数据里直接读出来，得结合站点周边的地表覆盖情况判断。
我在实际项目里踩过的坑是：拿平均风速做了高度订正，却忘了极大风速也要订，最后在结构抗风设计那一栏给了一个明显偏低的值。平均风和极大风的廓线指数并不相同，前者更接近中性层结条件下的经验值，后者受湍流和阵风强度影响，订正时应该保守一些。
4.4 机器学习建模时的特征工程注意点
拿这份数据做特征，有三个红线不能碰。
第一条，禁止未来信息泄漏。比如你要预测某日的降水，就不能把当日的相对湿度当成特征，因为那是同一时刻的观测。所有特征必须严格来自预测时刻之前，时间窗口的切分也要按时间顺序做，不能随机切分训练集和测试集——气象序列的强自相关性会让随机切分得到一个漂亮但虚假的验证分数。
第二条，缺测必须是显式特征。“这个站今天有没有观测”本身就是一条信息。很多站点的缺测不是随机的，而是集中在恶劣天气或设备故障期间，直接插补掉会破坏这个信号。我的做法是保留一个二值的is_missing列。
第三条，空间信息要用起来。相邻站点的同日观测高度相关，把周边站点的值作为特征往往能显著提升效果。但要注意站点分布不均——东部密集、西部稀疏，做空间最近邻时如果直接用经纬度欧氏距离，高原地区的“最近邻”实际地理距离可能是平原地区的十倍。更稳妥的做法是用海拔和气候区做分层，或者干脆把经纬度做等距投影后再算距离。
5. 踩坑实录与问题速查
5.1 数据读取阶段的五个高频坑
第一个坑是编码。原始 TXT 多为本地编码，用默认 utf-8 读会报UnicodeDecodeError，或者更糟——不报错但读出乱码。处理办法是先探测编码，或者在 open 时加errors="replace"先跑通流程，再回头修正。
第二个坑是行尾空白。部分文件每行末尾带一个多余的分隔符，用固定列名读会报“期望 13 列，实际 14 列”。解决方式是加usecols限定列索引，或者读完再切掉多余列。
第三个坑是列数不一致。早年份文件的字段比晚年份少（某些要素是后来才纳入观测的），强行用同一套列名合并会报错。我的处理是按年份分组，用不同的列名模板，最后统一到一套 schema 再 concat。
第四个坑是数值类型被推断错。区站号、年份这些看着像数字的字段，如果中间混进了非数字字符，pandas 会把整列变成 object，后续做算术直接失败。显式指定dtype是唯一可靠的解法。
第五个坑是内存。全量读入几十年的文本数据，很容易吃掉几个 G 内存。按年读取、逐年转换成 parquet、最后按需读取特定年份，是我现在固定的处理流程。
5.2 分析阶段的四个隐形坑
第一个隐形坑是站点集合变化导致的伪趋势。区域平均序列如果不管站点集合，早期只有几个站、后期有几十个站，平均值的变化就会混入“观测网络扩张”这个非气候因素。标准做法是只保留全时段连续观测的站点，或者用固定站点子集（比如要求每年至少有 330 天有效观测）算区域平均。
第二个隐形坑是0.1 毫米的降水阈值。日降水量小于 0.1 毫米的记录，在原始数据里可能是 0，也可能是微量降水被记为 0。做“降水日数”统计时，如果你用“降水量大于 0”做判据，结果会比用“大于等于 0.1 毫米”高出不少。哪个对？取决于你的业务定义，但必须明确写出来。
第三个隐形坑是日照时数与理论可照时数的混淆。日照时数最大不可能超过天文可照时数，但在高纬度地区，夏季可照时数接近 16 小时，而部分站点的记录能到 12 小时以上，看起来没问题，实际上可能是感应器灵敏度阈值造成的。做资源评估时，最好把 n/N 的比值画个分布图看一眼，正常应该集中在 0.3 到 0.7 之间，出现大量接近 1 的值就要怀疑。
第四个隐形坑是缺测率的时间不均匀。早年观测手段落后，缺测率高；近年自动化程度高，缺测率低。如果直接对全时段求平均，等于给早年数据更小的权重（因为有效样本少），这在统计上未必错，但解释时必须说清楚。
5.3 问题速查表
| 现象 | 最可能的原因 | 处理办法 |
|---|---|---|
| 均值明显偏离常识（如气温 3000 度） | 缺测码未替换 | 读取时用 na_values 纳入 32766 及 999xxx |
| merge 后行数为 0 | 区站号类型不一致或前导零丢失 | 统一转字符串并去空格 |
| 某些年份整年缺失 | 该年文件损坏或未下载 | 核对文件清单与文件大小 |
| 站点在地图上跑到境外 | 经纬度为度分格式未换算 | 按度分换算后做空间核查 |
| 气温序列出现台阶跳变 | 迁站 | 查元数据迁站记录，做均一性检验 |
| 极端值大量被剔除 | 质控策略过严 | 极端统计场景保留质控码 1 |
| 积温在 12 月 31 日断裂 | 按自然年聚合 | 改用农业年标签聚合 |
| 降水日数偏高 | 阈值用了大于 0 | 改为大于等于 0.1 毫米 |
| 读取速度极慢 | 每次重复解析文本 | 一次性转 parquet |
| 内存溢出 | 全量读入 | 按年分块处理 |
6. 几条个人经验
这份数据我用得最久，也栽得最多。最想说的是：别急着建模，先把数据本身摸清楚。我见过太多人拿到数据直接扔进模型，跑出一个 R² 0.9 的结果，然后发现测试集里混了训练年份的数据。气象序列的时间自相关性极强，任何一次切分失误都会被它放大成虚假的性能。
第二个体会是，元数据的重要性不亚于观测值本身。站点搬迁、仪器更换、观测时次调整，这些信息全在元数据里，而它们对长序列的影响远大于随机观测误差。我现在的习惯是，任何一份新拿到的气象数据，第一步永远是读元数据表，把站点数量、建站年份、迁站记录、缺测率画成四张图看一遍，看完再动观测值。
第三个体会是，质控码要当朋友，不要当敌人。新手看到一列满是 0 和 1 的码，第一反应是“这列没用，删了”。实际上这列是数据发布方留给你的一份“体检报告”，它告诉你哪些值需要多看一眼。真正专业的做法不是无脑剔除，而是分类处理，并且在报告里写清楚你的处理策略——这既是严谨性，也是让别人能复现你结果的前提。
最后一个技巧，如果你要处理多个要素、几十年的数据，先写一个能跑通单年单要素的最小脚本，跑通了再横向扩展。我吃过最亏的一次就是直接写了个大而全的管道，跑到一半报错，因为前面几百个文件已经处理过一部分，状态混乱，最后只能全部推倒重来。小步快跑，每步验证，在这类长周期数据处理任务上，是最省时间的做法。
