伏羲天气预报实操手册output.nc时间维度重采样与滑动平均平滑处理1. 引言从预报结果到实用数据当你第一次运行伏羲FuXi天气预报系统看到终端里滚动着“预报完成”的提示然后生成了一个output.nc文件时是不是既兴奋又有点迷茫兴奋的是这个强大的AI模型真的生成了未来15天的全球天气数据迷茫的是面对这个包含多个时间步、多个变量的NetCDF文件下一步该怎么做output.nc文件就像一盒未经切割的钻石原石——它包含了所有有价值的信息但直接使用可能会遇到几个实际问题时间分辨率不匹配伏羲默认每6小时输出一个预报结果但你的分析可能需要日平均、周平均数据预报“跳跃”现象单个时间步的预报可能存在随机波动影响趋势判断数据量太大15天预报360小时会产生60个时间步直接分析效率低下这就是为什么我们需要对output.nc进行后处理。本文将手把手教你两种最实用的处理方法时间维度重采样和滑动平均平滑。学完这篇教程你将能够将6小时间隔的预报数据转换为日、周等任意时间分辨率平滑预报曲线消除随机波动看清真实趋势提取关键气象要素的时间序列用于进一步分析和可视化无论你是气象专业的学生、气候研究人员还是对AI天气预报感兴趣的开发者这篇实操手册都能让你快速掌握伏羲预报结果的实用化处理技巧。2. 理解伏羲的输出数据结构在开始处理之前我们需要先搞清楚output.nc里面到底有什么。这就像你要整理一个仓库得先知道里面放了哪些东西、怎么放的。2.1 打开并查看output.nc文件首先让我们用Python打开这个NetCDF文件看看它的结构import xarray as xr # 打开伏羲生成的预报结果文件 file_path /root/fuxi2/output.nc # 你的文件路径 ds xr.open_dataset(file_path) # 查看数据集的基本信息 print(数据集维度信息) print(ds.dims) print(\n数据集包含的变量) for var_name in ds.data_vars: var ds[var_name] print(f- {var_name}: {var.dims}, 形状: {var.shape})运行这段代码你会看到类似这样的输出数据集维度信息 OrderedDict([(time, 60), (lat, 721), (lon, 1440)]) 数据集包含的变量 - Z: (time, lat, lon), 形状: (60, 721, 1440) - T: (time, lat, lon), 形状: (60, 721, 1440) - U: (time, lat, lon), 形状: (60, 721, 1440) - V: (time, lat, lon), 形状: (60, 721, 1440) - R: (time, lat, lon), 形状: (60, 721, 1440) - T2M: (time, lat, lon), 形状: (60, 721, 1440) - U10: (time, lat, lon), 形状: (60, 721, 1440) - V10: (time, lat, lon), 形状: (60, 721, 1440) - MSL: (time, lat, lon), 形状: (60, 721, 1440) - TP: (time, lat, lon), 形状: (60, 721, 1440)2.2 关键维度解读理解这些维度的含义很重要维度含义数值说明time时间维度60个时间步每步6小时共15天预报lat纬度维度721个格点从南纬90°到北纬90°分辨率0.25°lon经度维度1440个格点从东经0°到360°分辨率0.25°时间维度的特别说明 伏羲的预报从初始时刻开始每6小时输出一次结果。如果你在2024年1月1日00:00启动预报那么time[0] 对应 2024-01-01 00:00time[1] 对应 2024-01-01 06:00time[2] 对应 2024-01-01 12:00... 以此类推time[59] 对应 2024-01-15 18:002.3 变量含义速查表output.nc包含了10个气象变量为了方便后续处理这里列出它们的含义变量名全称单位物理意义Z位势高度m²/s²大气压力场的垂直分布反映天气系统强度T温度K大气温度13个垂直层次UU风分量m/s东西方向的风速VV风分量m/s南北方向的风速R相对湿度%空气湿度影响云和降水T2M2米温度K地面以上2米处的气温与人体感受最相关U1010米U风m/s地面以上10米处的东西风V1010米V风m/s地面以上10米处的南北风MSL海平面气压Pa修正到海平面的大气压力TP总降水量m6小时累积降水量现在我们对数据有了基本了解接下来进入实战环节。3. 方法一时间维度重采样重采样就像给数据“换一个时钟”。伏羲每6小时报一次天气但你可能更关心“明天最高温多少”或者“这周平均降水量多大”。重采样就是解决这个问题的。3.1 什么是时间重采样简单来说重采样就是改变数据的时间频率。比如降采样从高频到低频6小时→24小时升采样从低频到高频24小时→6小时对于天气预报结果我们主要用降采样因为减少数据量提高处理速度得到更有意义的时间统计日平均、周平均匹配其他数据源的时间分辨率3.2 基础重采样操作让我们从一个简单的例子开始计算北京地区北纬40°东经116°的日平均温度。import numpy as np import xarray as xr import pandas as pd # 1. 打开数据 ds xr.open_dataset(/root/fuxi2/output.nc) # 2. 设置时间坐标重要 # 假设预报开始时间是2024-01-01 00:00 start_time pd.Timestamp(2024-01-01 00:00:00) time_coords pd.date_range(startstart_time, periods60, freq6H) ds ds.assign_coords(timetime_coords) # 3. 提取北京附近的点近似位置 beijing_lat 40.0 # 北纬40度 beijing_lon 116.0 # 东经116度 # 找到最近的格点 lat_idx np.abs(ds.lat.values - beijing_lat).argmin() lon_idx np.abs(ds.lon.values - beijing_lon).argmin() # 4. 提取2米温度时间序列 t2m_beijing ds[T2M].isel(latlat_idx, lonlon_idx) print(北京2米温度原始数据前5个时间步) print(t2m_beijing[:5].values) print(f时间范围{t2m_beijing.time[0].values} 到 {t2m_beijing.time[-1].values})运行后你会看到类似这样的输出北京2米温度原始数据前5个时间步 [275.3 274.8 273.9 273.2 272.8] 时间范围2024-01-01T00:00:00.000000000 到 2024-01-03T06:00:00.0000000003.3 实现日平均重采样现在我们来计算日平均温度# 5. 重采样为日平均 # resample() 方法按时间频率分组然后mean()计算平均值 t2m_daily t2m_beijing.resample(time1D).mean() print(\n日平均温度) for i in range(min(5, len(t2m_daily))): # 显示前5天 date_str pd.Timestamp(t2m_daily.time[i].values).strftime(%Y-%m-%d) temp_k t2m_daily[i].values temp_c temp_k - 273.15 # 开尔文转摄氏度 print(f{date_str}: {temp_c:.1f}°C)输出示例日平均温度 2024-01-01: 2.5°C 2024-01-02: 1.8°C 2024-01-03: 1.2°C 2024-01-04: 0.9°C 2024-01-05: 0.5°C3.4 多种重采样频率示例不同的分析需求需要不同的时间分辨率。下面这个表格总结了常用的重采样方法分析目的重采样频率代码示例结果含义日变化分析6小时→24小时日平均.resample(time1D).mean()每天的平均状况周趋势分析6小时→7天周平均.resample(time7D).mean()每周的平均变化极端事件6小时→24小时日最大.resample(time1D).max()每天的最高值累积降水6小时→24小时日总和.resample(time1D).sum()每天的降水量月气候值6小时→月平均.resample(time1M).mean()月平均气候完整的多变量重采样示例def resample_forecast_data(input_path, output_path, freq1D, methodmean): 对伏羲预报结果进行时间重采样 参数 input_path: 输入nc文件路径 output_path: 输出nc文件路径 freq: 重采样频率如 1D日、7D周、1H小时 method: 统计方法mean平均、max最大、min最小、sum求和 # 打开数据集 ds xr.open_dataset(input_path) # 设置时间坐标如果还没有 if time not in ds.coords: start_time pd.Timestamp.now().replace(hour0, minute0, second0) time_coords pd.date_range(startstart_time, periodslen(ds.time), freq6H) ds ds.assign_coords(timetime_coords) # 根据method选择重采样方法 if method mean: ds_resampled ds.resample(timefreq).mean() elif method max: ds_resampled ds.resample(timefreq).max() elif method min: ds_resampled ds.resample(timefreq).min() elif method sum: ds_resampled ds.resample(timefreq).sum() else: raise ValueError(f不支持的method: {method}) # 保存结果 ds_resampled.to_netcdf(output_path) print(f重采样完成文件已保存到: {output_path}) print(f原始时间步数: {len(ds.time)}重采样后: {len(ds_resampled.time)}) return ds_resampled # 使用示例生成日平均预报 resampled_data resample_forecast_data( input_path/root/fuxi2/output.nc, output_path/root/fuxi2/output_daily_mean.nc, freq1D, methodmean )3.5 重采样的注意事项在实际使用中有几个细节需要注意时间坐标必须正确如果数据没有正确的时间坐标重采样会失败。确保你的time维度是datetime类型。边界处理重采样时时间窗口的边界会影响结果。xarray默认使用左闭右开区间即[00:00, 24:00)。缺失值处理如果某个时间窗口内全是缺失值重采样结果也会是缺失值。可以使用.resample().mean(skipnaTrue)跳过缺失值。内存考虑对全球数据进行重采样可能消耗大量内存。如果只需要特定区域可以先裁剪再重采样# 先裁剪中国区域再重采样节省内存 china_region ds.sel(latslice(15, 55), lonslice(70, 140)) china_daily china_region.resample(time1D).mean()4. 方法二滑动平均平滑处理如果你观察过原始的预报时间序列可能会发现曲线有些毛刺——温度忽高忽低气压时升时降。这些短期波动很多时候是预报的不确定性造成的而不是真实的天气变化。滑动平均可以帮助我们平滑这些波动看清真正的趋势。4.1 为什么需要滑动平均想象一下你看股票走势图每天的股价上上下下很难判断是涨是跌。但如果你看5日平均线、20日平均线趋势就清晰多了。天气预报数据也是同样的道理。滑动平均的三个主要作用滤除高频噪声消除6小时尺度上的随机波动凸显趋势让日变化、季节变化等长期趋势更明显提高可视化效果平滑的曲线更容易理解和解释4.2 滑动平均的基本原理滑动平均很简单对于时间序列上的每个点用它前后一段时间内的平均值来代替它自己。举个例子3点滑动平均原始数据 [10, 12, 11, 13, 14, 12, 15] 滑动平均 [ - , 11, 12, 12.67, 13, 13.67, - ]计算过程第二个点(12) (101211)/3 11第三个点(11) (121113)/3 12以此类推...4.3 实现滑动平均平滑让我们用Python实现滑动平均并比较平滑前后的效果import matplotlib.pyplot as plt def sliding_average(data, window_size5): 计算滑动平均 参数 data: 输入数据numpy数组或xarray DataArray window_size: 滑动窗口大小必须是奇数 返回 平滑后的数据两端(window_size//2)个点用NaN填充 if window_size % 2 0: window_size 1 # 确保是奇数 print(f窗口大小调整为奇数: {window_size}) half_window window_size // 2 smoothed np.full_like(data, np.nan) # 用NaN初始化 # 计算滑动平均 for i in range(half_window, len(data) - half_window): smoothed[i] np.mean(data[i-half_window : ihalf_window1]) return smoothed # 提取一段温度时间序列 temp_series ds[T2M].isel(latlat_idx, lonlon_idx).values # 应用不同窗口大小的滑动平均 window_sizes [3, 7, 15] # 对应18小时、42小时、90小时 smoothed_results {} for ws in window_sizes: smoothed_results[ws] sliding_average(temp_series, window_sizews) # 可视化比较 plt.figure(figsize(12, 6)) time_hours np.arange(len(temp_series)) * 6 # 转换为小时 plt.plot(time_hours, temp_series, b-, alpha0.3, label原始数据, linewidth1) colors [r-, g-, m-] for idx, ws in enumerate(window_sizes): plt.plot(time_hours, smoothed_results[ws], colors[idx], labelf{ws}点滑动平均 ({ws*6}小时), linewidth2) plt.xlabel(预报时间 (小时)) plt.ylabel(2米温度 (K)) plt.title(不同滑动窗口大小的平滑效果对比) plt.legend() plt.grid(True, alpha0.3) plt.tight_layout() plt.savefig(/root/fuxi2/smoothing_comparison.png, dpi150) plt.show()运行这段代码你会看到一张对比图清晰地展示不同窗口大小的平滑效果。4.4 xarray内置的滚动方法xarray提供了更便捷的.rolling()方法来实现滑动平均而且可以处理多维数据def smooth_forecast_data(ds, window_size5, dimtime): 使用xarray的rolling方法平滑预报数据 参数 ds: xarray Dataset 或 DataArray window_size: 滑动窗口大小 dim: 沿着哪个维度平滑通常是time 返回 平滑后的数据 # 使用rolling创建滑动窗口然后计算平均 # centerTrue表示窗口居中 ds_smoothed ds.rolling({dim: window_size}, centerTrue).mean() # rolling方法会在两端留下NaN我们可以选择 # 1. 保留NaN反映边界处置信度低 # 2. 用原始值填充简单处理 return ds_smoothed # 对整个数据集进行平滑 smoothed_ds smooth_forecast_data(ds, window_size7) # 比较平滑前后的北京温度 original_temp ds[T2M].isel(latlat_idx, lonlon_idx) smoothed_temp smoothed_ds[T2M].isel(latlat_idx, lonlon_idx) print(平滑前后对比前10个时间步) print(时间 原始温度(K) 平滑温度(K) 差值(K)) print(- * 50) for i in range(10): time_str pd.Timestamp(original_temp.time[i].values).strftime(%m-%d %H:%M) orig original_temp[i].values smooth smoothed_temp[i].values diff smooth - orig if not np.isnan(smooth) else np.nan print(f{time_str} {orig:7.2f} {smooth:7.2f} {diff:6.2f})4.5 滑动平均的高级技巧4.5.1 加权滑动平均有时候我们想让最近的数据权重更大。这就是加权滑动平均def weighted_sliding_average(data, window_size5): 使用三角形加权的滑动平均 weights np.arange(1, window_size//2 2) weights np.concatenate([weights, weights[-2::-1]]) # 创建对称权重 weights weights / weights.sum() # 归一化 half_window window_size // 2 smoothed np.full_like(data, np.nan) for i in range(half_window, len(data) - half_window): smoothed[i] np.sum(data[i-half_window : ihalf_window1] * weights) return smoothed4.5.2 处理缺失值实际数据中常有缺失值我们需要小心处理def sliding_average_with_nan(data, window_size5, min_valid3): 能处理NaN值的滑动平均 参数 min_valid: 窗口内至少需要多少个有效值才计算平均 half_window window_size // 2 smoothed np.full_like(data, np.nan) for i in range(half_window, len(data) - half_window): window_data data[i-half_window : ihalf_window1] valid_mask ~np.isnan(window_data) if np.sum(valid_mask) min_valid: smoothed[i] np.mean(window_data[valid_mask]) return smoothed4.5.3 选择最佳窗口大小窗口大小怎么选太大则过度平滑失去细节太小则平滑效果不足。这里有个经验法则分析目的推荐窗口大小对应时间尺度理由滤除随机噪声3-5点18-30小时保留日变化滤除短时波动分析日变化5-7点30-42小时平滑半天尺度波动分析趋势7-15点42-90小时看清多日变化趋势周尺度分析15-25点90-150小时滤除日变化凸显周趋势你可以用这段代码测试不同窗口大小def find_optimal_window(data, max_window25): 通过计算平滑度指标寻找最佳窗口大小 results [] for ws in range(3, max_window1, 2): # 只测试奇数 smoothed sliding_average(data, window_sizews) # 计算平滑度指标标准差越小越平滑 valid_data smoothed[~np.isnan(smoothed)] if len(valid_data) 0: smoothness 1 / np.std(valid_data) # 平滑度指标 results.append((ws, smoothness)) # 按平滑度排序 results.sort(keylambda x: x[1], reverseTrue) print(窗口大小推荐按平滑度排序) for ws, smoothness in results[:5]: # 显示前5个 print(f 窗口大小 {ws}点 ({ws*6}小时): 平滑度 {smoothness:.3f}) return results[0][0] if results else 5 # 使用示例 optimal_ws find_optimal_window(temp_series, max_window25) print(f\n推荐窗口大小: {optimal_ws}点 ({optimal_ws*6}小时))5. 实战案例完整处理流程现在我们把学到的两种方法结合起来完成一个完整的实战案例分析未来15天中国东部地区的降水趋势。5.1 案例背景与目标假设你是气象分析师需要获取伏羲对未来15天的降水预报将6小时间隔的数据转换为日累积降水平滑日降水序列识别干旱/多雨期生成可视化报告5.2 完整代码实现import xarray as xr import numpy as np import pandas as pd import matplotlib.pyplot as plt from matplotlib import cm import cartopy.crs as ccrs import cartopy.feature as cfeature def process_precipitation_forecast(input_path, regionNone, output_dir.): 完整处理降水预报的流程 参数 input_path: 伏羲输出文件路径 region: 区域范围 [lat_min, lat_max, lon_min, lon_max]默认中国东部 output_dir: 输出目录 # 1. 设置区域默认中国东部 if region is None: region [20, 45, 100, 125] # 中国东部 # 2. 打开数据并裁剪区域 print(步骤1: 加载数据并裁剪区域...) ds xr.open_dataset(input_path) # 设置时间坐标 start_time pd.Timestamp.now().replace(hour0, minute0, second0) time_coords pd.date_range(startstart_time, periodslen(ds.time), freq6H) ds ds.assign_coords(timetime_coords) # 裁剪区域 ds_region ds.sel( latslice(region[0], region[1]), lonslice(region[2], region[3]) ) # 3. 重采样计算日累积降水 print(步骤2: 计算日累积降水...) # TP是6小时累积降水单位是米转换为毫米 precipitation_mm ds_region[TP] * 1000 # 米→毫米 # 重采样为日累积求和 daily_precip precipitation_mm.resample(time1D).sum() # 4. 区域平均计算整个区域的平均日降水 print(步骤3: 计算区域平均...) # 使用面积加权平均考虑纬度变化 lat_weights np.cos(np.deg2rad(daily_precip.lat)) weighted_precip daily_precip.weighted(lat_weights) region_mean weighted_precip.mean(dim[lat, lon]) # 5. 滑动平均平滑 print(步骤4: 应用滑动平均平滑...) window_size 3 # 3天滑动平均 smoothed_mean region_mean.rolling(timewindow_size, centerTrue).mean() # 6. 保存处理结果 print(步骤5: 保存结果...) # 创建包含所有结果的Dataset result_ds xr.Dataset({ daily_precip: daily_precip, region_mean: region_mean, smoothed_mean: smoothed_mean }) output_path f{output_dir}/processed_precipitation.nc result_ds.to_netcdf(output_path) # 7. 生成可视化 print(步骤6: 生成可视化...) generate_precipitation_plots(result_ds, region, output_dir) print(f\n处理完成) print(f原始文件: {input_path}) print(f处理结果: {output_path}) return result_ds def generate_precipitation_plots(result_ds, region, output_dir): 生成降水分析的可视化图表 # 图1区域平均降水时间序列 plt.figure(figsize(12, 8)) # 子图1原始日降水与平滑曲线 plt.subplot(2, 2, 1) time_dates pd.to_datetime(result_ds[region_mean].time.values) plt.bar(time_dates, result_ds[region_mean].values, alpha0.5, label日降水, colorskyblue) plt.plot(time_dates, result_ds[smoothed_mean].values, r-, linewidth2, labelf{3}天滑动平均) plt.xlabel(日期) plt.ylabel(区域平均降水 (mm/天)) plt.title(中国东部区域平均日降水) plt.legend() plt.grid(True, alpha0.3) plt.xticks(rotation45) # 子图2累积降水 plt.subplot(2, 2, 2) cumulative_precip np.cumsum(result_ds[region_mean].values) plt.plot(time_dates, cumulative_precip, g-, linewidth2) plt.fill_between(time_dates, 0, cumulative_precip, alpha0.3, colorgreen) plt.xlabel(日期) plt.ylabel(累积降水 (mm)) plt.title(15天累积降水量) plt.grid(True, alpha0.3) plt.xticks(rotation45) # 子图3降水空间分布第7天示例 plt.subplot(2, 2, 3, projectionccrs.PlateCarree()) day_idx 6 # 第7天 precip_day result_ds[daily_precip].isel(timeday_idx) # 创建地图 ax plt.gca() ax.set_extent([region[2], region[3], region[0], region[1]], crsccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle:) ax.add_feature(cfeature.LAKES, alpha0.5) ax.add_feature(cfeature.RIVERS) # 绘制降水 im precip_day.plot(axax, transformccrs.PlateCarree(), cmapYlGnBu, add_colorbarFalse) plt.colorbar(im, axax, orientationhorizontal, pad0.05, label降水量 (mm)) date_str pd.Timestamp(time_dates[day_idx]).strftime(%Y-%m-%d) plt.title(f{date_str} 降水空间分布) # 子图4降水强度统计 plt.subplot(2, 2, 4) # 计算不同强度降水的天数 precip_values result_ds[region_mean].values bins [0, 1, 5, 10, 20, 50] # 降水强度分级 (mm/天) labels [微量(1), 小雨(1-5), 中雨(5-10), 大雨(10-20), 暴雨(20)] counts, _ np.histogram(precip_values, binsbins) colors [lightgray, lightblue, blue, darkblue, purple] bars plt.bar(range(len(counts)), counts, colorcolors) plt.xlabel(降水强度) plt.ylabel(天数) plt.title(降水强度分布) plt.xticks(range(len(labels)), labels, rotation45) # 在柱子上添加数值 for bar, count in zip(bars, counts): height bar.get_height() plt.text(bar.get_x() bar.get_width()/2., height, f{count}, hacenter, vabottom) plt.tight_layout() plt.savefig(f{output_dir}/precipitation_analysis.png, dpi150, bbox_inchestight) plt.close() # 生成文本报告 generate_text_report(result_ds, output_dir) def generate_text_report(result_ds, output_dir): 生成文本分析报告 mean_precip result_ds[region_mean].values smoothed_precip result_ds[smoothed_mean].values report f伏羲降水预报分析报告 生成时间: {pd.Timestamp.now().strftime(%Y-%m-%d %H:%M:%S)} 分析时段: {pd.Timestamp(result_ds.time[0].values).strftime(%Y-%m-%d)} 至 {pd.Timestamp(result_ds.time[-1].values).strftime(%Y-%m-%d)} 数据来源: FuXi 15天全球天气预报系统 【统计摘要】 - 平均日降水: {np.nanmean(mean_precip):.2f} mm/天 - 最大日降水: {np.nanmax(mean_precip):.2f} mm/天 (第{np.nanargmax(mean_precip)1}天) - 最小日降水: {np.nanmin(mean_precip):.2f} mm/天 (第{np.nanargmin(mean_precip)1}天) - 总累积降水: {np.nansum(mean_precip):.2f} mm - 降水天数(1mm): {np.sum(mean_precip 1)} 天 - 干旱天数(0.5mm): {np.sum(mean_precip 0.5)} 天 【趋势分析】 # 添加趋势判断 if len(smoothed_precip) 5: first_half np.nanmean(smoothed_precip[:5]) second_half np.nanmean(smoothed_precip[-5:]) if second_half first_half * 1.2: report - 降水呈增加趋势后期可能更湿润\n elif second_half first_half * 0.8: report - 降水呈减少趋势后期可能更干燥\n else: report - 降水趋势相对平稳\n # 识别关键特征日 report \n【关键特征日】\n # 找峰值和谷值 valid_indices ~np.isnan(smoothed_precip) if np.sum(valid_indices) 0: smooth_valid smoothed_precip[valid_indices] time_valid result_ds.time.values[valid_indices] # 找局部最大值峰值 from scipy.signal import find_peaks peaks, _ find_peaks(smooth_valid, heightnp.nanmean(smooth_valid)) for idx in peaks[:3]: # 显示前3个峰值 date_str pd.Timestamp(time_valid[idx]).strftime(%Y-%m-%d) report f- {date_str}: 降水峰值预计{smooth_valid[idx]:.1f} mm/天\n # 找局部最小值谷值 valleys, _ find_peaks(-smooth_valid) for idx in valleys[:3]: # 显示前3个谷值 date_str pd.Timestamp(time_valid[idx]).strftime(%Y-%m-%d) report f- {date_str}: 降水谷值预计{smooth_valid[idx]:.1f} mm/天\n report \n【使用建议】\n report 1. 本分析基于伏羲AI天气预报仅供参考\n report 2. 滑动平均帮助识别趋势但会平滑短期变化\n report 3. 建议结合其他预报资料进行综合判断\n # 保存报告 report_path f{output_dir}/precipitation_report.txt with open(report_path, w, encodingutf-8) as f: f.write(report) print(f分析报告已保存: {report_path}) # 运行完整流程 if __name__ __main__: result process_precipitation_forecast( input_path/root/fuxi2/output.nc, output_dir/root/fuxi2/results )5.3 案例输出解读运行上面的代码后你会得到processed_precipitation.nc包含处理后的NetCDF数据precipitation_analysis.png四合一可视化图表precipitation_report.txt文本分析报告报告示例内容伏羲降水预报分析报告 生成时间: 2024-01-15 14:30:00 分析时段: 2024-01-16 至 2024-01-30 【统计摘要】 - 平均日降水: 3.2 mm/天 - 最大日降水: 8.5 mm/天 (第7天) - 最小日降水: 0.2 mm/天 (第3天) - 总累积降水: 48.3 mm - 降水天数(1mm): 9天 - 干旱天数(0.5mm): 4天 【趋势分析】 - 降水呈增加趋势后期可能更湿润 【关键特征日】 - 2024-01-22: 降水峰值预计7.8 mm/天 - 2024-01-25: 降水峰值预计6.5 mm/天 - 2024-01-18: 降水谷值预计0.8 mm/天 【使用建议】 1. 本分析基于伏羲AI天气预报仅供参考 2. 滑动平均帮助识别趋势但会平滑短期变化 3. 建议结合其他预报资料进行综合判断6. 总结与最佳实践通过本文的学习你已经掌握了伏羲天气预报结果后处理的两大核心技能时间维度重采样和滑动平均平滑。让我们回顾一下关键要点6.1 技术要点回顾重采样改变时间频率将6小时间隔的预报转换为日、周等更有意义的统计量使用.resample()方法选择合适频率1D、7D等根据需求选择统计方法.mean()、.sum()、.max()等滑动平均平滑时间序列滤除随机波动凸显真实趋势使用.rolling()方法设置合适的窗口大小窗口大小通常选3-7点18-42小时平衡细节与平滑度考虑使用加权平均或处理缺失值的变体组合使用效果更佳先重采样降低数据频率再滑动平均平滑曲线是常见的最佳实践6.2 实用建议根据不同的应用场景这里有一些实用建议应用场景推荐处理流程理由气候趋势分析重采样为月平均 → 12月滑动平均滤除季节内波动看清年际变化农业气象服务重采样为日累积 → 7天滑动平均提供周尺度水分供需分析灾害预警保持6小时原始频率 → 3点滑动平均保留短时强降水信号仅滤除噪声能源需求预测重采样为日平均温度 → 5点滑动平均平滑温度曲线用于负荷预测航空天气预报保持原始频率必要时3点滑动平均需要高时间分辨率仅轻微平滑6.3 常见问题与解决Q: 处理后的数据还能用伏羲继续预报吗A: 不能。重采样和平滑会改变数据的时间结构和统计特性不适合作为下一轮预报的初始场。这些处理主要用于分析和可视化。Q: 滑动平均会导致预报滞后吗A: 使用centerTrue的滑动平均不会导致滞后因为它是用前后数据计算当前点的平均值。但边界处开头和结尾由于数据不足平滑效果会打折扣。Q: 如何处理全球数据的内存问题A: 三个策略先裁剪感兴趣的区域再处理使用Dask进行并行处理分变量处理保存中间结果Q: 有现成的工具可以直接用吗A: 你可以将本文的代码封装成函数库方便重复使用。这里提供一个简单的封装示例# fuxi_postprocess.py import xarray as xr import numpy as np import pandas as pd class FuXiPostProcessor: 伏羲预报后处理工具类 def __init__(self, data_path): self.ds xr.open_dataset(data_path) self._set_time_coords() def _set_time_coords(self): 设置时间坐标 if time not in self.ds.coords: start_time pd.Timestamp.now().replace(hour0, minute0, second0) time_coords pd.date_range(startstart_time, periodslen(self.ds.time), freq6H) self.ds self.ds.assign_coords(timetime_coords) def resample(self, freq1D, methodmean, variablesNone): 重采样数据 if variables is None: variables list(self.ds.data_vars) if method mean: return self.ds[variables].resample(timefreq).mean() elif method sum: return self.ds[variables].resample(timefreq).sum() # ... 其他方法 def smooth(self, window_size5, dimtime): 滑动平均平滑 return self.ds.rolling({dim: window_size}, centerTrue).mean() def process_region(self, region, variable, freq1D, window_size3): 完整的区域处理流程 # 裁剪区域 ds_region self.ds.sel( latslice(region[0], region[1]), lonslice(region[2], region[3]) ) # 重采样 ds_resampled ds_region.resample(timefreq).mean() # 区域平均 lat_weights np.cos(np.deg2rad(ds_resampled.lat)) weighted_data ds_resampled[variable].weighted(lat_weights) region_mean weighted_data.mean(dim[lat, lon]) # 平滑 smoothed region_mean.rolling(timewindow_size, centerTrue).mean() return { resampled: ds_resampled, region_mean: region_mean, smoothed: smoothed } # 使用示例 processor FuXiPostProcessor(/root/fuxi2/output.nc) result processor.process_region( region[20, 45, 100, 125], variableT2M, freq1D, window_size3 )6.4 下一步学习建议掌握了基础的后处理技能后你可以进一步探索多模型集成将伏羲预报与其他模型结果结合提高预报可靠性统计后处理使用MOSModel Output Statistics等技术校准预报偏差可视化进阶学习使用Cartopy、Metpy等专业气象可视化库自动化流程将处理流程封装成定时任务实现自动化预报产品生成天气预报的最终价值在于应用。通过对伏羲预报结果的合理后处理你可以将原始的数值输出转化为决策支持信息真正发挥AI天气预报的实用价值。获取更多AI镜像想探索更多AI镜像和应用场景访问 CSDN星图镜像广场提供丰富的预置镜像覆盖大模型推理、图像生成、视频生成、模型微调等多个领域支持一键部署。
伏羲天气预报实操手册：output.nc时间维度重采样与滑动平均平滑处理
伏羲天气预报实操手册output.nc时间维度重采样与滑动平均平滑处理1. 引言从预报结果到实用数据当你第一次运行伏羲FuXi天气预报系统看到终端里滚动着“预报完成”的提示然后生成了一个output.nc文件时是不是既兴奋又有点迷茫兴奋的是这个强大的AI模型真的生成了未来15天的全球天气数据迷茫的是面对这个包含多个时间步、多个变量的NetCDF文件下一步该怎么做output.nc文件就像一盒未经切割的钻石原石——它包含了所有有价值的信息但直接使用可能会遇到几个实际问题时间分辨率不匹配伏羲默认每6小时输出一个预报结果但你的分析可能需要日平均、周平均数据预报“跳跃”现象单个时间步的预报可能存在随机波动影响趋势判断数据量太大15天预报360小时会产生60个时间步直接分析效率低下这就是为什么我们需要对output.nc进行后处理。本文将手把手教你两种最实用的处理方法时间维度重采样和滑动平均平滑。学完这篇教程你将能够将6小时间隔的预报数据转换为日、周等任意时间分辨率平滑预报曲线消除随机波动看清真实趋势提取关键气象要素的时间序列用于进一步分析和可视化无论你是气象专业的学生、气候研究人员还是对AI天气预报感兴趣的开发者这篇实操手册都能让你快速掌握伏羲预报结果的实用化处理技巧。2. 理解伏羲的输出数据结构在开始处理之前我们需要先搞清楚output.nc里面到底有什么。这就像你要整理一个仓库得先知道里面放了哪些东西、怎么放的。2.1 打开并查看output.nc文件首先让我们用Python打开这个NetCDF文件看看它的结构import xarray as xr # 打开伏羲生成的预报结果文件 file_path /root/fuxi2/output.nc # 你的文件路径 ds xr.open_dataset(file_path) # 查看数据集的基本信息 print(数据集维度信息) print(ds.dims) print(\n数据集包含的变量) for var_name in ds.data_vars: var ds[var_name] print(f- {var_name}: {var.dims}, 形状: {var.shape})运行这段代码你会看到类似这样的输出数据集维度信息 OrderedDict([(time, 60), (lat, 721), (lon, 1440)]) 数据集包含的变量 - Z: (time, lat, lon), 形状: (60, 721, 1440) - T: (time, lat, lon), 形状: (60, 721, 1440) - U: (time, lat, lon), 形状: (60, 721, 1440) - V: (time, lat, lon), 形状: (60, 721, 1440) - R: (time, lat, lon), 形状: (60, 721, 1440) - T2M: (time, lat, lon), 形状: (60, 721, 1440) - U10: (time, lat, lon), 形状: (60, 721, 1440) - V10: (time, lat, lon), 形状: (60, 721, 1440) - MSL: (time, lat, lon), 形状: (60, 721, 1440) - TP: (time, lat, lon), 形状: (60, 721, 1440)2.2 关键维度解读理解这些维度的含义很重要维度含义数值说明time时间维度60个时间步每步6小时共15天预报lat纬度维度721个格点从南纬90°到北纬90°分辨率0.25°lon经度维度1440个格点从东经0°到360°分辨率0.25°时间维度的特别说明 伏羲的预报从初始时刻开始每6小时输出一次结果。如果你在2024年1月1日00:00启动预报那么time[0] 对应 2024-01-01 00:00time[1] 对应 2024-01-01 06:00time[2] 对应 2024-01-01 12:00... 以此类推time[59] 对应 2024-01-15 18:002.3 变量含义速查表output.nc包含了10个气象变量为了方便后续处理这里列出它们的含义变量名全称单位物理意义Z位势高度m²/s²大气压力场的垂直分布反映天气系统强度T温度K大气温度13个垂直层次UU风分量m/s东西方向的风速VV风分量m/s南北方向的风速R相对湿度%空气湿度影响云和降水T2M2米温度K地面以上2米处的气温与人体感受最相关U1010米U风m/s地面以上10米处的东西风V1010米V风m/s地面以上10米处的南北风MSL海平面气压Pa修正到海平面的大气压力TP总降水量m6小时累积降水量现在我们对数据有了基本了解接下来进入实战环节。3. 方法一时间维度重采样重采样就像给数据“换一个时钟”。伏羲每6小时报一次天气但你可能更关心“明天最高温多少”或者“这周平均降水量多大”。重采样就是解决这个问题的。3.1 什么是时间重采样简单来说重采样就是改变数据的时间频率。比如降采样从高频到低频6小时→24小时升采样从低频到高频24小时→6小时对于天气预报结果我们主要用降采样因为减少数据量提高处理速度得到更有意义的时间统计日平均、周平均匹配其他数据源的时间分辨率3.2 基础重采样操作让我们从一个简单的例子开始计算北京地区北纬40°东经116°的日平均温度。import numpy as np import xarray as xr import pandas as pd # 1. 打开数据 ds xr.open_dataset(/root/fuxi2/output.nc) # 2. 设置时间坐标重要 # 假设预报开始时间是2024-01-01 00:00 start_time pd.Timestamp(2024-01-01 00:00:00) time_coords pd.date_range(startstart_time, periods60, freq6H) ds ds.assign_coords(timetime_coords) # 3. 提取北京附近的点近似位置 beijing_lat 40.0 # 北纬40度 beijing_lon 116.0 # 东经116度 # 找到最近的格点 lat_idx np.abs(ds.lat.values - beijing_lat).argmin() lon_idx np.abs(ds.lon.values - beijing_lon).argmin() # 4. 提取2米温度时间序列 t2m_beijing ds[T2M].isel(latlat_idx, lonlon_idx) print(北京2米温度原始数据前5个时间步) print(t2m_beijing[:5].values) print(f时间范围{t2m_beijing.time[0].values} 到 {t2m_beijing.time[-1].values})运行后你会看到类似这样的输出北京2米温度原始数据前5个时间步 [275.3 274.8 273.9 273.2 272.8] 时间范围2024-01-01T00:00:00.000000000 到 2024-01-03T06:00:00.0000000003.3 实现日平均重采样现在我们来计算日平均温度# 5. 重采样为日平均 # resample() 方法按时间频率分组然后mean()计算平均值 t2m_daily t2m_beijing.resample(time1D).mean() print(\n日平均温度) for i in range(min(5, len(t2m_daily))): # 显示前5天 date_str pd.Timestamp(t2m_daily.time[i].values).strftime(%Y-%m-%d) temp_k t2m_daily[i].values temp_c temp_k - 273.15 # 开尔文转摄氏度 print(f{date_str}: {temp_c:.1f}°C)输出示例日平均温度 2024-01-01: 2.5°C 2024-01-02: 1.8°C 2024-01-03: 1.2°C 2024-01-04: 0.9°C 2024-01-05: 0.5°C3.4 多种重采样频率示例不同的分析需求需要不同的时间分辨率。下面这个表格总结了常用的重采样方法分析目的重采样频率代码示例结果含义日变化分析6小时→24小时日平均.resample(time1D).mean()每天的平均状况周趋势分析6小时→7天周平均.resample(time7D).mean()每周的平均变化极端事件6小时→24小时日最大.resample(time1D).max()每天的最高值累积降水6小时→24小时日总和.resample(time1D).sum()每天的降水量月气候值6小时→月平均.resample(time1M).mean()月平均气候完整的多变量重采样示例def resample_forecast_data(input_path, output_path, freq1D, methodmean): 对伏羲预报结果进行时间重采样 参数 input_path: 输入nc文件路径 output_path: 输出nc文件路径 freq: 重采样频率如 1D日、7D周、1H小时 method: 统计方法mean平均、max最大、min最小、sum求和 # 打开数据集 ds xr.open_dataset(input_path) # 设置时间坐标如果还没有 if time not in ds.coords: start_time pd.Timestamp.now().replace(hour0, minute0, second0) time_coords pd.date_range(startstart_time, periodslen(ds.time), freq6H) ds ds.assign_coords(timetime_coords) # 根据method选择重采样方法 if method mean: ds_resampled ds.resample(timefreq).mean() elif method max: ds_resampled ds.resample(timefreq).max() elif method min: ds_resampled ds.resample(timefreq).min() elif method sum: ds_resampled ds.resample(timefreq).sum() else: raise ValueError(f不支持的method: {method}) # 保存结果 ds_resampled.to_netcdf(output_path) print(f重采样完成文件已保存到: {output_path}) print(f原始时间步数: {len(ds.time)}重采样后: {len(ds_resampled.time)}) return ds_resampled # 使用示例生成日平均预报 resampled_data resample_forecast_data( input_path/root/fuxi2/output.nc, output_path/root/fuxi2/output_daily_mean.nc, freq1D, methodmean )3.5 重采样的注意事项在实际使用中有几个细节需要注意时间坐标必须正确如果数据没有正确的时间坐标重采样会失败。确保你的time维度是datetime类型。边界处理重采样时时间窗口的边界会影响结果。xarray默认使用左闭右开区间即[00:00, 24:00)。缺失值处理如果某个时间窗口内全是缺失值重采样结果也会是缺失值。可以使用.resample().mean(skipnaTrue)跳过缺失值。内存考虑对全球数据进行重采样可能消耗大量内存。如果只需要特定区域可以先裁剪再重采样# 先裁剪中国区域再重采样节省内存 china_region ds.sel(latslice(15, 55), lonslice(70, 140)) china_daily china_region.resample(time1D).mean()4. 方法二滑动平均平滑处理如果你观察过原始的预报时间序列可能会发现曲线有些毛刺——温度忽高忽低气压时升时降。这些短期波动很多时候是预报的不确定性造成的而不是真实的天气变化。滑动平均可以帮助我们平滑这些波动看清真正的趋势。4.1 为什么需要滑动平均想象一下你看股票走势图每天的股价上上下下很难判断是涨是跌。但如果你看5日平均线、20日平均线趋势就清晰多了。天气预报数据也是同样的道理。滑动平均的三个主要作用滤除高频噪声消除6小时尺度上的随机波动凸显趋势让日变化、季节变化等长期趋势更明显提高可视化效果平滑的曲线更容易理解和解释4.2 滑动平均的基本原理滑动平均很简单对于时间序列上的每个点用它前后一段时间内的平均值来代替它自己。举个例子3点滑动平均原始数据 [10, 12, 11, 13, 14, 12, 15] 滑动平均 [ - , 11, 12, 12.67, 13, 13.67, - ]计算过程第二个点(12) (101211)/3 11第三个点(11) (121113)/3 12以此类推...4.3 实现滑动平均平滑让我们用Python实现滑动平均并比较平滑前后的效果import matplotlib.pyplot as plt def sliding_average(data, window_size5): 计算滑动平均 参数 data: 输入数据numpy数组或xarray DataArray window_size: 滑动窗口大小必须是奇数 返回 平滑后的数据两端(window_size//2)个点用NaN填充 if window_size % 2 0: window_size 1 # 确保是奇数 print(f窗口大小调整为奇数: {window_size}) half_window window_size // 2 smoothed np.full_like(data, np.nan) # 用NaN初始化 # 计算滑动平均 for i in range(half_window, len(data) - half_window): smoothed[i] np.mean(data[i-half_window : ihalf_window1]) return smoothed # 提取一段温度时间序列 temp_series ds[T2M].isel(latlat_idx, lonlon_idx).values # 应用不同窗口大小的滑动平均 window_sizes [3, 7, 15] # 对应18小时、42小时、90小时 smoothed_results {} for ws in window_sizes: smoothed_results[ws] sliding_average(temp_series, window_sizews) # 可视化比较 plt.figure(figsize(12, 6)) time_hours np.arange(len(temp_series)) * 6 # 转换为小时 plt.plot(time_hours, temp_series, b-, alpha0.3, label原始数据, linewidth1) colors [r-, g-, m-] for idx, ws in enumerate(window_sizes): plt.plot(time_hours, smoothed_results[ws], colors[idx], labelf{ws}点滑动平均 ({ws*6}小时), linewidth2) plt.xlabel(预报时间 (小时)) plt.ylabel(2米温度 (K)) plt.title(不同滑动窗口大小的平滑效果对比) plt.legend() plt.grid(True, alpha0.3) plt.tight_layout() plt.savefig(/root/fuxi2/smoothing_comparison.png, dpi150) plt.show()运行这段代码你会看到一张对比图清晰地展示不同窗口大小的平滑效果。4.4 xarray内置的滚动方法xarray提供了更便捷的.rolling()方法来实现滑动平均而且可以处理多维数据def smooth_forecast_data(ds, window_size5, dimtime): 使用xarray的rolling方法平滑预报数据 参数 ds: xarray Dataset 或 DataArray window_size: 滑动窗口大小 dim: 沿着哪个维度平滑通常是time 返回 平滑后的数据 # 使用rolling创建滑动窗口然后计算平均 # centerTrue表示窗口居中 ds_smoothed ds.rolling({dim: window_size}, centerTrue).mean() # rolling方法会在两端留下NaN我们可以选择 # 1. 保留NaN反映边界处置信度低 # 2. 用原始值填充简单处理 return ds_smoothed # 对整个数据集进行平滑 smoothed_ds smooth_forecast_data(ds, window_size7) # 比较平滑前后的北京温度 original_temp ds[T2M].isel(latlat_idx, lonlon_idx) smoothed_temp smoothed_ds[T2M].isel(latlat_idx, lonlon_idx) print(平滑前后对比前10个时间步) print(时间 原始温度(K) 平滑温度(K) 差值(K)) print(- * 50) for i in range(10): time_str pd.Timestamp(original_temp.time[i].values).strftime(%m-%d %H:%M) orig original_temp[i].values smooth smoothed_temp[i].values diff smooth - orig if not np.isnan(smooth) else np.nan print(f{time_str} {orig:7.2f} {smooth:7.2f} {diff:6.2f})4.5 滑动平均的高级技巧4.5.1 加权滑动平均有时候我们想让最近的数据权重更大。这就是加权滑动平均def weighted_sliding_average(data, window_size5): 使用三角形加权的滑动平均 weights np.arange(1, window_size//2 2) weights np.concatenate([weights, weights[-2::-1]]) # 创建对称权重 weights weights / weights.sum() # 归一化 half_window window_size // 2 smoothed np.full_like(data, np.nan) for i in range(half_window, len(data) - half_window): smoothed[i] np.sum(data[i-half_window : ihalf_window1] * weights) return smoothed4.5.2 处理缺失值实际数据中常有缺失值我们需要小心处理def sliding_average_with_nan(data, window_size5, min_valid3): 能处理NaN值的滑动平均 参数 min_valid: 窗口内至少需要多少个有效值才计算平均 half_window window_size // 2 smoothed np.full_like(data, np.nan) for i in range(half_window, len(data) - half_window): window_data data[i-half_window : ihalf_window1] valid_mask ~np.isnan(window_data) if np.sum(valid_mask) min_valid: smoothed[i] np.mean(window_data[valid_mask]) return smoothed4.5.3 选择最佳窗口大小窗口大小怎么选太大则过度平滑失去细节太小则平滑效果不足。这里有个经验法则分析目的推荐窗口大小对应时间尺度理由滤除随机噪声3-5点18-30小时保留日变化滤除短时波动分析日变化5-7点30-42小时平滑半天尺度波动分析趋势7-15点42-90小时看清多日变化趋势周尺度分析15-25点90-150小时滤除日变化凸显周趋势你可以用这段代码测试不同窗口大小def find_optimal_window(data, max_window25): 通过计算平滑度指标寻找最佳窗口大小 results [] for ws in range(3, max_window1, 2): # 只测试奇数 smoothed sliding_average(data, window_sizews) # 计算平滑度指标标准差越小越平滑 valid_data smoothed[~np.isnan(smoothed)] if len(valid_data) 0: smoothness 1 / np.std(valid_data) # 平滑度指标 results.append((ws, smoothness)) # 按平滑度排序 results.sort(keylambda x: x[1], reverseTrue) print(窗口大小推荐按平滑度排序) for ws, smoothness in results[:5]: # 显示前5个 print(f 窗口大小 {ws}点 ({ws*6}小时): 平滑度 {smoothness:.3f}) return results[0][0] if results else 5 # 使用示例 optimal_ws find_optimal_window(temp_series, max_window25) print(f\n推荐窗口大小: {optimal_ws}点 ({optimal_ws*6}小时))5. 实战案例完整处理流程现在我们把学到的两种方法结合起来完成一个完整的实战案例分析未来15天中国东部地区的降水趋势。5.1 案例背景与目标假设你是气象分析师需要获取伏羲对未来15天的降水预报将6小时间隔的数据转换为日累积降水平滑日降水序列识别干旱/多雨期生成可视化报告5.2 完整代码实现import xarray as xr import numpy as np import pandas as pd import matplotlib.pyplot as plt from matplotlib import cm import cartopy.crs as ccrs import cartopy.feature as cfeature def process_precipitation_forecast(input_path, regionNone, output_dir.): 完整处理降水预报的流程 参数 input_path: 伏羲输出文件路径 region: 区域范围 [lat_min, lat_max, lon_min, lon_max]默认中国东部 output_dir: 输出目录 # 1. 设置区域默认中国东部 if region is None: region [20, 45, 100, 125] # 中国东部 # 2. 打开数据并裁剪区域 print(步骤1: 加载数据并裁剪区域...) ds xr.open_dataset(input_path) # 设置时间坐标 start_time pd.Timestamp.now().replace(hour0, minute0, second0) time_coords pd.date_range(startstart_time, periodslen(ds.time), freq6H) ds ds.assign_coords(timetime_coords) # 裁剪区域 ds_region ds.sel( latslice(region[0], region[1]), lonslice(region[2], region[3]) ) # 3. 重采样计算日累积降水 print(步骤2: 计算日累积降水...) # TP是6小时累积降水单位是米转换为毫米 precipitation_mm ds_region[TP] * 1000 # 米→毫米 # 重采样为日累积求和 daily_precip precipitation_mm.resample(time1D).sum() # 4. 区域平均计算整个区域的平均日降水 print(步骤3: 计算区域平均...) # 使用面积加权平均考虑纬度变化 lat_weights np.cos(np.deg2rad(daily_precip.lat)) weighted_precip daily_precip.weighted(lat_weights) region_mean weighted_precip.mean(dim[lat, lon]) # 5. 滑动平均平滑 print(步骤4: 应用滑动平均平滑...) window_size 3 # 3天滑动平均 smoothed_mean region_mean.rolling(timewindow_size, centerTrue).mean() # 6. 保存处理结果 print(步骤5: 保存结果...) # 创建包含所有结果的Dataset result_ds xr.Dataset({ daily_precip: daily_precip, region_mean: region_mean, smoothed_mean: smoothed_mean }) output_path f{output_dir}/processed_precipitation.nc result_ds.to_netcdf(output_path) # 7. 生成可视化 print(步骤6: 生成可视化...) generate_precipitation_plots(result_ds, region, output_dir) print(f\n处理完成) print(f原始文件: {input_path}) print(f处理结果: {output_path}) return result_ds def generate_precipitation_plots(result_ds, region, output_dir): 生成降水分析的可视化图表 # 图1区域平均降水时间序列 plt.figure(figsize(12, 8)) # 子图1原始日降水与平滑曲线 plt.subplot(2, 2, 1) time_dates pd.to_datetime(result_ds[region_mean].time.values) plt.bar(time_dates, result_ds[region_mean].values, alpha0.5, label日降水, colorskyblue) plt.plot(time_dates, result_ds[smoothed_mean].values, r-, linewidth2, labelf{3}天滑动平均) plt.xlabel(日期) plt.ylabel(区域平均降水 (mm/天)) plt.title(中国东部区域平均日降水) plt.legend() plt.grid(True, alpha0.3) plt.xticks(rotation45) # 子图2累积降水 plt.subplot(2, 2, 2) cumulative_precip np.cumsum(result_ds[region_mean].values) plt.plot(time_dates, cumulative_precip, g-, linewidth2) plt.fill_between(time_dates, 0, cumulative_precip, alpha0.3, colorgreen) plt.xlabel(日期) plt.ylabel(累积降水 (mm)) plt.title(15天累积降水量) plt.grid(True, alpha0.3) plt.xticks(rotation45) # 子图3降水空间分布第7天示例 plt.subplot(2, 2, 3, projectionccrs.PlateCarree()) day_idx 6 # 第7天 precip_day result_ds[daily_precip].isel(timeday_idx) # 创建地图 ax plt.gca() ax.set_extent([region[2], region[3], region[0], region[1]], crsccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle:) ax.add_feature(cfeature.LAKES, alpha0.5) ax.add_feature(cfeature.RIVERS) # 绘制降水 im precip_day.plot(axax, transformccrs.PlateCarree(), cmapYlGnBu, add_colorbarFalse) plt.colorbar(im, axax, orientationhorizontal, pad0.05, label降水量 (mm)) date_str pd.Timestamp(time_dates[day_idx]).strftime(%Y-%m-%d) plt.title(f{date_str} 降水空间分布) # 子图4降水强度统计 plt.subplot(2, 2, 4) # 计算不同强度降水的天数 precip_values result_ds[region_mean].values bins [0, 1, 5, 10, 20, 50] # 降水强度分级 (mm/天) labels [微量(1), 小雨(1-5), 中雨(5-10), 大雨(10-20), 暴雨(20)] counts, _ np.histogram(precip_values, binsbins) colors [lightgray, lightblue, blue, darkblue, purple] bars plt.bar(range(len(counts)), counts, colorcolors) plt.xlabel(降水强度) plt.ylabel(天数) plt.title(降水强度分布) plt.xticks(range(len(labels)), labels, rotation45) # 在柱子上添加数值 for bar, count in zip(bars, counts): height bar.get_height() plt.text(bar.get_x() bar.get_width()/2., height, f{count}, hacenter, vabottom) plt.tight_layout() plt.savefig(f{output_dir}/precipitation_analysis.png, dpi150, bbox_inchestight) plt.close() # 生成文本报告 generate_text_report(result_ds, output_dir) def generate_text_report(result_ds, output_dir): 生成文本分析报告 mean_precip result_ds[region_mean].values smoothed_precip result_ds[smoothed_mean].values report f伏羲降水预报分析报告 生成时间: {pd.Timestamp.now().strftime(%Y-%m-%d %H:%M:%S)} 分析时段: {pd.Timestamp(result_ds.time[0].values).strftime(%Y-%m-%d)} 至 {pd.Timestamp(result_ds.time[-1].values).strftime(%Y-%m-%d)} 数据来源: FuXi 15天全球天气预报系统 【统计摘要】 - 平均日降水: {np.nanmean(mean_precip):.2f} mm/天 - 最大日降水: {np.nanmax(mean_precip):.2f} mm/天 (第{np.nanargmax(mean_precip)1}天) - 最小日降水: {np.nanmin(mean_precip):.2f} mm/天 (第{np.nanargmin(mean_precip)1}天) - 总累积降水: {np.nansum(mean_precip):.2f} mm - 降水天数(1mm): {np.sum(mean_precip 1)} 天 - 干旱天数(0.5mm): {np.sum(mean_precip 0.5)} 天 【趋势分析】 # 添加趋势判断 if len(smoothed_precip) 5: first_half np.nanmean(smoothed_precip[:5]) second_half np.nanmean(smoothed_precip[-5:]) if second_half first_half * 1.2: report - 降水呈增加趋势后期可能更湿润\n elif second_half first_half * 0.8: report - 降水呈减少趋势后期可能更干燥\n else: report - 降水趋势相对平稳\n # 识别关键特征日 report \n【关键特征日】\n # 找峰值和谷值 valid_indices ~np.isnan(smoothed_precip) if np.sum(valid_indices) 0: smooth_valid smoothed_precip[valid_indices] time_valid result_ds.time.values[valid_indices] # 找局部最大值峰值 from scipy.signal import find_peaks peaks, _ find_peaks(smooth_valid, heightnp.nanmean(smooth_valid)) for idx in peaks[:3]: # 显示前3个峰值 date_str pd.Timestamp(time_valid[idx]).strftime(%Y-%m-%d) report f- {date_str}: 降水峰值预计{smooth_valid[idx]:.1f} mm/天\n # 找局部最小值谷值 valleys, _ find_peaks(-smooth_valid) for idx in valleys[:3]: # 显示前3个谷值 date_str pd.Timestamp(time_valid[idx]).strftime(%Y-%m-%d) report f- {date_str}: 降水谷值预计{smooth_valid[idx]:.1f} mm/天\n report \n【使用建议】\n report 1. 本分析基于伏羲AI天气预报仅供参考\n report 2. 滑动平均帮助识别趋势但会平滑短期变化\n report 3. 建议结合其他预报资料进行综合判断\n # 保存报告 report_path f{output_dir}/precipitation_report.txt with open(report_path, w, encodingutf-8) as f: f.write(report) print(f分析报告已保存: {report_path}) # 运行完整流程 if __name__ __main__: result process_precipitation_forecast( input_path/root/fuxi2/output.nc, output_dir/root/fuxi2/results )5.3 案例输出解读运行上面的代码后你会得到processed_precipitation.nc包含处理后的NetCDF数据precipitation_analysis.png四合一可视化图表precipitation_report.txt文本分析报告报告示例内容伏羲降水预报分析报告 生成时间: 2024-01-15 14:30:00 分析时段: 2024-01-16 至 2024-01-30 【统计摘要】 - 平均日降水: 3.2 mm/天 - 最大日降水: 8.5 mm/天 (第7天) - 最小日降水: 0.2 mm/天 (第3天) - 总累积降水: 48.3 mm - 降水天数(1mm): 9天 - 干旱天数(0.5mm): 4天 【趋势分析】 - 降水呈增加趋势后期可能更湿润 【关键特征日】 - 2024-01-22: 降水峰值预计7.8 mm/天 - 2024-01-25: 降水峰值预计6.5 mm/天 - 2024-01-18: 降水谷值预计0.8 mm/天 【使用建议】 1. 本分析基于伏羲AI天气预报仅供参考 2. 滑动平均帮助识别趋势但会平滑短期变化 3. 建议结合其他预报资料进行综合判断6. 总结与最佳实践通过本文的学习你已经掌握了伏羲天气预报结果后处理的两大核心技能时间维度重采样和滑动平均平滑。让我们回顾一下关键要点6.1 技术要点回顾重采样改变时间频率将6小时间隔的预报转换为日、周等更有意义的统计量使用.resample()方法选择合适频率1D、7D等根据需求选择统计方法.mean()、.sum()、.max()等滑动平均平滑时间序列滤除随机波动凸显真实趋势使用.rolling()方法设置合适的窗口大小窗口大小通常选3-7点18-42小时平衡细节与平滑度考虑使用加权平均或处理缺失值的变体组合使用效果更佳先重采样降低数据频率再滑动平均平滑曲线是常见的最佳实践6.2 实用建议根据不同的应用场景这里有一些实用建议应用场景推荐处理流程理由气候趋势分析重采样为月平均 → 12月滑动平均滤除季节内波动看清年际变化农业气象服务重采样为日累积 → 7天滑动平均提供周尺度水分供需分析灾害预警保持6小时原始频率 → 3点滑动平均保留短时强降水信号仅滤除噪声能源需求预测重采样为日平均温度 → 5点滑动平均平滑温度曲线用于负荷预测航空天气预报保持原始频率必要时3点滑动平均需要高时间分辨率仅轻微平滑6.3 常见问题与解决Q: 处理后的数据还能用伏羲继续预报吗A: 不能。重采样和平滑会改变数据的时间结构和统计特性不适合作为下一轮预报的初始场。这些处理主要用于分析和可视化。Q: 滑动平均会导致预报滞后吗A: 使用centerTrue的滑动平均不会导致滞后因为它是用前后数据计算当前点的平均值。但边界处开头和结尾由于数据不足平滑效果会打折扣。Q: 如何处理全球数据的内存问题A: 三个策略先裁剪感兴趣的区域再处理使用Dask进行并行处理分变量处理保存中间结果Q: 有现成的工具可以直接用吗A: 你可以将本文的代码封装成函数库方便重复使用。这里提供一个简单的封装示例# fuxi_postprocess.py import xarray as xr import numpy as np import pandas as pd class FuXiPostProcessor: 伏羲预报后处理工具类 def __init__(self, data_path): self.ds xr.open_dataset(data_path) self._set_time_coords() def _set_time_coords(self): 设置时间坐标 if time not in self.ds.coords: start_time pd.Timestamp.now().replace(hour0, minute0, second0) time_coords pd.date_range(startstart_time, periodslen(self.ds.time), freq6H) self.ds self.ds.assign_coords(timetime_coords) def resample(self, freq1D, methodmean, variablesNone): 重采样数据 if variables is None: variables list(self.ds.data_vars) if method mean: return self.ds[variables].resample(timefreq).mean() elif method sum: return self.ds[variables].resample(timefreq).sum() # ... 其他方法 def smooth(self, window_size5, dimtime): 滑动平均平滑 return self.ds.rolling({dim: window_size}, centerTrue).mean() def process_region(self, region, variable, freq1D, window_size3): 完整的区域处理流程 # 裁剪区域 ds_region self.ds.sel( latslice(region[0], region[1]), lonslice(region[2], region[3]) ) # 重采样 ds_resampled ds_region.resample(timefreq).mean() # 区域平均 lat_weights np.cos(np.deg2rad(ds_resampled.lat)) weighted_data ds_resampled[variable].weighted(lat_weights) region_mean weighted_data.mean(dim[lat, lon]) # 平滑 smoothed region_mean.rolling(timewindow_size, centerTrue).mean() return { resampled: ds_resampled, region_mean: region_mean, smoothed: smoothed } # 使用示例 processor FuXiPostProcessor(/root/fuxi2/output.nc) result processor.process_region( region[20, 45, 100, 125], variableT2M, freq1D, window_size3 )6.4 下一步学习建议掌握了基础的后处理技能后你可以进一步探索多模型集成将伏羲预报与其他模型结果结合提高预报可靠性统计后处理使用MOSModel Output Statistics等技术校准预报偏差可视化进阶学习使用Cartopy、Metpy等专业气象可视化库自动化流程将处理流程封装成定时任务实现自动化预报产品生成天气预报的最终价值在于应用。通过对伏羲预报结果的合理后处理你可以将原始的数值输出转化为决策支持信息真正发挥AI天气预报的实用价值。获取更多AI镜像想探索更多AI镜像和应用场景访问 CSDN星图镜像广场提供丰富的预置镜像覆盖大模型推理、图像生成、视频生成、模型微调等多个领域支持一键部署。
