简介一份面向数据分析初学者及气象数据爱好者的完整项目资料包基于中国天气网某城市历史天气数据进行全流程分析。项目提供Python爬虫源代码可自动抓取气温、湿度、风力和空气质量等字段并支持在Jupyter Notebook中直接运行后续通过数据清洗与特征统计生成雷达图、条形图等可视化图表对城市天气变化规律及指标间关系进行直观解读。压缩包共39个文件包含20张PNG图表输出以及15个XML和4个rels文件——它们构成实验报告Word文档的内部结构整体仅1.69MB轻便易获取。值得注意的是资料中附带整理完成的实验报告详细记录从爬虫采集到图表分析的每一步思路读者可对照源码复现结果亦可借此掌握天气数据获取、Pandas处理与Matplotlib绘图的完整链路。当前已有1853人学习下载尤其适合用于数据分析课程设计、期末实验或入门实战参考。1. 气象分析到底在分析什么一份气温数据能挖出什么拿到十年的逐日气象数据第一件事不是画图而是先问一句这份数据能不能支撑我要下的结论。气象分析数据分析这个方向在从业者手里通常指把观测站或再分析资料变成可量化的气候结论——温度趋势、降水变化、极端事件频率、不同站点之间的差异而不是做天气预报。两者的差别很实际预报看未来三天分析看过去十年到三十年发生了什么、正在怎么变。它能解决的具体问题分布在很多行业里。农业上算积温和无霜期决定播种窗口能源行业算制冷度日和采暖度日用来估夏季用电峰值物流和保险行业则关心极端降水、大风天气的出现频率用来定风险敞口。适合做这件事的人一类是数据分析师想往行业方向深挖另一类是气象相关业务人员想把观测数据转成决策依据。我一般把这类项目拆成五步走数据源选型、清洗、时间序列与空间统计、可视化、结论验证。前两步决定后面所有分析的地基后面每一步都有实实在在的坑。这篇就按这条路径展开每步都用能直接跑的代码说明白。2. 先解决数据问题气象数据有哪些来源CSV与NetCDF格式怎么选2.1 三种常见数据源观测站、再分析资料、气象API做气象分析数据源的选择几乎决定了分析的上限。观测站数据来自国家气象信息中心、NOAA GHCN这类机构特点是站点的实测记录精度最高但站点分布不均——东部密、西部稀而且部分站点序列有断点。再分析资料最常见的是欧洲中期天气预报中心的ERA5和NCEP再分析资料它们是模式与观测同化出来的全球网格数据空间上覆盖完整、时间序列长适合做区域尺度的气候诊断但单点数值和实测有偏差。第三类是气象API比如Open-Meteo这类免费接口适合快速做原型验证不需要下载文件但历史深度和字段完整度有限。数据源空间粒度时间覆盖适合场景获取成本观测站数据站点点位数十年逐日/逐时单站或多站对比、趋势分析需注册申请审核后下载再分析资料全球网格0.25°~1°数十年逐小时/逐日区域空间分析、缺测补全、气候诊断免费数据量大气象API站点或网格受接口限制原型验证、快速出图免费额度有限流选型逻辑很简单做地方尺度的分析优先观测站做空间分布和区域对比用再分析资料项目初期摸流程先用API把代码跑通。常见做法是混用——观测站数据为主再分析资料做交叉验证这样结论才经得起追问。2.2 拿到数据第一步清洗缺失值与统一时间格式气象原始数据永远不会是干净的直接可用状态。以中国气象数据网的站点逐日数据为例文件里常见三类坑缺失值用特定数值填充32744、32766这类气温单位是0.1℃降水单位是0.1mm。如果不知道这些约定直接读后面的趋势分析全是错的。import pandas as pd # 读取站点日值数据na_values把特定数值识别为缺失 df pd.read_csv( station_daily.csv, na_values[32744, 32766, 9999, -9999], usecols[station, date, tem, pre, wind], ) # 中国站点数据中气温和降水的原始单位是0.1℃和0.1mm df[tem] df[tem] / 10.0 df[pre] df[pre] / 10.0 # 统一时间列构建DatetimeIndex df[date] pd.to_datetime(df[date], format%Y%m%d) df df.set_index(date).sort_index() print(df.head()) print(df.isna().sum())这段代码做了三件事缺测识别、单位换算、时间索引构建。na_values参数是关键它把气象数据里常见的填充缺测值统一映射为NaN后面所有统计天然忽略缺失。单位换算藏得很深——如果没有看过数据说明文档做出来的温度曲线会是真实值的10倍这种错误在数据量大时很难肉眼看出来。to_datetime的时间格式参数%Y%m%d对应当天日期列的原始字符串如果日期列是20240101这种八位数字这个格式刚好匹配。2.3 把NetCDF网格数据变成可分析的表格再分析资料的NetCDF文件没法直接用Pandas打开需要一个转换步骤。常见做法是用xarray读取后再选点或区域聚合转成DataFrame做后续分析。这一步的难点不在API调用而在理解NetCDF的维度结构——通常三维是(time, lat, lon)但不同数据集的纬度顺序和单位不一样ERA5的经纬度单位是度温度单位是开尔文这些都要先确认。import xarray as xr import pandas as pd # 打开再分析资料NetCDF文件 ds xr.open_dataset(era5_single_level.nc) print(ds.coords) # 确认维度名称和范围 # 用最近邻方法提取目标站点所在格点北京约39.9N, 116.4E lat_target, lon_target 39.9, 116.4 # 低温转摄氏开尔文减273.15 ds[t2m_c] ds[t2m] - 273.15 # 选择时间范围示例取2020年 df_series ( ds[t2m_c] .sel(timeslice(2020-01-01, 2020-12-31), latlat_target, lonlon_target, methodnearest) .to_dataframe() .drop(columns[lat, lon]) .dropna() ) print(df_series.head())这段代码的核心是sel加methodnearest它能自动找到离目标经纬度最近的网格点避免手工找索引。需要留意的是再分析资料的默认时间是UTC和本地观测站的时间系统差8小时如果你要对比日值数据得先做时区平移第5章有详细说明。dropna用于把海陆边界上落在海里的网格点剔除这些位置没有有效观测同化数值是填充值。3. 用Pandas做温度时间序列分析气候态、距平与趋势3.1 构建时间序列与重采样日值转月值、年值温度数据按日存储时噪声很大直接看趋势线会被逐日波动淹没。标准做法是先按月和按年重采样用月均值和年均值来观察气候尺度的变化。Pandas的resample在这里有两个必须注意的细节新版本2.2推荐用字符串 ME 表示按月取月末旧版的 M 已经被标记弃用按年聚合用 YE同样是为了消歧。# 日值转月均值和年均值 monthly_mean df[tem].resample(ME).mean() annual_mean df[tem].resample(YE).mean() # 查看重采样前后的数据量变化 print(f原始日值数量: {df[tem].shape[0]}) print(f月均值数量: {monthly_mean.shape[0]}) print(f年均值数量: {annual_mean.shape[0]}) # 把年均值保存为DataFrame annual_df annual_mean.to_frame(nameannual_tem) annual_df[year] annual_df.index.yearresample(ME)的语义是月末Month End聚合它会把每个月所有日值取平均同时自动处理各月天数不同的问题——2月只有28天也不会被当成30天算。这个细节比groupby(df.index.month)要可靠得多后者不会自动对齐日历边界遇到跨年数据时容易出错。重采样之后的数据量大幅缩减后续回归和可视化都基于月均或年均趋势信号才不会被压制。3.2 气候态均值与距平为什么拿30年做基线气候分析里有个基础概念叫气候态均值简单说就是一段足够长的历史时期的平均状态。气象上通常用30年作为标准基线WMO推荐的经典时段是1981-2010现在越来越多人切换到1991-2020。用这个基线算出来的差值叫距平距平为正代表偏暖为负代表偏冷。比起绝对温度距平更适合做趋势分析因为它剔除了站点海拔和纬度带来的本底差异。# 以1991-2020为气候态基线计算逐日气候态均值 clim_start, clim_end 1991, 2020 clim_mask (df.index clim_start) (df.index clim_end) clim_data df.loc[clim_mask, tem] # 按一年中的第几天分组计算逐日气候态 daily_clim clim_data.groupby(clim_data.index.dayofyear).mean() # 计算全序列每日距平 df[anomaly] df[tem] - df.index.dayofyear.map(daily_clim) # 月距平和年距平 monthly_anom df[anomaly].resample(ME).mean() annual_anom df[anomaly].resample(YE).mean() print(annual_anom.tail())关键点在于用dayofyear分组计算逐日气候态这样能保留季节内每一天的基准值而不是每个月用一个粗粒度平均值。如果只用月气候态1月上旬和下旬的温差会被抹平距平序列里会混入季节残留信号。基线期的选择直接影响距平符号——同一年的1月用1981-2010算可能是正距平切到1991-2020可能就变成负距平这是正常现象报告里必须写清楚用哪个时段。3.3 趋势分析与可视化斜率怎么算才不被夏季波动干扰算温度趋势最直接的方法是线性回归以年份为自变量、年距平为因变量斜率就是每十年的变化量。但这里有个统计陷阱如果直接用月距平做回归相邻月份之间的自相关会让趋势的显著性被高估。常见做法是回归前先聚合成年均值牺牲一些样本量换来自变量的独立性。from scipy.stats import linregress import matplotlib.pyplot as plt # 准备回归数据 analysis_df annual_anom.dropna().to_frame(nameanom) analysis_df[time_idx] range(len(analysis_df)) # 线性回归 slope, intercept, r_value, p_value, std_err linregress( analysis_df[time_idx], analysis_df[anom] ) # 输出每十年变化量 years_per_decade 10 trend_per_decade slope * years_per_decade print(f趋势: {trend_per_decade:.2f} ℃/10年, p值: {p_value:.4f}) # 绘图 plt.figure(figsize(10, 5)) plt.plot(analysis_df.index, analysis_df[anom], label年距平, colorblack, lw0.8) plt.plot(analysis_df.index, intercept slope * analysis_df[time_idx], labelf线性趋势 {trend_per_decade:.2f} ℃/10年, colorred, lw2) plt.axhline(0, colorgray, ls--, lw0.5) plt.legend() plt.xlabel(年份) plt.ylabel(温度距平 (℃)) plt.title(年均温度距平与线性趋势) plt.show()linregress返回的slope是年平均距平随序号变化的速率乘10就是每十年的变化量。p_value用于判断趋势是否统计显著小于0.05通常被认为不是随机波动造成的。绘图时把距平序列和回归线叠在一起能直观看出趋势是否稳定——如果序列前30年平稳、后10年急剧上升单一线性趋势会把这段转折掩盖掉。遇到这种形态更好的做法是分时段回归而不是依赖全序列一个斜率。4. 降水与极端事件分析为什么不能照搬温度的套路4.1 降水的统计特性零膨胀与偏态分布降水数据和温度数据在统计性质上完全不是一回事。温度接近正态分布均值有物理意义而降水是零膨胀的偏态分布——大部分日子不下雨数值为0少数日子出现大值均值会被几场极端暴雨拉高不能代表典型状态。如果把降水当温度那样求均值、画趋势结论会很离谱。# 查看降水的分布特征 pre df[pre].dropna() # 基本统计量 stats pre.describe() print(stats) # 中位数和众数 median_pre pre.median() zero_ratio (pre 0).mean() print(f降水日比例0.1mm: {(pre 0.1).mean():.2%}) print(f中位数: {median_pre:.1f} mm) # 日降水直方图对数y轴更直观 import matplotlib.pyplot as plt fig, ax plt.subplots() ax.hist(pre[pre 0.1], bins50, logTrue, colorsteelblue) ax.set_xlabel(日降水量 (mm)) ax.set_ylabel(天数 (对数坐标)) ax.set_title(降水日降水量分布剔除无雨日) plt.show()describe给出的均值和中位数如果差距很大就已经提示数据是偏态的。zero_ratio告诉你无雨日占比中国北方很多站点超过70%。画直方图时用了logTrue否则大部分柱子都挤在0-5mm区间根本看不出极端降水尾巴的形状。对这类数据描述统计应改用降水日数、百分位数、超过阈值的频次而不是均值。4.2 降水日与极端降水阈值0.1mm和95百分位气象行业对降水日有明确定义日降水量大于等于0.1mm算一个降水日。极端降水事件则常用百分位阈值来界定——把历史所有降水日的降水量从小到大排序取第95百分位作为极端阈值凡是超过这个值的日子就算极端降水事件。这个方法比固定阈值比如50mm更合理因为不同气候区的降水底数差异巨大西北站点年降水总量可能不如东南一场暴雨用统一阈值会漏掉区域特征。# 定义降水日与极端降水阈值 pre_days pre[pre 0.1] # 95百分位阈值仅基于降水日 threshold_95 pre_days.quantile(0.95) print(f极端降水阈值95百分位: {threshold_95:.1f} mm) # 统计每年极端降水日数 df[is_extreme] df[pre] threshold_95 annual_extreme_days df[is_extreme].resample(YE).sum() # 年降水总量 annual_pre df[pre].resample(YE).sum() print(annual_extreme_days.tail())用quantile(0.95)代替固定阈值可以让不同站的极端事件定义在统计上可比。需要留意的是计算百分位时只用降水日样本而不是全部天数否则0值会把阈值压得很低导致「极端事件」里混入普通降雨。输出每年的极端降水日数后可以看到它的年际波动通常比温度大得多——降水本来就是高变异变量一年的极端日数从3天跳到10天不罕见。4.3 连续干期与雨季集中度用滑动窗口量化干旱与汛期除了极端降水事件农业和水利上更关心的是无雨持续时间和降水集中在哪个季节。连续无雨日是抗旱决策的直接指标而降水集中度决定了水库调度策略。两者都可以用简单统计实现不需要复杂的干旱指数模型。# 计算每年最大连续无雨日数 pre_binary (df[pre] 0.1).astype(int) # 无雨日标记为1 max_dry_spell {} for year, group in pre_binary.groupby(pre_binary.index.year): # 按连续为1的块分组取最大块长度 spell group.groupby((group ! group.shift()).cumsum()).sum() max_dry_spell[year] spell.max() if len(spell) 0 else 0 # 计算雨季5-9月降水占全年比例 seasonal df.groupby(df.index.month)[pre].sum() rainy_season_ratio ( seasonal.loc[5:9].sum() / seasonal.sum() ) print(f汛期5-9月降水占比: {rainy_season_ratio:.1%})连续无雨日的计算用了shift加cumsum的分组技巧——当天的无雨状态与前一日不同时分组编号加1从而把连续的无雨日块切分出来。max_dry_spell的结果直接反映季节性干旱风险。雨季集中度则简单得多按月汇总后算5到9月占比华北很多站点这个比例超过70%意味着汛期之外的半年几乎无降水。5. 气象数据分析避坑五个高频翻车点与排查方法5.1 时区陷阱UTC与北京时间的8小时偏移现象把再分析资料的日降水和观测站数据放在同一张图里对比发现降水事件总是对不上甚至日期都差了一天。原因再分析资料ERA5、NCEP时间基准是UTC而国内观测站数据通常是北京时间两者相差8小时。日值数据如果按自然日聚合UTC的一天从北京时间早8点开始夜间到早晨的降水会被划分到错误的日期。解决在读取网格数据后先把时间索引转成目标时区再按日重采样。# 把UTC时间转为北京时间UTC8再按日聚合 df_utc df_era5.copy() df_utc.index df_utc.index.tz_localize(UTC).tz_convert(Asia/Shanghai) df_daily_utc df_utc.resample(D).mean()tz_localize先给无时区的时间戳标记为UTCtz_convert再转到北京时间。如果原始数据本身有时区信息直接tz_convert即可。这个转换必须在所有聚合操作之前做顺序反了会造成无法挽回的偏移。5.2 缺测值不一定是NaN32744、9999、-9999现象读进来的数据没有NaN但画趋势图时出现一个离谱的尖峰比如7月气温突然变成9999℃。原因气象数据文件早期用特定数值表示缺测不同机构约定不同常见的有32744、32766、9999、-9999这些值没有在CSV里标记为缺失被当成真实观测读入。解决读文件前先打开原始数据说明把缺测值通过na_values传入或者读进来后用条件筛选统一替换。# 检查有没有未识别的缺测值 for col in df.columns: for bad_val in [32744, 32766, 9999, -9999]: n_bad (df[col] bad_val).sum() if n_bad 0: print(f列 {col} 中发现缺测值 {bad_val}数量 {n_bad}) # 统一替换为NaN df df.replace([32744, 32766, 9999, -9999], np.nan)这段排查代码建议在清洗流程里固定跑一次特别是拿到新数据源时。气象站数据的说明文档经常藏在下载包的readme里里面会写明缺测值编码和单位信息花两分钟读一遍能省一整天排错时间。5.3 闰年与2月29日重采样的隐蔽多一天现象按年对比各月均值时每年2月的平均值出现小幅但系统的偏差某几年数值明显不同。原因2月有28天和29天两种长度如果按自然月重采样日值数量不同导致月均值权重变化更隐蔽的是dayofyear在闰年3月之后会比平年多一天导致逐日气候态对齐错位。解决月重采样用resample(ME)可以避免多数问题计算逐日气候态时对闰年的2月29日做单独处理。# 剔除2月29日保证逐日气候态对齐 df_clean df[ ~((df.index.month 2) (df.index.day 29)) ] # 重新按dayofyear算气候态 daily_clim df_clean.groupby(df_clean.index.dayofyear).mean()去掉2月29日只损失一个样本对月均值和气候态的影响微乎其微但能消除每年多一天带来的对齐偏差。如果做的是逐日气候态产品发布这个步骤特别关键否则3月1日之后每一天的气候态基准都会错位一天。5.4 网格数据取点最近邻不是「取整」现象用sel(lat39.9, lon116.4)从NetCDF里取北京的气温和站点实测对比相关系数很高但绝对偏差一直存在而且偏差在不同季节大小不一。原因ERA5的经纬度网格是0.25°间隔格点坐标可能是39.75、40.0这种值直接传入的39.9并不落在格点上如果手写取整逻辑round到最近的0.25可能选到距离目标几十公里外的格点山地和沿海站点偏差尤其明显。解决用sel(..., methodnearest)它内部按距离选择最近格点更严格的场景先算经纬度距离矩阵再选最近邻。# 正确方式最近邻选取 ts_point ds[t2m].sel( lat39.9, lon116.4, methodnearest ) # 严谨方式先确认实际选中的格点坐标 print(f选中格点: lat{ts_point.lat.values:.2f}, lon{ts_point.lon.values:.2f})打印选中格点的坐标是个好习惯——它告诉你实际用的数据来自哪个网格点而不是想当然认为就是目标位置。地形复杂区域如果最近格点与实际站点海拔差超过500米气温偏差会达到3℃以上这时候应该改用双线性插值而不是最近邻。5.5 单位陷阱气温0.1℃、降水0.1mm现象画出的温度曲线在30到40℃之间剧烈震荡年降水量动辄上万毫米明显超出常识。原因国内很多气象数据为节省存储空间采用缩小10倍的整数编码——气温原始值代表0.1℃降水代表0.1mm风速代表0.1m/s。没有按说明除以10所有分析结果都会系统性放大10倍。解决在数据加载层统一做单位换算不要在使用时临时处理。# 加载时统一单位换算 unit_map {tem: 0.1, pre: 0.1, wind: 0.1} for col, factor in unit_map.items(): if col in df.columns: df[col] df[col] * factor这个坑几乎每个用过中国气象站点数据的人都踩过。我现在的习惯是写一个标准加载函数数据一进内存就完成单位换算和缺失值处理任何下游分析都不会再碰到原始值。另一个值得做的检查是画图后先看一眼纵轴范围——如果温度范围在-30到45℃之外第一反应不是分析问题而是单位或缺测处理出了问题。6. 进阶让气象分析结论更稳的三个实操技巧6.1 滑动平均去噪窗口宽度怎么定逐日温度曲线噪声大逐日降水曲线噪声更大。滑动平均是简单有效的去噪方式但窗口宽度要根据你想保留的信号尺度来选择7天窗口能去掉天气尺度波动保留一周以上的天气过程30天窗口适合看月际背景如果想看年际信号直接用年均值比任何滑动平均都干净。窗口宽度没有统一标准我在报告里通常同时画原始序列和两种窗口的滑动平均线让读者自己判断。df[tem_smooth_7d] df[tem].rolling(window7, centerTrue).mean() df[tem_smooth_30d] df[tem].rolling(window30, centerTrue).mean()centerTrue让滑动平均窗口以当前天为中心而不是只取过去这样平滑后的序列不会整体滞后。注意到窗口开头和结尾会有NaN这是边界效应绘图时可选择从有值的位置开始。6.2 相关分析的样本量陷阱气象时间序列的自相关统计两个气象变量比如温度和湿度的Pearson相关系数时如果直接把逐日数据全丢进去n值看起来很大十年有3650个样本但气象时间序列存在强自相关——今天的温度和昨天几乎一样有效样本量远小于实际样本量。直接按n3650查显著性表会把不相关的变量判成显著相关。修正办法是用有效样本量计算自由度。def effective_n(series1, series2): 计算考虑自相关后的有效样本量 n len(series1) r1 series1.autocorr(lag1) r2 series2.autocorr(lag1) if abs(r1) 1 or abs(r2) 1: return n n_eff n * (1 - r1 * r2) / (1 r1 * r2) return int(n_eff)使用autocorr(lag1)得到滞后1阶自相关系数代入修正公式。常见替代方案是先把数据聚合成月距平再做相关月距平序列的自相关大幅降低有效样本量的问题自然缓解。在报告里标注有效样本量是严谨性的加分项也是很多论文质检点。6.3 验证结论的两个办法分段时间对比与交叉数据任何趋势或相关结论在发布前都应该经过验证。第一个办法是分段稳定性检验把序列按时间切成前后两半分别计算趋势如果两半的趋势符号不同说明全序列的趋势可能是气候变率噪声而不是稳定趋势。第二个办法是交叉数据验证用再分析资料算同一个站点的趋势或极端事件频率与观测站结果对比符号和量级一致才说明结论不是单一数据源的系统误差造成的。这两个验证方法不需要额外代码复用前面章节的回归和统计函数换数据切片就行。我的习惯是把所有分析封装成以DataFrame为输入的函数验证时只需要传入不同时段或不同来源的数据。有一次我在一个站点上算出了显著的增温趋势分段检验发现前半段几乎无趋势、后半段急剧上升这个信息比单一斜率值重要得多——它说明了变化发生的时段而不是笼统一句「在变暖」。气象数据分析做久了会发现大部分翻车不是统计模型不够高级而是数据在进入模型之前就错了。我现在的项目里固定保留一个数据校验层单位、时区、缺测、闰年四道检查全部跑完才开始分析这套习惯的养成比任何算法都值钱。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/10 19:55:40 · 阅读 1.2k Python气象数据分析实战：从数据清洗到温度与降水趋势提取
简介一份面向数据分析初学者及气象数据爱好者的完整项目资料包基于中国天气网某城市历史天气数据进行全流程分析。项目提供Python爬虫源代码可自动抓取气温、湿度、风力和空气质量等字段并支持在Jupyter Notebook中直接运行后续通过数据清洗与特征统计生成雷达图、条形图等可视化图表对城市天气变化规律及指标间关系进行直观解读。压缩包共39个文件包含20张PNG图表输出以及15个XML和4个rels文件——它们构成实验报告Word文档的内部结构整体仅1.69MB轻便易获取。值得注意的是资料中附带整理完成的实验报告详细记录从爬虫采集到图表分析的每一步思路读者可对照源码复现结果亦可借此掌握天气数据获取、Pandas处理与Matplotlib绘图的完整链路。当前已有1853人学习下载尤其适合用于数据分析课程设计、期末实验或入门实战参考。1. 气象分析到底在分析什么一份气温数据能挖出什么拿到十年的逐日气象数据第一件事不是画图而是先问一句这份数据能不能支撑我要下的结论。气象分析数据分析这个方向在从业者手里通常指把观测站或再分析资料变成可量化的气候结论——温度趋势、降水变化、极端事件频率、不同站点之间的差异而不是做天气预报。两者的差别很实际预报看未来三天分析看过去十年到三十年发生了什么、正在怎么变。它能解决的具体问题分布在很多行业里。农业上算积温和无霜期决定播种窗口能源行业算制冷度日和采暖度日用来估夏季用电峰值物流和保险行业则关心极端降水、大风天气的出现频率用来定风险敞口。适合做这件事的人一类是数据分析师想往行业方向深挖另一类是气象相关业务人员想把观测数据转成决策依据。我一般把这类项目拆成五步走数据源选型、清洗、时间序列与空间统计、可视化、结论验证。前两步决定后面所有分析的地基后面每一步都有实实在在的坑。这篇就按这条路径展开每步都用能直接跑的代码说明白。2. 先解决数据问题气象数据有哪些来源CSV与NetCDF格式怎么选2.1 三种常见数据源观测站、再分析资料、气象API做气象分析数据源的选择几乎决定了分析的上限。观测站数据来自国家气象信息中心、NOAA GHCN这类机构特点是站点的实测记录精度最高但站点分布不均——东部密、西部稀而且部分站点序列有断点。再分析资料最常见的是欧洲中期天气预报中心的ERA5和NCEP再分析资料它们是模式与观测同化出来的全球网格数据空间上覆盖完整、时间序列长适合做区域尺度的气候诊断但单点数值和实测有偏差。第三类是气象API比如Open-Meteo这类免费接口适合快速做原型验证不需要下载文件但历史深度和字段完整度有限。数据源空间粒度时间覆盖适合场景获取成本观测站数据站点点位数十年逐日/逐时单站或多站对比、趋势分析需注册申请审核后下载再分析资料全球网格0.25°~1°数十年逐小时/逐日区域空间分析、缺测补全、气候诊断免费数据量大气象API站点或网格受接口限制原型验证、快速出图免费额度有限流选型逻辑很简单做地方尺度的分析优先观测站做空间分布和区域对比用再分析资料项目初期摸流程先用API把代码跑通。常见做法是混用——观测站数据为主再分析资料做交叉验证这样结论才经得起追问。2.2 拿到数据第一步清洗缺失值与统一时间格式气象原始数据永远不会是干净的直接可用状态。以中国气象数据网的站点逐日数据为例文件里常见三类坑缺失值用特定数值填充32744、32766这类气温单位是0.1℃降水单位是0.1mm。如果不知道这些约定直接读后面的趋势分析全是错的。import pandas as pd # 读取站点日值数据na_values把特定数值识别为缺失 df pd.read_csv( station_daily.csv, na_values[32744, 32766, 9999, -9999], usecols[station, date, tem, pre, wind], ) # 中国站点数据中气温和降水的原始单位是0.1℃和0.1mm df[tem] df[tem] / 10.0 df[pre] df[pre] / 10.0 # 统一时间列构建DatetimeIndex df[date] pd.to_datetime(df[date], format%Y%m%d) df df.set_index(date).sort_index() print(df.head()) print(df.isna().sum())这段代码做了三件事缺测识别、单位换算、时间索引构建。na_values参数是关键它把气象数据里常见的填充缺测值统一映射为NaN后面所有统计天然忽略缺失。单位换算藏得很深——如果没有看过数据说明文档做出来的温度曲线会是真实值的10倍这种错误在数据量大时很难肉眼看出来。to_datetime的时间格式参数%Y%m%d对应当天日期列的原始字符串如果日期列是20240101这种八位数字这个格式刚好匹配。2.3 把NetCDF网格数据变成可分析的表格再分析资料的NetCDF文件没法直接用Pandas打开需要一个转换步骤。常见做法是用xarray读取后再选点或区域聚合转成DataFrame做后续分析。这一步的难点不在API调用而在理解NetCDF的维度结构——通常三维是(time, lat, lon)但不同数据集的纬度顺序和单位不一样ERA5的经纬度单位是度温度单位是开尔文这些都要先确认。import xarray as xr import pandas as pd # 打开再分析资料NetCDF文件 ds xr.open_dataset(era5_single_level.nc) print(ds.coords) # 确认维度名称和范围 # 用最近邻方法提取目标站点所在格点北京约39.9N, 116.4E lat_target, lon_target 39.9, 116.4 # 低温转摄氏开尔文减273.15 ds[t2m_c] ds[t2m] - 273.15 # 选择时间范围示例取2020年 df_series ( ds[t2m_c] .sel(timeslice(2020-01-01, 2020-12-31), latlat_target, lonlon_target, methodnearest) .to_dataframe() .drop(columns[lat, lon]) .dropna() ) print(df_series.head())这段代码的核心是sel加methodnearest它能自动找到离目标经纬度最近的网格点避免手工找索引。需要留意的是再分析资料的默认时间是UTC和本地观测站的时间系统差8小时如果你要对比日值数据得先做时区平移第5章有详细说明。dropna用于把海陆边界上落在海里的网格点剔除这些位置没有有效观测同化数值是填充值。3. 用Pandas做温度时间序列分析气候态、距平与趋势3.1 构建时间序列与重采样日值转月值、年值温度数据按日存储时噪声很大直接看趋势线会被逐日波动淹没。标准做法是先按月和按年重采样用月均值和年均值来观察气候尺度的变化。Pandas的resample在这里有两个必须注意的细节新版本2.2推荐用字符串 ME 表示按月取月末旧版的 M 已经被标记弃用按年聚合用 YE同样是为了消歧。# 日值转月均值和年均值 monthly_mean df[tem].resample(ME).mean() annual_mean df[tem].resample(YE).mean() # 查看重采样前后的数据量变化 print(f原始日值数量: {df[tem].shape[0]}) print(f月均值数量: {monthly_mean.shape[0]}) print(f年均值数量: {annual_mean.shape[0]}) # 把年均值保存为DataFrame annual_df annual_mean.to_frame(nameannual_tem) annual_df[year] annual_df.index.yearresample(ME)的语义是月末Month End聚合它会把每个月所有日值取平均同时自动处理各月天数不同的问题——2月只有28天也不会被当成30天算。这个细节比groupby(df.index.month)要可靠得多后者不会自动对齐日历边界遇到跨年数据时容易出错。重采样之后的数据量大幅缩减后续回归和可视化都基于月均或年均趋势信号才不会被压制。3.2 气候态均值与距平为什么拿30年做基线气候分析里有个基础概念叫气候态均值简单说就是一段足够长的历史时期的平均状态。气象上通常用30年作为标准基线WMO推荐的经典时段是1981-2010现在越来越多人切换到1991-2020。用这个基线算出来的差值叫距平距平为正代表偏暖为负代表偏冷。比起绝对温度距平更适合做趋势分析因为它剔除了站点海拔和纬度带来的本底差异。# 以1991-2020为气候态基线计算逐日气候态均值 clim_start, clim_end 1991, 2020 clim_mask (df.index clim_start) (df.index clim_end) clim_data df.loc[clim_mask, tem] # 按一年中的第几天分组计算逐日气候态 daily_clim clim_data.groupby(clim_data.index.dayofyear).mean() # 计算全序列每日距平 df[anomaly] df[tem] - df.index.dayofyear.map(daily_clim) # 月距平和年距平 monthly_anom df[anomaly].resample(ME).mean() annual_anom df[anomaly].resample(YE).mean() print(annual_anom.tail())关键点在于用dayofyear分组计算逐日气候态这样能保留季节内每一天的基准值而不是每个月用一个粗粒度平均值。如果只用月气候态1月上旬和下旬的温差会被抹平距平序列里会混入季节残留信号。基线期的选择直接影响距平符号——同一年的1月用1981-2010算可能是正距平切到1991-2020可能就变成负距平这是正常现象报告里必须写清楚用哪个时段。3.3 趋势分析与可视化斜率怎么算才不被夏季波动干扰算温度趋势最直接的方法是线性回归以年份为自变量、年距平为因变量斜率就是每十年的变化量。但这里有个统计陷阱如果直接用月距平做回归相邻月份之间的自相关会让趋势的显著性被高估。常见做法是回归前先聚合成年均值牺牲一些样本量换来自变量的独立性。from scipy.stats import linregress import matplotlib.pyplot as plt # 准备回归数据 analysis_df annual_anom.dropna().to_frame(nameanom) analysis_df[time_idx] range(len(analysis_df)) # 线性回归 slope, intercept, r_value, p_value, std_err linregress( analysis_df[time_idx], analysis_df[anom] ) # 输出每十年变化量 years_per_decade 10 trend_per_decade slope * years_per_decade print(f趋势: {trend_per_decade:.2f} ℃/10年, p值: {p_value:.4f}) # 绘图 plt.figure(figsize(10, 5)) plt.plot(analysis_df.index, analysis_df[anom], label年距平, colorblack, lw0.8) plt.plot(analysis_df.index, intercept slope * analysis_df[time_idx], labelf线性趋势 {trend_per_decade:.2f} ℃/10年, colorred, lw2) plt.axhline(0, colorgray, ls--, lw0.5) plt.legend() plt.xlabel(年份) plt.ylabel(温度距平 (℃)) plt.title(年均温度距平与线性趋势) plt.show()linregress返回的slope是年平均距平随序号变化的速率乘10就是每十年的变化量。p_value用于判断趋势是否统计显著小于0.05通常被认为不是随机波动造成的。绘图时把距平序列和回归线叠在一起能直观看出趋势是否稳定——如果序列前30年平稳、后10年急剧上升单一线性趋势会把这段转折掩盖掉。遇到这种形态更好的做法是分时段回归而不是依赖全序列一个斜率。4. 降水与极端事件分析为什么不能照搬温度的套路4.1 降水的统计特性零膨胀与偏态分布降水数据和温度数据在统计性质上完全不是一回事。温度接近正态分布均值有物理意义而降水是零膨胀的偏态分布——大部分日子不下雨数值为0少数日子出现大值均值会被几场极端暴雨拉高不能代表典型状态。如果把降水当温度那样求均值、画趋势结论会很离谱。# 查看降水的分布特征 pre df[pre].dropna() # 基本统计量 stats pre.describe() print(stats) # 中位数和众数 median_pre pre.median() zero_ratio (pre 0).mean() print(f降水日比例0.1mm: {(pre 0.1).mean():.2%}) print(f中位数: {median_pre:.1f} mm) # 日降水直方图对数y轴更直观 import matplotlib.pyplot as plt fig, ax plt.subplots() ax.hist(pre[pre 0.1], bins50, logTrue, colorsteelblue) ax.set_xlabel(日降水量 (mm)) ax.set_ylabel(天数 (对数坐标)) ax.set_title(降水日降水量分布剔除无雨日) plt.show()describe给出的均值和中位数如果差距很大就已经提示数据是偏态的。zero_ratio告诉你无雨日占比中国北方很多站点超过70%。画直方图时用了logTrue否则大部分柱子都挤在0-5mm区间根本看不出极端降水尾巴的形状。对这类数据描述统计应改用降水日数、百分位数、超过阈值的频次而不是均值。4.2 降水日与极端降水阈值0.1mm和95百分位气象行业对降水日有明确定义日降水量大于等于0.1mm算一个降水日。极端降水事件则常用百分位阈值来界定——把历史所有降水日的降水量从小到大排序取第95百分位作为极端阈值凡是超过这个值的日子就算极端降水事件。这个方法比固定阈值比如50mm更合理因为不同气候区的降水底数差异巨大西北站点年降水总量可能不如东南一场暴雨用统一阈值会漏掉区域特征。# 定义降水日与极端降水阈值 pre_days pre[pre 0.1] # 95百分位阈值仅基于降水日 threshold_95 pre_days.quantile(0.95) print(f极端降水阈值95百分位: {threshold_95:.1f} mm) # 统计每年极端降水日数 df[is_extreme] df[pre] threshold_95 annual_extreme_days df[is_extreme].resample(YE).sum() # 年降水总量 annual_pre df[pre].resample(YE).sum() print(annual_extreme_days.tail())用quantile(0.95)代替固定阈值可以让不同站的极端事件定义在统计上可比。需要留意的是计算百分位时只用降水日样本而不是全部天数否则0值会把阈值压得很低导致「极端事件」里混入普通降雨。输出每年的极端降水日数后可以看到它的年际波动通常比温度大得多——降水本来就是高变异变量一年的极端日数从3天跳到10天不罕见。4.3 连续干期与雨季集中度用滑动窗口量化干旱与汛期除了极端降水事件农业和水利上更关心的是无雨持续时间和降水集中在哪个季节。连续无雨日是抗旱决策的直接指标而降水集中度决定了水库调度策略。两者都可以用简单统计实现不需要复杂的干旱指数模型。# 计算每年最大连续无雨日数 pre_binary (df[pre] 0.1).astype(int) # 无雨日标记为1 max_dry_spell {} for year, group in pre_binary.groupby(pre_binary.index.year): # 按连续为1的块分组取最大块长度 spell group.groupby((group ! group.shift()).cumsum()).sum() max_dry_spell[year] spell.max() if len(spell) 0 else 0 # 计算雨季5-9月降水占全年比例 seasonal df.groupby(df.index.month)[pre].sum() rainy_season_ratio ( seasonal.loc[5:9].sum() / seasonal.sum() ) print(f汛期5-9月降水占比: {rainy_season_ratio:.1%})连续无雨日的计算用了shift加cumsum的分组技巧——当天的无雨状态与前一日不同时分组编号加1从而把连续的无雨日块切分出来。max_dry_spell的结果直接反映季节性干旱风险。雨季集中度则简单得多按月汇总后算5到9月占比华北很多站点这个比例超过70%意味着汛期之外的半年几乎无降水。5. 气象数据分析避坑五个高频翻车点与排查方法5.1 时区陷阱UTC与北京时间的8小时偏移现象把再分析资料的日降水和观测站数据放在同一张图里对比发现降水事件总是对不上甚至日期都差了一天。原因再分析资料ERA5、NCEP时间基准是UTC而国内观测站数据通常是北京时间两者相差8小时。日值数据如果按自然日聚合UTC的一天从北京时间早8点开始夜间到早晨的降水会被划分到错误的日期。解决在读取网格数据后先把时间索引转成目标时区再按日重采样。# 把UTC时间转为北京时间UTC8再按日聚合 df_utc df_era5.copy() df_utc.index df_utc.index.tz_localize(UTC).tz_convert(Asia/Shanghai) df_daily_utc df_utc.resample(D).mean()tz_localize先给无时区的时间戳标记为UTCtz_convert再转到北京时间。如果原始数据本身有时区信息直接tz_convert即可。这个转换必须在所有聚合操作之前做顺序反了会造成无法挽回的偏移。5.2 缺测值不一定是NaN32744、9999、-9999现象读进来的数据没有NaN但画趋势图时出现一个离谱的尖峰比如7月气温突然变成9999℃。原因气象数据文件早期用特定数值表示缺测不同机构约定不同常见的有32744、32766、9999、-9999这些值没有在CSV里标记为缺失被当成真实观测读入。解决读文件前先打开原始数据说明把缺测值通过na_values传入或者读进来后用条件筛选统一替换。# 检查有没有未识别的缺测值 for col in df.columns: for bad_val in [32744, 32766, 9999, -9999]: n_bad (df[col] bad_val).sum() if n_bad 0: print(f列 {col} 中发现缺测值 {bad_val}数量 {n_bad}) # 统一替换为NaN df df.replace([32744, 32766, 9999, -9999], np.nan)这段排查代码建议在清洗流程里固定跑一次特别是拿到新数据源时。气象站数据的说明文档经常藏在下载包的readme里里面会写明缺测值编码和单位信息花两分钟读一遍能省一整天排错时间。5.3 闰年与2月29日重采样的隐蔽多一天现象按年对比各月均值时每年2月的平均值出现小幅但系统的偏差某几年数值明显不同。原因2月有28天和29天两种长度如果按自然月重采样日值数量不同导致月均值权重变化更隐蔽的是dayofyear在闰年3月之后会比平年多一天导致逐日气候态对齐错位。解决月重采样用resample(ME)可以避免多数问题计算逐日气候态时对闰年的2月29日做单独处理。# 剔除2月29日保证逐日气候态对齐 df_clean df[ ~((df.index.month 2) (df.index.day 29)) ] # 重新按dayofyear算气候态 daily_clim df_clean.groupby(df_clean.index.dayofyear).mean()去掉2月29日只损失一个样本对月均值和气候态的影响微乎其微但能消除每年多一天带来的对齐偏差。如果做的是逐日气候态产品发布这个步骤特别关键否则3月1日之后每一天的气候态基准都会错位一天。5.4 网格数据取点最近邻不是「取整」现象用sel(lat39.9, lon116.4)从NetCDF里取北京的气温和站点实测对比相关系数很高但绝对偏差一直存在而且偏差在不同季节大小不一。原因ERA5的经纬度网格是0.25°间隔格点坐标可能是39.75、40.0这种值直接传入的39.9并不落在格点上如果手写取整逻辑round到最近的0.25可能选到距离目标几十公里外的格点山地和沿海站点偏差尤其明显。解决用sel(..., methodnearest)它内部按距离选择最近格点更严格的场景先算经纬度距离矩阵再选最近邻。# 正确方式最近邻选取 ts_point ds[t2m].sel( lat39.9, lon116.4, methodnearest ) # 严谨方式先确认实际选中的格点坐标 print(f选中格点: lat{ts_point.lat.values:.2f}, lon{ts_point.lon.values:.2f})打印选中格点的坐标是个好习惯——它告诉你实际用的数据来自哪个网格点而不是想当然认为就是目标位置。地形复杂区域如果最近格点与实际站点海拔差超过500米气温偏差会达到3℃以上这时候应该改用双线性插值而不是最近邻。5.5 单位陷阱气温0.1℃、降水0.1mm现象画出的温度曲线在30到40℃之间剧烈震荡年降水量动辄上万毫米明显超出常识。原因国内很多气象数据为节省存储空间采用缩小10倍的整数编码——气温原始值代表0.1℃降水代表0.1mm风速代表0.1m/s。没有按说明除以10所有分析结果都会系统性放大10倍。解决在数据加载层统一做单位换算不要在使用时临时处理。# 加载时统一单位换算 unit_map {tem: 0.1, pre: 0.1, wind: 0.1} for col, factor in unit_map.items(): if col in df.columns: df[col] df[col] * factor这个坑几乎每个用过中国气象站点数据的人都踩过。我现在的习惯是写一个标准加载函数数据一进内存就完成单位换算和缺失值处理任何下游分析都不会再碰到原始值。另一个值得做的检查是画图后先看一眼纵轴范围——如果温度范围在-30到45℃之外第一反应不是分析问题而是单位或缺测处理出了问题。6. 进阶让气象分析结论更稳的三个实操技巧6.1 滑动平均去噪窗口宽度怎么定逐日温度曲线噪声大逐日降水曲线噪声更大。滑动平均是简单有效的去噪方式但窗口宽度要根据你想保留的信号尺度来选择7天窗口能去掉天气尺度波动保留一周以上的天气过程30天窗口适合看月际背景如果想看年际信号直接用年均值比任何滑动平均都干净。窗口宽度没有统一标准我在报告里通常同时画原始序列和两种窗口的滑动平均线让读者自己判断。df[tem_smooth_7d] df[tem].rolling(window7, centerTrue).mean() df[tem_smooth_30d] df[tem].rolling(window30, centerTrue).mean()centerTrue让滑动平均窗口以当前天为中心而不是只取过去这样平滑后的序列不会整体滞后。注意到窗口开头和结尾会有NaN这是边界效应绘图时可选择从有值的位置开始。6.2 相关分析的样本量陷阱气象时间序列的自相关统计两个气象变量比如温度和湿度的Pearson相关系数时如果直接把逐日数据全丢进去n值看起来很大十年有3650个样本但气象时间序列存在强自相关——今天的温度和昨天几乎一样有效样本量远小于实际样本量。直接按n3650查显著性表会把不相关的变量判成显著相关。修正办法是用有效样本量计算自由度。def effective_n(series1, series2): 计算考虑自相关后的有效样本量 n len(series1) r1 series1.autocorr(lag1) r2 series2.autocorr(lag1) if abs(r1) 1 or abs(r2) 1: return n n_eff n * (1 - r1 * r2) / (1 r1 * r2) return int(n_eff)使用autocorr(lag1)得到滞后1阶自相关系数代入修正公式。常见替代方案是先把数据聚合成月距平再做相关月距平序列的自相关大幅降低有效样本量的问题自然缓解。在报告里标注有效样本量是严谨性的加分项也是很多论文质检点。6.3 验证结论的两个办法分段时间对比与交叉数据任何趋势或相关结论在发布前都应该经过验证。第一个办法是分段稳定性检验把序列按时间切成前后两半分别计算趋势如果两半的趋势符号不同说明全序列的趋势可能是气候变率噪声而不是稳定趋势。第二个办法是交叉数据验证用再分析资料算同一个站点的趋势或极端事件频率与观测站结果对比符号和量级一致才说明结论不是单一数据源的系统误差造成的。这两个验证方法不需要额外代码复用前面章节的回归和统计函数换数据切片就行。我的习惯是把所有分析封装成以DataFrame为输入的函数验证时只需要传入不同时段或不同来源的数据。有一次我在一个站点上算出了显著的增温趋势分段检验发现前半段几乎无趋势、后半段急剧上升这个信息比单一斜率值重要得多——它说明了变化发生的时段而不是笼统一句「在变暖」。气象数据分析做久了会发现大部分翻车不是统计模型不够高级而是数据在进入模型之前就错了。我现在的项目里固定保留一个数据校验层单位、时区、缺测、闰年四道检查全部跑完才开始分析这套习惯的养成比任何算法都值钱。希望帮到你。本文还有配套的精品资源点击获取