简介面向数据分析课程与气象场景实训的资源包适合高校学生、数据竞赛新手快速上手爬虫与可视化分析。内容以中国天气网某城市天气数据为对象涵盖爬虫抓取、数据整理、雷达图与条形图生成、分析说明完整呈现从数据采集到结论输出的流程。压缩包共39个文件大小1.69MB包含20张PNG分析图表、15个XML文档组件、4个rels关联文件其中图表用于直接查看可视化结果XML与rels构成实验报告Word文档的结构化组件源代码可在Jupyter Notebook中直接运行复现。目前已有1853人浏览学习可用于课程作业、实验报告参考也能帮助学习者独立完成一次天气数据分析闭环并借鉴其中的数据清洗思路与可视化表达。1. 气象分析不是画曲线它先得把天气拆成能算的数据前一阵某业务平台要复盘“为什么今年园区夏季异常高温”拿到气象数据包的人第一反应是打开站点文件画一条月平均气温曲线。结果半天就被卡住时间戳既有本地时间又有世界时站点记录里埋着大量缺测值个别时次还算出 50℃ 的“气象异常”——这些异常信号全是脏数据的产物。气象分析数据分析说的正是这么一回事它不是看云图而是用数据工程手段把大气观测和模式输出整理成能进统计检验的干净面板再从中读出趋势和极端事件。这几年我做过不少类似项目站点观测、区域再分析格点场、雷达回波都碰过。我的体感是气象分析的技术栈本身不复杂真正吃时间的是数据形态切换、时空对齐、质量控制和口径定义这四件事通常占掉七成工作量。下面按我平时处理这类问题的顺序来讲怎么读入、怎么对齐、怎么算统计量、坑在哪最后给一个能直接上手的验证习惯。适合手里有气象数据要做分析报告的从业者也适合想转行做气象数据分析的新手。2. 数据读入与格式识别站点、格点、雷达三类气象数据怎么进内存气象数据不能用“一个 CSV 装天下”来理解。动手前先把数据形态认清比急着写 pandas 重要得多因为后续所有代码都建立在对结构的判断上。2.1 先分清数据形态站点观测、格点再分析、雷达与卫星反演数据形态典型载体基本结构分析用途常用读取工具站点观测CSV/文本/数据库时间戳 站号 若干要素定点实况、质量检验、真值对比pandas格点再分析NetCDF/GRIB数组time/lat/lon区域气候背景、格局分析xarray雷达基数据二进制/HDF5极坐标扫描体含反射率/径向速度强对流识别、短临降水估测pyart 或设备专用接口卫星反演HDF/NetCDF扫描条带常带云污染标识云宏观、降水估计、大尺度环流satpy 等专用工具这几类数据在落地分析里的角色不一样。再分析格点场是“模式与观测融合后的产物”空间完整、时间连续适合做气候态和趋势这类需要大样本的统计站点观测则零散但更接近真实适合做验证基准。雷达数据时间分辨率高但覆盖范围窄、质量受地物杂波干扰严重通常只在强对流个例分析里单独用很少直接进入长期趋势计算。我拿到一个数据包的第一动作是先确认它是哪一类再确认每类形态对应的元数据字段。比如 NetCDF 文件里的维度名可能是 lon/lat也可能是 longitude/latitude、x/y变量名可能是 t2m也可能是 tem/tmp/TT。不要凭文件名猜结构最好直接打开看坐标和变量。2.2 用 xarray 打开 NetCDF 格点场的最小闭环格点数据我最常碰到的就是 NetCDF。一个典型的业务场景拿到的区域再分析文件是逐小时或逐三小时输出大小从几百 MB 到几十 GB 不等。对这种文件第一步不是写转换脚本而是先打印坐标和变量名确认时间步长、经纬度范围和单位。import xarray as xr import pandas as pd ds xr.open_dataset(region_grid_sample.nc) # 打印维坐标和变量名确认时间步长、经纬度范围与单位 print(ds.coords) print(ds.data_vars) # 按经纬度切片以某研究区域为例经度 98~105纬度 35~40 # lat 倒序时slice 边界要大值在前 region ds.sel(lonslice(98, 105), latslice(40, 35)) # 把区域场转成长表方便和站点数据合并做后续分析 df_grid region[t2m].to_dataframe().reset_index() # 再分析和模式产品通常用开尔文存温度业务习惯用摄氏度 df_grid[t2m_c] df_grid[t2m] - 273.15 print(df_grid.head())逻辑说明to_dataframe会把多维数组展平成“多维索引 数值”的行形式等价于一张时间—经纬度—温度的长表这是和站点数据合并的基础。但要注意坐标范围较大时“网格点数 × 时次数”会迅速膨胀直接全部读进内存可能不够长序列建议先sel切片再load()。几个参数值得单独说。lonslice(98, 105)是闭区间左右边界都会被包含latslice(40, 35)下界比上界大是因为常见格点文件的纬度由南向北保存反着写才能得到目标区域正着写会直接得到空数组。- 273.15是开尔文与摄氏零度的差但前提是确认源单位确实是 K有的区域产品直接存摄氏度反复“转换”是这类代码最常见的翻车点。提示open_dataset默认是惰性加载只有真正取值时才读磁盘。大文件不要立刻load()完整份先做范围筛选。2.3 站点 CSV 的读入习惯与雷达数据的“非表格式”处理站点观测通常是多站共用一张表。拿到这种表我第一件事是确认两样东西站号是不是唯一标识时间列有没有时区信息。很多历史资料里同一站号在不同阶段代表不同位置比如站点迁移后站号保留、数据继续累积直接按站号分组算年均值相当于把两段不同地理坐标的记录混成一个序列。处理思路是给站号加版本后缀再判断是否拼接。站点 CSV 的读入本身不复杂但有几个默认习惯要改obs pd.read_csv(obs_sample.csv, parse_dates[time]) obs obs.replace(-9999, pd.NA) # 先把缺测标记换掉防止混进统计 obs[year] obs[time].dt.year obs[station_id] obs[station_id].astype(str) # 站号可能以 0 开头别读成数字这段代码看起来普通实际价值在最后两行把站点编号保留为字符串是为了防止“一号站点被强转成整数后丢前导零”这是多源站点合并时非常隐蔽的坑。雷达基数据则完全不同。它通常以极坐标形式存储无法直接用 pandas 读取我的习惯是先调用配套接口把极坐标转成经纬度网格并落成标准 NetCDF再回到 xarray 流程。这样做的理由是分析代码只依赖一个统一的数据接口后期换数据源时不用重写逻辑。到这里三类数据链路已经理清站点靠 pandas格点靠 xarray雷达转一次格式后并入 xarray。接下来要解决的是时间与空间的坐标统一问题这是整个分析链条里代价最高的部分。3. 时空对齐与插值把乱序观测整理成一张能分析的标准表无论数据来自哪里最终想要的分析表形状是一致的每一行是一个站点在某时刻的某个物理量时间按固定分辨率对齐。原始资料往往不满足这个形状所以要先处理时间和空间两个对齐轴。3.1 时间清洗时区换算与整点规整气象资料的通用时间标准是 UTC而不少站点业务系统存的是当地标准时间。直接pd.to_datetime不会自动识别时区只会把时间当作“无时区对象”处理和 UTC 数据合并时就会整体偏移数小时。# 站点文件给的是业务当地时本例为东八区模式资料通常是 UTC obs[time_utc] ( obs[time] .dt.tz_localize(Etc/GMT-8) # 先声明这是当地时 .dt.tz_convert(UTC) # 再转成世界时 .dt.tz_localize(None) # 去掉时区标记方便与 pandas 普通时间列对齐 ) # 把分钟规整到整点避免 10:31、10:52 这类记录造成多余时层 obs[time_hour] obs[time_utc].dt.floor(h) # 与格点资料按最近时次拼接同一站内匹配超过 30 分钟不匹配则丢弃 merged pd.merge_asof( obs.sort_values(time_hour), df_grid.sort_values(time), left_ontime_hour, right_ontime, bystation_id, tolerancepd.Timedelta(30min), )这里有两个函数很容易搞混。tz_localize只是给无时区的数组贴上一个时区标记并不改变时间数字本身tz_convert才真正调整钟表时刻。最后的tz_localize(None)把带时区时间重新转回朴素时间对象消除与 pandas 普通时间列合并时的类型冲突。merge_asof适合处理“观测是分钟级、模式输出是整点级”的拼接它不强制观测对齐到整点而是保留真实时间并按最近时次取值。参数bystation_id保证合并只在同一站内进行tolerance30min规定偏差超过半小时的观测直接丢弃避免把前一天的几乎零点数据错配到后一天。这个容忍度需要按业务频率调整逐小时观测的资料可以放宽分钟级资料可以收紧。3.2 空间对齐站点与格点互相靠拢的两种常用做法空间对齐要回答两个问题把格点值解释到站点位置或者把站点观测外推到格点上。第一种常用双线性插值第二种常用反距离加权IDW。from scipy.spatial import cKDTree import numpy as np def idw_at(lon_t, lat_t, pts, k4, p2): pts: n×3 数组列为 lon, lat, value 返回目标点的反距离加权插值结果 coords pts[:, :2] tree cKDTree(coords) dist, idx tree.query([lon_t, lat_t], kk) # 目标点正好落在站点上时直接取该点值避免权重除零 if np.isscalar(dist): return pts[idx, 2] w 1.0 / (dist ** p 1e-6) return np.dot(w, pts[idx, 2]) / w.sum()这个函数有三个参数要理解。邻居数 k 决定插值场的光滑程度k1 会得到一块块“马赛克”k 太大又会让远处站点穿透局部特征区域气象资料我一般从 k4 开始调。幂次 p 控制距离衰减速度p 越大结果越接近最近邻p2 是常用默认值。用 xarray 自带插值直接从格点取站点值是更贴近业务的路径# 对每个站点坐标在格点场上做双线性插值 site_values region[t2m].interp( lonxr.DataArray(site_lons, dimssite), latxr.DataArray(site_lats, dimssite), methodlinear, )methodlinear是双线性methodnearest是最近邻。对气温这类连续变量我默认用 linear对降水这种大量 0 值的物理量反而用 nearest因为线性插值会把两个晴雨不同的站中间“插出小雨”这在业务上很难接受。3.3 单位换算与质量控制的“物理上限”多源资料最常见的混乱是单位不统一温度有开尔文和摄氏度相对湿度有 0~1 和 0~100 两套写法降水有 mm/h 和 mm/day。在做任何统计之前把单位统一成目标量纲并给每个要素设物理上限。def qc_temperature(series): # 典型地面气温物理范围超出区间直接设缺测 return series.where((series -70.0) (series 60.0)) def qc_rh(series): # 相对湿度有的存 0~1有的存 0~100先统一到 0~100 再过滤 if series.max() 1.0: series series * 100.0 return series.where((series 0.0) (series 100.0))这两段没有玄学每条边界都要有依据气温低于 -70℃ 基本只可能出现在极地冰盖高空普通业务站点不会出现相对湿度超过 100% 通常是传感器结冰或雨滴击中感应端导致这类记录的处理方式不是截断到 100而是整条置缺测因为超出量程本身就说明设备状态不可信。还有一个容易漏掉的检查只要分析用到日均值或月均值就必须返回去看当天原始记录到底有多少条。某站一天只有 2 条观测和正常 24 条观测算出来的日均值方差完全不同。我习惯在汇总前额外生成一列daily_count低于阈值的日期不参与后续统计。4. 统计与极端事件气候态基准、趋势检验和极端天气怎么界定数据整理干净之后气象分析的核心问题就从“怎么读”变成“怎么下结论”。这一章把最常用的三件事讲清楚气候态基准、趋势检验、极端事件判别。4.1 30年基准与距平把“今年热不热”变成标准化数字“今年夏天比常年热”的“常年”不是任意历史平均而是固定的 30 年气候基准通常取最近三个完整十年。先用基准期算日序气候态再逐日做差得到距平这是气象分析里最基础也最好用的一招。# 只取基准期并剔除 2 月 29 日避免闰年单独拉高某天的样本量 base df[df[time].dt.year.between(1991, 2020)] base base[~((base[time].dt.month 2) (base[time].dt.day 29))] # 多年平均的日序气候态例如 07-15 这天的多年平均气温 clim_doy base.groupby(base[time].dt.strftime(%m-%d))[t2m_c].mean() # 逐年逐日距平 当日实测值 - 当天气候态值 df[anom] df[t2m_c] - df[time].dt.strftime(%m-%d).map(clim_doy) # 看某年盛夏的距平分布中位数和分位数能说明问题 year_view df[df[time].dt.year 2024] print(year_view[year_view[time].dt.month.isin([6, 7, 8])][anom].describe())strftime(%m-%d)把日期转成“07-15”这种字符串再作为索引 map 回气候态比手动把月份和日期拆开再合并少出错。还有个细节基准期里剔除 2 月 29 日是因为闰日会让该日序的样本量不一致气候态均值会被少数年份拉偏。如果想进一步比较不同区域可以把距平标准化即除以基准期内该日序的标准差。标准化距平的物理意义是“偏离常年几个标准差”严寒地区的两度距平可能是强信号热带地区的两度距平只是日常噪音标准化能让两个区域放进同一张图里比较。4.2 趋势检验线性斜率为什么不能只看一条直线“这个区域在变暖吗”是高频需求。直接对历年平均气温做线性回归很容易被个别极端年份拉偏。气象里更常用的是 Mann-Kendall 检验判断趋势是否存在Sen 斜率估计趋势强度两者对异常值都不敏感。from scipy.stats import kendalltau, theilslopes annual df.groupby(df[time].dt.year)[t2m_c].mean().dropna() # 趋势显著性检验p 0.05 可认为趋势不是随机波动 tau, p_val kendalltau(annual.index.values, annual.values) # Sen 斜率估计单位是每时间步的变化量 slope, low, high theilslopes(annual.values, annual.index) print(fp{p_val:.3f}, trend{slope*10:.2f}°C/10年)kendalltau返回的tau是排序相关性大小本身不直接代表升温幅度真正写进结论的是 p 值。p 大于 0.05 时即使散点图看着在升也应写成“趋势未通过显著性检验”而不是“有升温趋势”。业务复盘里因为“看着像”写错结论的案例不少见。theilslopes第一组输出是斜率第二组是置信区间。代码里slope*10是固定换算因为annual.index的单位是“年”一个时间步就是一年乘 10 就是“十年升温量”。这样报数可以避免口头换算错位。这个检验假设序列没有强自相关。年序列的自相关通常不严重但如果换成月序列必须先去掉季节循环否则冬夏的周期波动会污染秩检验。4.3 极端事件判别百分位阈值与连续天数业务上界定极端高温不只看“哪天超过 37℃”更常用的是“超过某百分位阈值并持续多少天”。百分位法天然消除地区差异北方 35℃ 是极端南方 35℃ 可能只是普通夏日用绝对阈值会漏掉局部极端事件。# 以暖季样本算第 95 百分位作为高温阈值 warm df[df[time].dt.month.isin([5, 6, 7, 8, 9])] thr warm[t2m_c].quantile(0.95) # 找出超过阈值的日序列并给每天标记连续性分组 hot (df[t2m_c] thr).astype(int) df[run_id] (hot.diff() ! 0).cumsum() # 统计每次极端事件的起止时间和持续天数 events (df[hot 1] .groupby(run_id) .agg(start(time, min), end(time, max), dur(time, size)) .reset_index()) events events[events[dur] 3] # 连续 3 天及以上才记为一次热浪 print(events)生成run_id那两行是核心技巧对 0/1 序列做差分差分不等于 0 说明状态翻转累积和递增于是每个连续段得到唯一编号。这种方法在处理连阴雨、持续雾天等场景同样通用。有两个反直觉点要提前处理。一是跨年时事件会被切成两段12 月 31 日开始的极端事件到次年 1 月会断掉需要按“自然年 run_id”共同分组再判断连续性。二是百分位阈值本身随样本变化温度尽量用暖季样本降水则用全年样本并叠加绝对阈值比如“日降水量 ≥ 50mm 且超过当地 99 分位”这种叠加定义在汇报时解释成本最低。5. 气象数据分析避坑时间边界、空间尺度与数据源的 5 个坑这一章集中讲我实际项目中踩过的坑。每条按现象、原因、解决来写遇到类似情况可以直接对照排查。5.1 跨年冬季被截断或重复计算现象分析 12 月到次年 2 月的冬季平均气温输出结果明明有这几个月却少了一段平均值整体偏暖。 原因切片通常写成“月份 12 或月份 2”但没有同时按年份分组同一个 12 月被前后两个“冬季”重复计入。 解决先构造气候年字段把 12 月归属到下一年df[season_year] df[time].dt.year.where( df[time].dt.month 3, df[time].dt.year - 1 )然后按season_year分组。不同业务对“冬季”的起止定义有差别换项目时先问清楚再算。5.2 2月29日让气候态均值出现细小的系统偏差现象多年平均的日序曲线在 2 月末出现一个小峰单看不明显但和对比资料放一起就很扎眼。 原因闰日四年一遇样本量只有平年日期的四分之一均值容易被少数极端年份拉动。 解决基准期先剔除 2 月 29 日。还可以把平年 2 月 28 日复制一份作为闰日的气候替代值保证当年逐日都有气候态可查算距平时不缺口。5.3 缺测标记被当成真实观测极值统计全线失真现象统计某区域极端高温时多个站在同一时次出现完全相同的高值而且该数值在记录里反复出现。 原因很多观测文件用全局特殊值标记缺测比如 -9999、9999按数值解析时把这些标记当成了真实气温。 解决读入后先替换再过滤df df.replace(-9999, pd.NA) df[t2m_c] qc_temperature(df[t2m_c])如果文件里带质量控制码列所有标明“可疑”的记录一律按缺测处理。假数据是抬高极端阈值的元凶只靠物理范围过滤不够。5.4 沿海插值出现海陆数据打架现象把稀疏站点用 IDW 插值到沿海格点后陆地一侧出现明显温度斑块和海面再分析数据对比像两个世界。 原因IDW 不区分陆地和海洋海岸线附近的站点把海洋网格同化成陆地温度场模式的海面温度又反过来污染了陆侧插值。 解决先对格点场做海陆掩膜插值后把落在海上的网格置缺测或把站点有效插值半径限制在 30~50km。山区站点的影响范围往往比平原小得多限制半径虽然是经验值但对控制跨界污染很有效。5.5 时间戳归整方向不一致滑动平均整体偏移现象用小时数据算 3 小时滑动平均结果和基准产品对比始终差一个时次峰值整体左移。 原因有人用floor(3h)有人用round(3h)还有人把“02/05/08 点”直接当成整点时间但资料实际记录的是时段累计值。错位之下平均值自然对不上。 解决先明确时间标签是“时段开始”还是“时段结束”。如果是结束时间要做一次归位df[time_adj] df[time] - pd.Timedelta(hours3) # 存的是结束时间改为起始时刻这个坑最常出现在降水、辐射累计量和雷达估测产品里气温这类瞬时值基本没有此问题。多源合并项目务必把“瞬时量”和“累计量”分开列。6. 用泰勒图给分析结果做“体检”一个能直接复现的验证习惯整理出距平、趋势和极端事件表后我通常不直接写报告而是先做一轮结果体检。体检对象是插值或重采样后的格点序列对照基准是原始站点观测用三个统计量来打分相关系数、归一化标准差、中心化均方根误差这三个量组成泰勒图的骨架。为什么非做不可气象数据经过插值、再分析、极端阈值计算后人的直觉经常失灵——曲线看着挺像标准差可能已经缩水 30%极端事件统计量随之失真。我只信数字。import numpy as np def taylor_stats(obs, est): obs np.asarray(obs, dtypefloat) est np.asarray(est, dtypefloat) # 扣除各自均值后比较波动形态避免平均偏差掩盖问题 obs_c obs - obs.mean() est_c est - est.mean() cc np.corrcoef(obs_c, est_c)[0, 1] rmse_c np.sqrt(np.mean((obs_c - est_c) ** 2)) ratio np.std(est_c) / np.std(obs_c) return cc, rmse_c, ratio三个指标各有分工相关系数反映时序波动是否一致标准差比反映插值把振幅磨平还是夸大中心化 RMSE 把前两者综合成一个损失指标。对气温来说插值结果与站点观测的相关系数在 0.95 以上才适合业务使用标准差比落在 0.9~1.1 之间可接受低于 0.85 说明插值过度平滑。多组对比用表格列示会更直观对比项相关系数标准差比中心化RMSE结论格点再分析 vs 站点0.970.941.1℃可用IDW插值 vs 站点0.910.781.9℃振幅被压慎用雷达雨量 vs 雨量计0.821.253.5mm系统偏强需订正做完整套体检后我还有个习惯性动作每年底把当年的观测均值、距平、极端事件计数和过去五年同一口径的统计结果列在一起对比。如果某站“极端高温天数”从 5 天突然跳到 30 天先不要写“气候突变”而是回头查站点是否搬迁、仪器是否更换、质量控制码是否变化。这套回查逻辑比任何高级统计都省钱也能避免在报告里闹笑话。气象分析这门手艺的大部分价值不在算法先进而在口径稳定、流程可复现、结果可验证。我自己的习惯是先验后算、再出结论一套能解释三个统计量的固定模板让我少跑偏很多次。希望这份从读入、对齐、统计、避坑到验证的完整链条能在你处理气象数据时真正帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/21 22:13:07 · 阅读 1.2k 气象数据分析实战：从数据清洗到极端事件判别的完整流程
简介面向数据分析课程与气象场景实训的资源包适合高校学生、数据竞赛新手快速上手爬虫与可视化分析。内容以中国天气网某城市天气数据为对象涵盖爬虫抓取、数据整理、雷达图与条形图生成、分析说明完整呈现从数据采集到结论输出的流程。压缩包共39个文件大小1.69MB包含20张PNG分析图表、15个XML文档组件、4个rels关联文件其中图表用于直接查看可视化结果XML与rels构成实验报告Word文档的结构化组件源代码可在Jupyter Notebook中直接运行复现。目前已有1853人浏览学习可用于课程作业、实验报告参考也能帮助学习者独立完成一次天气数据分析闭环并借鉴其中的数据清洗思路与可视化表达。1. 气象分析不是画曲线它先得把天气拆成能算的数据前一阵某业务平台要复盘“为什么今年园区夏季异常高温”拿到气象数据包的人第一反应是打开站点文件画一条月平均气温曲线。结果半天就被卡住时间戳既有本地时间又有世界时站点记录里埋着大量缺测值个别时次还算出 50℃ 的“气象异常”——这些异常信号全是脏数据的产物。气象分析数据分析说的正是这么一回事它不是看云图而是用数据工程手段把大气观测和模式输出整理成能进统计检验的干净面板再从中读出趋势和极端事件。这几年我做过不少类似项目站点观测、区域再分析格点场、雷达回波都碰过。我的体感是气象分析的技术栈本身不复杂真正吃时间的是数据形态切换、时空对齐、质量控制和口径定义这四件事通常占掉七成工作量。下面按我平时处理这类问题的顺序来讲怎么读入、怎么对齐、怎么算统计量、坑在哪最后给一个能直接上手的验证习惯。适合手里有气象数据要做分析报告的从业者也适合想转行做气象数据分析的新手。2. 数据读入与格式识别站点、格点、雷达三类气象数据怎么进内存气象数据不能用“一个 CSV 装天下”来理解。动手前先把数据形态认清比急着写 pandas 重要得多因为后续所有代码都建立在对结构的判断上。2.1 先分清数据形态站点观测、格点再分析、雷达与卫星反演数据形态典型载体基本结构分析用途常用读取工具站点观测CSV/文本/数据库时间戳 站号 若干要素定点实况、质量检验、真值对比pandas格点再分析NetCDF/GRIB数组time/lat/lon区域气候背景、格局分析xarray雷达基数据二进制/HDF5极坐标扫描体含反射率/径向速度强对流识别、短临降水估测pyart 或设备专用接口卫星反演HDF/NetCDF扫描条带常带云污染标识云宏观、降水估计、大尺度环流satpy 等专用工具这几类数据在落地分析里的角色不一样。再分析格点场是“模式与观测融合后的产物”空间完整、时间连续适合做气候态和趋势这类需要大样本的统计站点观测则零散但更接近真实适合做验证基准。雷达数据时间分辨率高但覆盖范围窄、质量受地物杂波干扰严重通常只在强对流个例分析里单独用很少直接进入长期趋势计算。我拿到一个数据包的第一动作是先确认它是哪一类再确认每类形态对应的元数据字段。比如 NetCDF 文件里的维度名可能是 lon/lat也可能是 longitude/latitude、x/y变量名可能是 t2m也可能是 tem/tmp/TT。不要凭文件名猜结构最好直接打开看坐标和变量。2.2 用 xarray 打开 NetCDF 格点场的最小闭环格点数据我最常碰到的就是 NetCDF。一个典型的业务场景拿到的区域再分析文件是逐小时或逐三小时输出大小从几百 MB 到几十 GB 不等。对这种文件第一步不是写转换脚本而是先打印坐标和变量名确认时间步长、经纬度范围和单位。import xarray as xr import pandas as pd ds xr.open_dataset(region_grid_sample.nc) # 打印维坐标和变量名确认时间步长、经纬度范围与单位 print(ds.coords) print(ds.data_vars) # 按经纬度切片以某研究区域为例经度 98~105纬度 35~40 # lat 倒序时slice 边界要大值在前 region ds.sel(lonslice(98, 105), latslice(40, 35)) # 把区域场转成长表方便和站点数据合并做后续分析 df_grid region[t2m].to_dataframe().reset_index() # 再分析和模式产品通常用开尔文存温度业务习惯用摄氏度 df_grid[t2m_c] df_grid[t2m] - 273.15 print(df_grid.head())逻辑说明to_dataframe会把多维数组展平成“多维索引 数值”的行形式等价于一张时间—经纬度—温度的长表这是和站点数据合并的基础。但要注意坐标范围较大时“网格点数 × 时次数”会迅速膨胀直接全部读进内存可能不够长序列建议先sel切片再load()。几个参数值得单独说。lonslice(98, 105)是闭区间左右边界都会被包含latslice(40, 35)下界比上界大是因为常见格点文件的纬度由南向北保存反着写才能得到目标区域正着写会直接得到空数组。- 273.15是开尔文与摄氏零度的差但前提是确认源单位确实是 K有的区域产品直接存摄氏度反复“转换”是这类代码最常见的翻车点。提示open_dataset默认是惰性加载只有真正取值时才读磁盘。大文件不要立刻load()完整份先做范围筛选。2.3 站点 CSV 的读入习惯与雷达数据的“非表格式”处理站点观测通常是多站共用一张表。拿到这种表我第一件事是确认两样东西站号是不是唯一标识时间列有没有时区信息。很多历史资料里同一站号在不同阶段代表不同位置比如站点迁移后站号保留、数据继续累积直接按站号分组算年均值相当于把两段不同地理坐标的记录混成一个序列。处理思路是给站号加版本后缀再判断是否拼接。站点 CSV 的读入本身不复杂但有几个默认习惯要改obs pd.read_csv(obs_sample.csv, parse_dates[time]) obs obs.replace(-9999, pd.NA) # 先把缺测标记换掉防止混进统计 obs[year] obs[time].dt.year obs[station_id] obs[station_id].astype(str) # 站号可能以 0 开头别读成数字这段代码看起来普通实际价值在最后两行把站点编号保留为字符串是为了防止“一号站点被强转成整数后丢前导零”这是多源站点合并时非常隐蔽的坑。雷达基数据则完全不同。它通常以极坐标形式存储无法直接用 pandas 读取我的习惯是先调用配套接口把极坐标转成经纬度网格并落成标准 NetCDF再回到 xarray 流程。这样做的理由是分析代码只依赖一个统一的数据接口后期换数据源时不用重写逻辑。到这里三类数据链路已经理清站点靠 pandas格点靠 xarray雷达转一次格式后并入 xarray。接下来要解决的是时间与空间的坐标统一问题这是整个分析链条里代价最高的部分。3. 时空对齐与插值把乱序观测整理成一张能分析的标准表无论数据来自哪里最终想要的分析表形状是一致的每一行是一个站点在某时刻的某个物理量时间按固定分辨率对齐。原始资料往往不满足这个形状所以要先处理时间和空间两个对齐轴。3.1 时间清洗时区换算与整点规整气象资料的通用时间标准是 UTC而不少站点业务系统存的是当地标准时间。直接pd.to_datetime不会自动识别时区只会把时间当作“无时区对象”处理和 UTC 数据合并时就会整体偏移数小时。# 站点文件给的是业务当地时本例为东八区模式资料通常是 UTC obs[time_utc] ( obs[time] .dt.tz_localize(Etc/GMT-8) # 先声明这是当地时 .dt.tz_convert(UTC) # 再转成世界时 .dt.tz_localize(None) # 去掉时区标记方便与 pandas 普通时间列对齐 ) # 把分钟规整到整点避免 10:31、10:52 这类记录造成多余时层 obs[time_hour] obs[time_utc].dt.floor(h) # 与格点资料按最近时次拼接同一站内匹配超过 30 分钟不匹配则丢弃 merged pd.merge_asof( obs.sort_values(time_hour), df_grid.sort_values(time), left_ontime_hour, right_ontime, bystation_id, tolerancepd.Timedelta(30min), )这里有两个函数很容易搞混。tz_localize只是给无时区的数组贴上一个时区标记并不改变时间数字本身tz_convert才真正调整钟表时刻。最后的tz_localize(None)把带时区时间重新转回朴素时间对象消除与 pandas 普通时间列合并时的类型冲突。merge_asof适合处理“观测是分钟级、模式输出是整点级”的拼接它不强制观测对齐到整点而是保留真实时间并按最近时次取值。参数bystation_id保证合并只在同一站内进行tolerance30min规定偏差超过半小时的观测直接丢弃避免把前一天的几乎零点数据错配到后一天。这个容忍度需要按业务频率调整逐小时观测的资料可以放宽分钟级资料可以收紧。3.2 空间对齐站点与格点互相靠拢的两种常用做法空间对齐要回答两个问题把格点值解释到站点位置或者把站点观测外推到格点上。第一种常用双线性插值第二种常用反距离加权IDW。from scipy.spatial import cKDTree import numpy as np def idw_at(lon_t, lat_t, pts, k4, p2): pts: n×3 数组列为 lon, lat, value 返回目标点的反距离加权插值结果 coords pts[:, :2] tree cKDTree(coords) dist, idx tree.query([lon_t, lat_t], kk) # 目标点正好落在站点上时直接取该点值避免权重除零 if np.isscalar(dist): return pts[idx, 2] w 1.0 / (dist ** p 1e-6) return np.dot(w, pts[idx, 2]) / w.sum()这个函数有三个参数要理解。邻居数 k 决定插值场的光滑程度k1 会得到一块块“马赛克”k 太大又会让远处站点穿透局部特征区域气象资料我一般从 k4 开始调。幂次 p 控制距离衰减速度p 越大结果越接近最近邻p2 是常用默认值。用 xarray 自带插值直接从格点取站点值是更贴近业务的路径# 对每个站点坐标在格点场上做双线性插值 site_values region[t2m].interp( lonxr.DataArray(site_lons, dimssite), latxr.DataArray(site_lats, dimssite), methodlinear, )methodlinear是双线性methodnearest是最近邻。对气温这类连续变量我默认用 linear对降水这种大量 0 值的物理量反而用 nearest因为线性插值会把两个晴雨不同的站中间“插出小雨”这在业务上很难接受。3.3 单位换算与质量控制的“物理上限”多源资料最常见的混乱是单位不统一温度有开尔文和摄氏度相对湿度有 0~1 和 0~100 两套写法降水有 mm/h 和 mm/day。在做任何统计之前把单位统一成目标量纲并给每个要素设物理上限。def qc_temperature(series): # 典型地面气温物理范围超出区间直接设缺测 return series.where((series -70.0) (series 60.0)) def qc_rh(series): # 相对湿度有的存 0~1有的存 0~100先统一到 0~100 再过滤 if series.max() 1.0: series series * 100.0 return series.where((series 0.0) (series 100.0))这两段没有玄学每条边界都要有依据气温低于 -70℃ 基本只可能出现在极地冰盖高空普通业务站点不会出现相对湿度超过 100% 通常是传感器结冰或雨滴击中感应端导致这类记录的处理方式不是截断到 100而是整条置缺测因为超出量程本身就说明设备状态不可信。还有一个容易漏掉的检查只要分析用到日均值或月均值就必须返回去看当天原始记录到底有多少条。某站一天只有 2 条观测和正常 24 条观测算出来的日均值方差完全不同。我习惯在汇总前额外生成一列daily_count低于阈值的日期不参与后续统计。4. 统计与极端事件气候态基准、趋势检验和极端天气怎么界定数据整理干净之后气象分析的核心问题就从“怎么读”变成“怎么下结论”。这一章把最常用的三件事讲清楚气候态基准、趋势检验、极端事件判别。4.1 30年基准与距平把“今年热不热”变成标准化数字“今年夏天比常年热”的“常年”不是任意历史平均而是固定的 30 年气候基准通常取最近三个完整十年。先用基准期算日序气候态再逐日做差得到距平这是气象分析里最基础也最好用的一招。# 只取基准期并剔除 2 月 29 日避免闰年单独拉高某天的样本量 base df[df[time].dt.year.between(1991, 2020)] base base[~((base[time].dt.month 2) (base[time].dt.day 29))] # 多年平均的日序气候态例如 07-15 这天的多年平均气温 clim_doy base.groupby(base[time].dt.strftime(%m-%d))[t2m_c].mean() # 逐年逐日距平 当日实测值 - 当天气候态值 df[anom] df[t2m_c] - df[time].dt.strftime(%m-%d).map(clim_doy) # 看某年盛夏的距平分布中位数和分位数能说明问题 year_view df[df[time].dt.year 2024] print(year_view[year_view[time].dt.month.isin([6, 7, 8])][anom].describe())strftime(%m-%d)把日期转成“07-15”这种字符串再作为索引 map 回气候态比手动把月份和日期拆开再合并少出错。还有个细节基准期里剔除 2 月 29 日是因为闰日会让该日序的样本量不一致气候态均值会被少数年份拉偏。如果想进一步比较不同区域可以把距平标准化即除以基准期内该日序的标准差。标准化距平的物理意义是“偏离常年几个标准差”严寒地区的两度距平可能是强信号热带地区的两度距平只是日常噪音标准化能让两个区域放进同一张图里比较。4.2 趋势检验线性斜率为什么不能只看一条直线“这个区域在变暖吗”是高频需求。直接对历年平均气温做线性回归很容易被个别极端年份拉偏。气象里更常用的是 Mann-Kendall 检验判断趋势是否存在Sen 斜率估计趋势强度两者对异常值都不敏感。from scipy.stats import kendalltau, theilslopes annual df.groupby(df[time].dt.year)[t2m_c].mean().dropna() # 趋势显著性检验p 0.05 可认为趋势不是随机波动 tau, p_val kendalltau(annual.index.values, annual.values) # Sen 斜率估计单位是每时间步的变化量 slope, low, high theilslopes(annual.values, annual.index) print(fp{p_val:.3f}, trend{slope*10:.2f}°C/10年)kendalltau返回的tau是排序相关性大小本身不直接代表升温幅度真正写进结论的是 p 值。p 大于 0.05 时即使散点图看着在升也应写成“趋势未通过显著性检验”而不是“有升温趋势”。业务复盘里因为“看着像”写错结论的案例不少见。theilslopes第一组输出是斜率第二组是置信区间。代码里slope*10是固定换算因为annual.index的单位是“年”一个时间步就是一年乘 10 就是“十年升温量”。这样报数可以避免口头换算错位。这个检验假设序列没有强自相关。年序列的自相关通常不严重但如果换成月序列必须先去掉季节循环否则冬夏的周期波动会污染秩检验。4.3 极端事件判别百分位阈值与连续天数业务上界定极端高温不只看“哪天超过 37℃”更常用的是“超过某百分位阈值并持续多少天”。百分位法天然消除地区差异北方 35℃ 是极端南方 35℃ 可能只是普通夏日用绝对阈值会漏掉局部极端事件。# 以暖季样本算第 95 百分位作为高温阈值 warm df[df[time].dt.month.isin([5, 6, 7, 8, 9])] thr warm[t2m_c].quantile(0.95) # 找出超过阈值的日序列并给每天标记连续性分组 hot (df[t2m_c] thr).astype(int) df[run_id] (hot.diff() ! 0).cumsum() # 统计每次极端事件的起止时间和持续天数 events (df[hot 1] .groupby(run_id) .agg(start(time, min), end(time, max), dur(time, size)) .reset_index()) events events[events[dur] 3] # 连续 3 天及以上才记为一次热浪 print(events)生成run_id那两行是核心技巧对 0/1 序列做差分差分不等于 0 说明状态翻转累积和递增于是每个连续段得到唯一编号。这种方法在处理连阴雨、持续雾天等场景同样通用。有两个反直觉点要提前处理。一是跨年时事件会被切成两段12 月 31 日开始的极端事件到次年 1 月会断掉需要按“自然年 run_id”共同分组再判断连续性。二是百分位阈值本身随样本变化温度尽量用暖季样本降水则用全年样本并叠加绝对阈值比如“日降水量 ≥ 50mm 且超过当地 99 分位”这种叠加定义在汇报时解释成本最低。5. 气象数据分析避坑时间边界、空间尺度与数据源的 5 个坑这一章集中讲我实际项目中踩过的坑。每条按现象、原因、解决来写遇到类似情况可以直接对照排查。5.1 跨年冬季被截断或重复计算现象分析 12 月到次年 2 月的冬季平均气温输出结果明明有这几个月却少了一段平均值整体偏暖。 原因切片通常写成“月份 12 或月份 2”但没有同时按年份分组同一个 12 月被前后两个“冬季”重复计入。 解决先构造气候年字段把 12 月归属到下一年df[season_year] df[time].dt.year.where( df[time].dt.month 3, df[time].dt.year - 1 )然后按season_year分组。不同业务对“冬季”的起止定义有差别换项目时先问清楚再算。5.2 2月29日让气候态均值出现细小的系统偏差现象多年平均的日序曲线在 2 月末出现一个小峰单看不明显但和对比资料放一起就很扎眼。 原因闰日四年一遇样本量只有平年日期的四分之一均值容易被少数极端年份拉动。 解决基准期先剔除 2 月 29 日。还可以把平年 2 月 28 日复制一份作为闰日的气候替代值保证当年逐日都有气候态可查算距平时不缺口。5.3 缺测标记被当成真实观测极值统计全线失真现象统计某区域极端高温时多个站在同一时次出现完全相同的高值而且该数值在记录里反复出现。 原因很多观测文件用全局特殊值标记缺测比如 -9999、9999按数值解析时把这些标记当成了真实气温。 解决读入后先替换再过滤df df.replace(-9999, pd.NA) df[t2m_c] qc_temperature(df[t2m_c])如果文件里带质量控制码列所有标明“可疑”的记录一律按缺测处理。假数据是抬高极端阈值的元凶只靠物理范围过滤不够。5.4 沿海插值出现海陆数据打架现象把稀疏站点用 IDW 插值到沿海格点后陆地一侧出现明显温度斑块和海面再分析数据对比像两个世界。 原因IDW 不区分陆地和海洋海岸线附近的站点把海洋网格同化成陆地温度场模式的海面温度又反过来污染了陆侧插值。 解决先对格点场做海陆掩膜插值后把落在海上的网格置缺测或把站点有效插值半径限制在 30~50km。山区站点的影响范围往往比平原小得多限制半径虽然是经验值但对控制跨界污染很有效。5.5 时间戳归整方向不一致滑动平均整体偏移现象用小时数据算 3 小时滑动平均结果和基准产品对比始终差一个时次峰值整体左移。 原因有人用floor(3h)有人用round(3h)还有人把“02/05/08 点”直接当成整点时间但资料实际记录的是时段累计值。错位之下平均值自然对不上。 解决先明确时间标签是“时段开始”还是“时段结束”。如果是结束时间要做一次归位df[time_adj] df[time] - pd.Timedelta(hours3) # 存的是结束时间改为起始时刻这个坑最常出现在降水、辐射累计量和雷达估测产品里气温这类瞬时值基本没有此问题。多源合并项目务必把“瞬时量”和“累计量”分开列。6. 用泰勒图给分析结果做“体检”一个能直接复现的验证习惯整理出距平、趋势和极端事件表后我通常不直接写报告而是先做一轮结果体检。体检对象是插值或重采样后的格点序列对照基准是原始站点观测用三个统计量来打分相关系数、归一化标准差、中心化均方根误差这三个量组成泰勒图的骨架。为什么非做不可气象数据经过插值、再分析、极端阈值计算后人的直觉经常失灵——曲线看着挺像标准差可能已经缩水 30%极端事件统计量随之失真。我只信数字。import numpy as np def taylor_stats(obs, est): obs np.asarray(obs, dtypefloat) est np.asarray(est, dtypefloat) # 扣除各自均值后比较波动形态避免平均偏差掩盖问题 obs_c obs - obs.mean() est_c est - est.mean() cc np.corrcoef(obs_c, est_c)[0, 1] rmse_c np.sqrt(np.mean((obs_c - est_c) ** 2)) ratio np.std(est_c) / np.std(obs_c) return cc, rmse_c, ratio三个指标各有分工相关系数反映时序波动是否一致标准差比反映插值把振幅磨平还是夸大中心化 RMSE 把前两者综合成一个损失指标。对气温来说插值结果与站点观测的相关系数在 0.95 以上才适合业务使用标准差比落在 0.9~1.1 之间可接受低于 0.85 说明插值过度平滑。多组对比用表格列示会更直观对比项相关系数标准差比中心化RMSE结论格点再分析 vs 站点0.970.941.1℃可用IDW插值 vs 站点0.910.781.9℃振幅被压慎用雷达雨量 vs 雨量计0.821.253.5mm系统偏强需订正做完整套体检后我还有个习惯性动作每年底把当年的观测均值、距平、极端事件计数和过去五年同一口径的统计结果列在一起对比。如果某站“极端高温天数”从 5 天突然跳到 30 天先不要写“气候突变”而是回头查站点是否搬迁、仪器是否更换、质量控制码是否变化。这套回查逻辑比任何高级统计都省钱也能避免在报告里闹笑话。气象分析这门手艺的大部分价值不在算法先进而在口径稳定、流程可复现、结果可验证。我自己的习惯是先验后算、再出结论一套能解释三个统计量的固定模板让我少跑偏很多次。希望这份从读入、对齐、统计、避坑到验证的完整链条能在你处理气象数据时真正帮到你。本文还有配套的精品资源点击获取