简介面向数据分析课程实验与气象数据入门探索的一份完整资料包。作者通过爬虫采集中国天气网某城市的天气信息完成数据清洗、统计分析与可视化并配有整理好的实验报告适合正在学习爬虫、数据分析与图表呈现的读者参考。压缩包共三十九个文件以二十张结果图、十五个文档组件和四个关联关系文件为主整体仅约一点六九兆字节结构清晰便于解压后直接打开源码在交互式编程环境中运行验证复现所有图表结果。目前已有一千八百五十二人浏览学习说明其在同类实验资源中具备参考价值。读者可获得完整爬虫源代码、分析脚本、雷达图与条形图等可视化成果及实验报告既能对照实验流程理解数据处理思路也能修改城市参数开展新的分析练习。1. 气象分析在解决什么问题从站点观测到可决策的气象事实过去做气候服务最常被问的一句话是“明天到底下不下雨”后来在能源和物流项目里问题变成了“这个区域未来三天的风电出力会掉多少”“连续高温会不会触发限电预案”。这些问题看似不同落到底层都是同一件事把站点观测、数值预报、再分析格点这些气象数据经过清洗、统计、插值和可视化变成一张能直接支撑决策的结论表。这就是气象分析数据分析的核心工作。它听起来离普通数据分析很远实际上技术栈高度重合只是多了时空概念、气象编码和物理单位的约束。适合做生产调度、农业灾害评估、气象服务系统、能源预测的工程师。一个反直觉的结论是这个方向真正的成本不在算力而在数据口径——同一份降水数据统计时段不同结论可能完全相反。2. 气象数据从哪来、长什么样NetCDF、GRIB 和站点 CSV 的读法差别做气象分析的第一步不是写算法而是把数据打开、看清它是哪个口径。我常接触到的气象数据源有三类地面站点观测、数值预报产品、再分析格点数据。它们回答的问题完全不同读取方式也完全不同先把这个分清楚后面才不会拿着 CSV 的思路去读二进制格点。数据类别典型来源常见分辨率常见格式延迟适合做什么站点观测国家地面站、区域自动站、NOAA ISD逐小时/逐日CSV、文本、NetCDF小时级到天级实况统计、极端事件复盘、校准预报数值预报GFS、ECMWF、区域模式0.25° 到几公里GRIB2、NetCDF分钟到小时级未来 1-15 天逐时要素预报再分析ERA5、ERA-Interim0.25° 至 0.1°NetCDF、GRIB几天到几个月历史长序列研究、缺少站点时的背景场三类数据里站点观测是“事实”但空间不连续再分析是“用同化算法把观测填到网格上的结果”空间连续但不等于实测数值预报是“从当前状态外推的未来”延迟最低但误差随时间增长。实际项目里这三者经常混用比如用站点观测做检验、用再分析做长时间序列统计、用 GFS 做未来几天的预警输入。2.1 再分析、站点观测、数值预报三种口径决定分析能回答的问题选数据源时先问自己要回答什么。做历史气候趋势分析直接用站点观测最可靠但站点分布不均山区可能几百公里没有一个站这时候 ERA5 这类再分析能补上空间覆盖但要注意它是模式同化结果高海拔地区的温度可能和站点实测差出好几度不能拿格点值当实测值用。做短期预警必须用数值预报但 GFS 在北京时间凌晨 02:30 左右才能拿到 00 时次早上 8 点要出结论就得接受一定时效损失。读取方式上CSV、NetCDF、GRIB 三类格式几乎是三个世界。CSV 是表格思维pandas 一把梭NetCDF 和 GRIB 是自描述二进制格点格式坐标轴、单位、变量名都写在文件里需要用 xarray 这类库去读。很多人拿到 GRIB 文件直接 pd.read_csv读出来一堆乱码这不是代码问题是格式认知问题。我的习惯是拿到数据先花两分钟看文件头确认变量名、单位、网格范围和坐标顺序再用对应的库去打开。2.2 用 xarray 打开 ERA5 的 NetCDF区域与时间切片的写法ERA5 的单层变量 NetCDF 是最适合入门的气象格点数据变量名直观、带单位、坐标齐全。读它不需要安装额外解码库xarray 直接支持。我一般会用xr.open_dataset打开然后用sel按时间和区域切片避免把全量数据加载进内存。import xarray as xr # 打开 ERA5 单层 NetCDF比如 2m 温度 ds xr.open_dataset(era5_single_202407.nc) # 按时间范围切 2024 年 7 月 1-2 日中国区域按经纬度切 ds_sub ds.sel( timeslice(2024-07-01, 2024-07-02), latitudeslice(60, 15), # 北半球从北往南切注意值是降序 longitudeslice(70, 140) # 经度要确认是 0-360 还是 -180-180 ) print(ds_sub)这里有两个易错点。其一ERA5 的latitude坐标是降序排列slice(60, 15)才能把北纬 60 度到 15 度这段切出来写反了就得到空数组。其二longitude的切片范围取决于文件本身如果文件经度是 0-360slice(70, 140)没问题如果是 -180-180就得写成slice(70, 140)照样能切但想切东经 120 度附近时两个坐标系的表达式不一样。拿到文件先print(ds.longitude.values)看一眼这个动作能省半小时。2.3 用 cfgrib 打开 GFS 的 GRIB2依赖和常见报错GFS 的预报产品以 GRIB2 为主这个格式比 NetCDF 麻烦xarray 需要借助 cfgrib 引擎才能读。cfgrib 底层依赖 ECMWF 的 eccodes 解码库安装和报错是新手最容易卡住的地方。# 推荐用 conda 装会自动带齐二进制依赖 conda install -c conda-forge cfgrib eccodesimport xarray as xr # 读取 GFS 单个时次的 GRIB2 文件 # 文件名类似 gfs.t00z.pgrb2.0p25.f000 ds xr.open_dataset(gfs.t00z.pgrb2.0p25.f000, enginecfgrib) print(ds)如果打开报错90% 是eccodes没装对或者 GRIB 文件里含多个变量但气压层不一致。解决方式是用backend_kwargs{errors: ignore}跳过暂时不需要的变量或者先用grib_ls命令行工具看一眼文件里到底有哪些字段。还要注意 GFS 的经度坐标默认是 0-360做中国区域分析时要先转成 -180-180这一步我在第五节专门讲。2.4 站点观测 CSV 的读取站号要当字符串缺测值先转 NaN站点观测数据最常见是 CSV 或文本格式看起来普通但有两个坑经常让人翻车。一个是站号是字符串直接用 pandas 默认读会变成 int前导零全部丢失站号 01001 变成 1001后续关联站点表直接失败。另一个是缺测值不同数据源写法不一样常见有 -999、9999、32744如果不处理做平均、求和时这些假数字会混进结果。import pandas as pd # 站点观测 CSV站号、时间(UTC)、纬度、经度、2m温度、降水、风速 df pd.read_csv( obs_202407.csv, dtype{station_id: str}, # 站号保留前导零 parse_dates[time], # 时间列解析成 datetime na_values[-999, 9999, 32744] # 一次性把缺测值转成 NaN ) print(df.head()) print(df.dtypes)na_values这个参数值得多说一句。气象数据的缺测码并不统一有些老资料用 9999有些用 -999.9有些用 3274416 位整型最大值附近。把它们统一转成NaN之后后续所有统计运算都走 pandas 的缺失值逻辑省得每写一个函数都去判断一次缺测码。这个习惯我在所有气象项目里都会保留成本极低收益极大。3. 把原始观测变成干净的时间序列物理范围检查、缺测填补与重采样数据能打开之后真正耗时间的部分才开始。气象观测数据的脏不是普通的脏是那种“看起来挺干净一算就出妖”的脏。比如夏天某站气温出现 81 度比如降水累计值出来一个负数比如同一个站在同一时刻出现两条记录。这些问题靠肉眼扫不出来必须在管道里前置一套清洗逻辑。3.1 先做物理合理性检查再做插值降水和温度的处理完全不同气象数据的清洗顺序有讲究不是上来就fillna而是先做物理合理性检查再做时间填补。物理合理性检查就是把每个气象要素放到一个现实范围内超出范围直接置为缺失。温度在 -80 到 60 摄氏度之间逐小时降水在 0 到 100 毫米之间日降水在 0 到 800 毫米之间风速在 0 到 100 米每秒之间——这是我自己常用的初始范围具体阈值要根据气候区调整。import numpy as np import pandas as pd def clean_obs(df): df df.copy() # 物理范围检查超出合理范围视作观测异常置为 NaN df.loc[(df[t2m] -80) | (df[t2m] 60), t2m] np.nan df.loc[(df[precip] 0) | (df[precip] 800), precip] np.nan # 温度可以用时间插值补短缺口limit 限制最多补 3 个时次 df df.sort_values(time).set_index(time) df[t2m] df[t2m].interpolate(methodtime, limit3) df df.reset_index() # 降水不做插值缺测保持 NaN统计时再决定如何处理 return df温度的interpolate(methodtime)是按时间间隔做线性插值两个时次间隔 6 小时和间隔 1 小时插出来的权重不一样比普通interpolate()更合理。limit3是为了防止大段缺测被“填出”一条平滑的假曲线。降水不插值的原因很直接降水是累积量相邻两个小时一个下了 50 毫米一个下 0 毫米是正常天气插值会凭空制造一场不存在的雨。宁可让统计时段缺一个值也不能编一个值。3.2 把 UTC 时间转成北京时再做逐日统计20-20 时口径怎么对齐气象观测资料的标准时间一般是 UTC中国气象台的发布口径是北京时UTC8。如果把 UTC 时间直接按自然日统计逐日最高温会整体错位更典型的是日降水量和气象台发布值差出一大截。这里有个必须知道的气象行业惯例中国的地面气象观测日界是北京时间 20 时也就是说“今天”的日降水量统计的是今天 20 时到明天 20 时的累积值不是 0 点到 0 点。# 先把无时区的 UTC 时间转成北京时 df[time_beijing] df[time].dt.tz_localize(UTC).dt.tz_convert(Asia/Shanghai) # 北京时 20-20 时气象日把时间轴往前挪 20 小时再按自然日分组 df[day_20] (df[time_beijing] - pd.Timedelta(hours20)).dt.date precip_day df.groupby([station_id, day_20])[precip].sum(min_count1) # 普通自然日0-0时用于温度极值统计 df[day_nat] df[time_beijing].dt.date temp_day df.groupby([station_id, day_nat])[t2m].agg([max, min, mean])代码里的tz_localize(UTC)是在告诉 pandas“这个时间列是无时区的 UTC 时间”tz_convert才是真正转时区顺序不能反。min_count1是降水统计的关键参数如果某一天所有时次都是缺测sum()会返回 0看似正常实则骗人加上min_count1后缺测日返回 NaN至少提醒你这天数据不可信。至于 20-20 时口径把北京时往前挪 20 小时再取日期等价于把 20 点当作新一天的开始这是一劳永逸的做法。3.3 从离散站点到网格场IDW 插值与搜索半径两个参数站点观测是离散的画等值线、算区域平均、和格点数据做对比时都需要把站点插值到规则网格上。气象上插值方法很多克里金、自然邻域、IDW 都很常见但工程上我大部分时候先用反距离权重插值因为它透明、可解释、没有黑匣子。import numpy as np def idw_interp(lon_s, lat_s, val_s, lon_g, lat_g, power2.0, radius1.2): 站点到网格的反距离权重插值。 lon_s/lat_s/val_s: 站点经度、纬度、值一维数组 lon_g/lat_g: 目标网格经度/纬度二维 meshgrid power: 距离权重指数越大越受近点影响 radius: 搜索半径度超过该距离的站点不参与该格点计算 out np.full(lon_g.shape, np.nan) for i in range(lon_g.shape[0]): for j in range(lon_g.shape[1]): d np.sqrt((lon_s - lon_g[i, j]) ** 2 (lat_s - lat_g[i, j]) ** 2) ok np.where(d radius)[0] # 只取半径内的站点 if len(ok) 0: continue w 1.0 / (d[ok] ** power 1e-8) out[i, j] np.sum(w * val_s[ok]) / np.sum(w) return out这个版本是教学级的几百个站点足够用上千个站点建议改用scipy.spatial.cKDTree做近邻搜索不然双重循环会很慢。两个参数要格外上心。power越大插值结果越贴近近处站点但容易出现“牛眼”即站点周围一圈和周围区域差别过分突兀radius是挡住远处站点影响的唯一防线气象上站点距常达几十公里不设半径的话西藏的站点都可能影响到华北格点。后面第五节我会展开讲一个血泪教训不设 radius 插出来的高温中心在实况里根本找不到对应站点。4. 气象指标怎么算风场换算、温湿派生量与滑动窗口预警清洗干净之后真正开始算气象指标。这个阶段的项目价值最高因为下游业务——无论是预警、复盘还是特征工程——都用这里产出的指标。气象指标计算有两个特点一是约定俗成的物理概念多风向风速、露点、相对湿度各有严格定义二是很多指标不是单点时次的而是滑动窗口和阈值的组合这部分最容易写出“能跑但算错”的代码。4.1 从 U/V 分量求风速风向气象风是来向不是去向数值预报和再分析产品里风场一般用 U/V 分量表示U 是东西方向分量正值表示空气向东运动即西风V 是南北方向分量正值表示空气向北运动即南风。很多第一次接触格点风场的人会困惑“为什么没有直接给风速和风向”其实换算只差一行公式真正的坑在方向上。import numpy as np def u_v_to_wind(u, v): u/v 分量转风速和气象风向。 气象风向定义0 度为正北顺时针增加表示风的来向。 例如返回 270 表示风从西边来即西风。 spd np.sqrt(u**2 v**2) # arctan2(v, u) 得到数学角x 轴正向为 0 度逆时针 # 270 - 角度 转成气象方位角%360 防止出现负值 wdir (270 - np.degrees(np.arctan2(v, u))) % 360 return spd, wdir270 - arctan2这个表达式几乎就是气象风场换算的“标准答案”但前提是搞清楚来向和去向。民用领域习惯说“风向东南风”指来向航海和部分工程领域用去向两者相差 180 度。如果用错了风玫瑰图上会出现两个“北”整个图没法看。另外注意气象文件里有时直接提供wind_from_direction和wind_speed两个字段比如大部分地面站点 CSV这时候就别绕 U/V 了直接用现成字段。4.2 用 MetPy 算露点与相对湿度单位要先量化湿度相关指标是农业、仓储、电力负荷分析的常客。单有温度和相对湿度还不够很多场景需要露点温度它比相对湿度更稳定地反映空气中实际水汽含量。MetPy 是气象计算库里的标准件但它对单位极其敏感这是新手最容易踩的坑。import metpy.calc as mpcalc from metpy.units import units # 温度 31.5℃相对湿度 65% temperature 31.5 * units.degC rh 65 * units.percent # 由温度和相对湿度求露点温度 td mpcalc.dewpoint_from_relative_humidity(temperature, rh) print(td.to(degC))MetPy 的底层是 pint 单位库传入的每一个数都必须带单位。直接传31.5会报错或者给出一个默认单位结果完全不可信。处理批量数据时先把 xarray 数据集用.metpy.quantify()统一加上单位再传给 MetPy 的计算函数比在循环里逐个* units.degC干净得多。算出来的露点温度同样带单位输出时.to(degC)转一下养成这个习惯就不会出现“单位对不上”的玄学问题。4.3 高温与强降水预警滑动窗口和最小有效时次怎么设业务预警指标很少只看单个时次连续高温、24 小时累积强降水这类事件本质上都是滑动窗口计算。这里的核心参数有两个窗口长度和最小有效时次。import pandas as pd # 场景逐小时站点数据算连续 3 天日最高温 35℃ 的高温过程 # 先按站排序再按天聚合出日最高温 df df.sort_values([station_id, time_beijing]) daily df.groupby([station_id, day_nat])[t2m].max().rename(tmax).reset_index() # 判断单日是否高温 daily[hot] daily[tmax] 35 # 连续 3 日滑动窗口窗口内高温日数 daily[hot_3d] daily.groupby(station_id)[hot].rolling(3, min_periods3).sum().reset_index(dropTrue) daily[heat_warning] daily[hot_3d] 1 # 3 天内至少 1 天高温 # 强降水逐小时降水 24 小时滑动累积 50mm df[p24] df.groupby(station_id)[precip].rolling(24, min_periods20).sum().reset_index(dropTrue) df[heavy_rain] df[p24] 50两个参数要重视。min_periods是滑动窗口里最少需要多少个有效观测值比如 24 小时窗口允许缺 4 个小时少于 20 个有效值就返回 NaN避免用半个窗口去判断强降水。rolling的min_periods设太严缺一个值整段作废设太松计算结果又会偏小。我一般按窗口的 80% 作为默认值。另外滚动计算必须按站点分组后再做否则两个站的记录会相互串窗口这种错是最难查的。5. 气象分析的五个常见坑时区、缺测、插值与投影的排查记录这一节全是真金白银的踩坑记录。气象数据分析项目里算法模型通常不是瓶颈数据口径和地理投影才是反复翻车的地方。下面五个问题我都在不同项目里遇到过按“现象、原因、解决”写出来照着排查能省不少时间。5.1 温度曲线整体平移 8 小时本地时区和 UTC 没对齐现象画出的逐小时温度曲线最高温出现在北京时间 20 点最低温出现在上午 10 点和实际天气过程完全对不上。原因数据里的时间是 UTC没有转成北京时。UTC 比北京时慢 8 小时最高温实际出现在北京时 14-16 点但原始数据里那个时刻记录的是 UTC 6-8 点画图时没转换曲线整体向右偏了 8 小时。解决读入时间列后立刻做时区处理不要等画图再处理。先把无时区时间tz_localize(UTC)再tz_convert(Asia/Shanghai)然后才做重采样和画图。如果一个项目里既有 UTC 数据又有北京时数据建议所有中间变量统一用北京时只在数据入口和出口做转换。5.2 等值线图出现“不存在”的高温中心IDW 半径没限制现象用站点温度插值出网格场画等值线后发现山区出现一个 40℃ 的高温中心但周边几十公里内根本没有超过 35℃ 的站点。原因IDW 插值时没有设搜索半径网格上每个点都会用到全区域所有站点。当某个网格点周围站点稀疏时几百公里外的站点会被赋予较高权重插出一个实况中不存在的极值中心。这是我早期做插值吃过的大亏一度认为是插值算法“玄学”。解决回到第三节的idw_interp给radius赋一个合理的度数值。我自己在区域分析中常用 1.2-1.5 度作为默认搜索半径站点密集地区可以缩小到 0.5 度。另外在插值前把计算区域裁剪到站点凸包内凸包外的格点直接置空避免外推到无站点区域。5.3 GRIB 文件打不开pandas 读不了二进制气象格式现象拿到 GFS 的 GRIB2 文件习惯性执行pd.read_csv(gfs.t00z.pgrb2.0p25.f000)结果全是乱码或者直接报ParserError。原因GRIB2 是二进制自描述格式不是文本文件pandas 不是干这个的。GRIB 文件内部包含多层嵌套的消息需要 eccodes 库解码后才能读到变量数组。解决用 xarray 配 cfgrib 引擎打开如第二节所示。安装层面优先conda install -c conda-forge cfgrib eccodes把整个解码环境一次性装齐。如果拿到的是老式 GRIB1先把后缀确认清楚cfgrib 对 GRIB1 支持有限可以先用wgrib2转成 NetCDF 再继续加工。5.4 日累计降水与发布值对不上日切时刻和缺测处理现象自己按自然日统计的日降水量和气象台发布的日值总是差几个毫米有时差出十几毫米且差值没有固定规律。原因两个问题叠加。第一是日切时刻不一致气象台按北京时 20-20 时统计自己按自然日 0-0 时统计一场从昨天 19 点下到 21 点的雨会被划分到不同日期。第二是缺测值没处理干净个别站点在降水时段上报 -999直接求和处理时把整日累计拉成负数或者拉低平均值。解决把时间轴往前挪 20 小时再按自然日切代码在第三节已经给出。同时检查降水列的缺测码是否全部转成了 NaN求和时加min_count1并在输出里附带“当日有效时次”字段。有效时次不足 80% 的日期宁可标成可疑也不要放进统计结论。5.5 风玫瑰出现两个“北”来向与去向、数学角与气象角混用现象用同一份风数据画风玫瑰图用 matplotlib 自带极坐标画出来有两个大向前的“北”风向频率分布看起来完全不对。原因两个约定混在一起。第一是风向定义气象风是“来向”即图中指向北的方向代表风从北方吹来而很多绘图库默认使用“去向”。第二是角度系统matplotlib 极坐标默认 0 度在右边东方逆时针增加气象上 0 度在北、顺时针增加。两套方向系统不统一图就废了。解决绘图前先用u_v_to_wind把风场统一成气象来向然后在 matplotlib 极坐标里做两个设置set_theta_offset(np.pi / 2)让 0 度转到正北set_theta_direction(-1)让角度顺时针增加。如果是用 windrose 库它对气象风向已有约定直接传标准气象风向角0-360北为 0即可但传入前同样要确认数据是来向而不是去向。6. 从分析到常规运行回算验证、定时任务和可回溯的数据习惯分析脚本跑通只是第一步气象项目最终都要进入常规运行每天自动拉数据、跑管道、出结论、推送下游。这个阶段我在意三件事结果验证、定时调度、可回溯性。6.1 用一段历史回算做回归验证RMSE 和命中率一起看任何一个气象指标管道上线之前我都会用过去 30-60 天的数据做一次回算把管道算出的结果和气象台发布值对比。只算 RMSE 不够业务预警更关心命中率两个指标一起看。import numpy as np import pandas as pd # result 是管道回算的逐日最高温published 是气象台发布的实况值 # 实际使用时从 DataFrame 里取两列 result pd.Series({... 日最高温: ...}) published pd.Series({... 日最高温: ...}) # RMSE反映整体偏差 rmse float(np.sqrt(((result - published) ** 2).mean())) # 高温命中率阈值为 35℃两边同时触发的比例 hit_rate float(((result 35) (published 35)).mean()) print(fRMSE{rmse:.1f}℃ 高温命中率{hit_rate:.0%})回算还有一个额外价值它能暴露数据管道里隐含的时间偏移。如果 RMSE 数值不大但每日偏差在固定方向先检查时区如果命中率低得离谱检查日切时刻。我遇到过回算 RMSE 只有 1.5℃ 但高温命中率不到 60% 的情况最后查出来是官方口径用的日最高温统计时段和我差了 4 个小时。6.2 定时任务与产物输出文件命名带上数据时间和处理版本气象数据是强时效性的GFS 每天四个时次站点观测每天零点后更新运行计划要卡着数据可用时间。凌晨 2 点半到 3 点是常见窗口前一天观测基本齐全当天最新预报也出了 00 时次。# crontab 示例每天凌晨 2:30 运行气象分析管道 30 2 * * * cd /opt/weather python run_pipeline.py logs/run_$(date \%F).log 21产物文件的命名我坚持带两个信息数据时间和管道版本。比如forecast_20240701_t00_v2.3.csv后面谁用这份文件都能一眼看出它是哪个预报时次、用什么代码跑的。管道脚本一定要进 git每次改阈值、改插值参数都提交一次。这两个习惯合在一起相当于给每一次分析买了后悔药——出问题能精确回到某个版本的代码和数据上。6.3 一个值得养成的习惯原始文件别删处理脚本进 git我见过太多次这种场景分析结果对不上想回溯数据结果原始下载文件已经被清理脚本删了只剩处理后的中间表中间表又没记来源整个数据链变成黑匣子。我现在的做法是原始文件下载后放raw/目录按数据来源和日期分目录存放只追加不删除处理脚本每一个阶段一个文件输出文件名带输入文件的时间范围管道跑完自动写一个meta.json记录输入文件、运行时间、脚本版本和关键参数。这个习惯在长周期项目里救过我很多次。气象分析这件事技术难度往往不在模型而在数据全链路能不能说清楚。先保证每一步可回溯再谈优化算法——把原始数据、处理脚本、版本记录都留住你做的每一版分析才对下游有意义。希望帮到你。本文还有配套的精品资源点击获取