文章目录
前言
评估是模式预报和资料同化的重要一步,这里主要介绍一下用一般的bias和RMSE来检验模拟的可靠性,虽然结果不是很乐观,但不影响评估方法。
一、实验介绍
接下来简单介绍一下实验所用的数据和方法
1.数据介绍
1.1 gfs(Global Forecast System)数据(用于驱动WRF模式)
GFS数据是由全球预报系统(GFS)模式gfs官网生成的全球数值天气预报数据。该数据覆盖全球,从地表到高层大气,提供多种气象变量的信息,如温度、风速、湿度、降水量、气压等。GFS数据以较高的空间分辨率(通常为0.25度或0.5度)和时间分辨率(1小时或3小时)提供,适用于天气预报、气候分析、环境监测等领域。GFS数据每天更新四次(00Z、06Z、12Z、18Z),并提供长达16天的预报信息。
1.2 gdas数据
GDAS (Global Data Assimilation System,全球数据同化系统), GDAS数据是由美国国家环境预报中心(NCEP)开发的一个综合全球观测数据的系统生成的数据产品。GDAS通过将来自卫星、地面站、气象气球、飞机和船舶等多种观测源的数据整合到数值天气预报模式中,以生成一个最佳的全球大气状态分析。
GDAS数据提供了全球大气的三维结构,包括温度、湿度、风速、气压等变量,通常以较高的空间分辨率(例如1度或0.25度)和时间分辨率(6小时)发布。GDAS的分析结果是GFS等数值预报模式的初始条件,因此在提高预报精度方面起着关键作用。GDAS数据广泛应用于天气预报、气候研究、环境监测等领域。
1.3 seawinds数据
seawinds是指由SeaWinds仪器收集的海洋表面风场数据。SeaWinds是一种微波散射计,搭载在NASA的QuikSCAT卫星和NASA/JAXA的ADEOS-II卫星上。通过发射微波信号并测量其返回的散射强度,SeaWinds能够精确地测量海洋表面的风速和风向。
Seawinds数据以高分辨率提供全球海洋表面风场信息,通常网格分辨率为0.25度,时间分辨率为每日多次。由于这些数据涵盖全球范围内的海洋表面,特别是在风暴、飓风等极端天气事件的监测和预报中起到了至关重要的作用。Seawinds数据广泛应用于气象学、海洋学、气候研究以及天气预报的改进,帮助科学家更好地理解和模拟大气与海洋的相互作用。
2.模式介绍
2.1WRF
WRF (Weather Research and Forecasting)WRF模式是一种广泛应用的数值天气预报和大气模拟工具。它由多个美国机构联合开发,适用于从区域到全球的气象研究和预报。WRF具有灵活的模块化设计,支持多种物理参数化方案,能够模拟大气中的各种复杂过程,常用于天气预报、气候研究、环境影响评估等领域。其高分辨率和可扩展性使得WRF成为科研和操作预报中的重要工具。
2.1.1模式设置
2.2 GSI同化系统
GSI(Gridpoint Statistical Interpolation,网格点统计插值)是一个三维变分(3DVAR)数据同化系统,广泛用于数值天气预报模型中。由美国国家环境预报中心(NCEP)开发,GSI将各种观测数据(如卫星、雷达、气象站等)与数值预报模型结合,通过最优插值方法生成更接近真实大气状态的初始条件。这些改进的初始条件用于驱动数值天气预报模型,从而提高预报的准确性。GSI系统灵活性强,适用于全球和区域模式中的多种物理和观测系统,是气象预报和研究的重要工具。
GSI与很多海洋数据同化系统不同,其物理约束更为严格。
同化方案就不仔细介绍了。
二、模式评估
1.插值
利用seawinds数据作为真值检验模式预报效果时,需要将两者数据插值到同一时空分辨率,
1.1 seawinds数据的提取
seawinds数据时每天一个nc文件,所以评估多天不同时间段需要将seawinds数据提取出来,并且存入新的nc文件。
程序如下
import numpy as np
from netCDF4 import Dataset
# 定义输入文件路径
seawind_file_29 = "/seawind/20240628-0630/NBSv02_wind_6hourly_20240629.nc"
seawind_file_30 = "/seawind/20240628-0630/NBSv02_wind_6hourly_20240630.nc"
# 定义目标经纬度范围
lon_min, lon_max = 121.5, 125
lat_min, lat_max = 33, 36
# 打开 Seawind NetCDF 文件
f29 = Dataset(seawind_file_29)
f30 = Dataset(seawind_file_30)
# 从 f29 文件中获取经纬度坐标
lons_f29 = f29.variables['lon'][:]
lats_f29 = f29.variables['lat'][:]
# 从 f30 文件中获取经纬度坐标
lons_f30 = f30.variables['lon'][:]
lats_f30 = f30.variables['lat'][:]
# 获取目标经纬度范围内的数据的索引
def get_indices(lons, lats, lon_min, lon_max, lat_min, lat_max):
lon_indices = np.where((lons >= lon_min) & (lons <= lon_max))[0]
lat_indices = np.where((lats >= lat_min) & (lats <= lat_max))[0]
return lon_indices, lat_indices
# 从 f29 文件中获取目标范围内的数据的索引
lon_indices_f29, lat_indices_f29 = get_indices(lons_f29, lats_f29, lon_min, lon_max, lat_min, lat_max)
# 从 f30 文件中获取目标范围内的数据的索引
lon_indices_f30, lat_indices_f30 = get_indices(lons_f30, lats_f30, lon_min, lon_max, lat_min, lat_max)
# 目标时间点索引
time_indices_29 = [0, 1, 2, 3] # 代表29日的四个时间点
time_indices_30 = [0, 1, 2] # 代表30日的0、6、12时
time_indices = time_indices_29 + time_indices_30
# 创建目标文件路径
output_file_path = "/DatadiskExt/xychen/seawind/20240628-0630/seawind_extracted.nc"
# 创建目标文件
with Dataset(output_file_path, 'w', format='NETCDF4') as ncfile:
# 创建维度
ncfile.createDimension('time', len(time_indices))
ncfile.createDimension('lat', len(lat_indices_f29))
ncfile.createDimension('lon', len(lon_indices_f29))
# 创建变量
times = ncfile.createVariable('time', np.float64, ('time',)) # 使用float64类型存储时间索引
lats = ncfile.createVariable('lat', np.float32, ('lat',))
lons = ncfile.createVariable('lon', np.float32, ('lon',))
u_wind_var = ncfile.createVariable('u_wind', np.float32, ('time', 'lat', 'lon'))
v_wind_var = ncfile.createVariable('v_wind', np.float32, ('time', 'lat', 'lon'))
# 写入数据
lats[:] = lats_f29[lat_indices_f29] # 提取目标纬度范围内的数据
lons[:] = lons_f29[lon_indices_f29] # 提取目标经度范围内的数据
times[:] = [0, 6, 12, 18, 0, 6, 12] # 时间索引
# 提取29日四个时刻的数据
for idx, t in enumerate(time_indices_29):
u_wind = f29.variables['u_wind'][t, 0, :, :]
v_wind = f29.variables['v_wind'][t, 0, :, :]
u_wind_var[idx, :, :] = u_wind[lat_indices_f29, :][:, lon_indices_f29]
v_wind_var[idx, :, :] = v_wind[lat_indices_f29, :][:, lon_indices_f29]
# 提取30日0、6、12时的数据
for idx, t in enumerate(time_indices_30):
u_wind = f30.variables['u_wind'][t, 0, :, :]
v_wind = f30.variables['v_wind'][t, 0, :, :]
u_wind_var[len(time_indices_29) + idx, :, :] = u_wind[lat_indices_f30, :][:, lon_indices_f30]
v_wind_var[len(time_indices_29) + idx, :, :] = v_wind[lat_indices_f30, :][:, lon_indices_f30]
print(f"seawind_extracted.nc 文件已保存到 {output_file_path}")
1.2 WRF输出数据的提取和插值
WRF输出的数据的水平分布是在其自己的质量网格点,往往都不是整数,不像seawind数据那么工整,所以需要对其插值后才能进行评估。
import numpy as np
from netCDF4 import Dataset
from scipy.interpolate import griddata
# 定义输入和输出文件路径
f1_path = "/2024062812/bkg/wrfout_d01_2024-06-28_12:00:00"#初始模拟
f2_path = "/wrfout/2024062812/wrfout_before_update_bc" #assimilation without updata_bc
f3_path = "/wrfout/2024062812/wrfout_after_update_bc" #assimilation with updata_bc
output_dir = "/DatadiskExt/xychen/wrfout/2024062812/"
# 定义目标经纬度网格点
lon_target = np.arange(121.5, 125.25, 0.25)
lat_target = np.arange(33, 36.25, 0.25)
target_lons, target_lats = np.meshgrid(lon_target, lat_target)
# 打开 NetCDF 文件
f1 = Dataset(f1_path)
f2 = Dataset(f2_path)
f3 = Dataset(f3_path)
# 获取时间索引
time_indices = [0, 6, 12, 18, 24, 30, 36, 42, 48]
# 创建插值结果的空列表
u_f1_interp_all = []
v_f1_interp_all = []
u_f2_interp_all = []
v_f2_interp_all = []
u_f3_interp_all = []
v_f3_interp_all = []
# 对每个时间点进行插值
for idx in time_indices:
# 从 f1 文件中获取时间索引处的风速数据
u_f1 = f1.variables['U10'][idx, :, :]
v_f1 = f1.variables['V10'][idx, :, :]
# 从 f1 文件中获取经纬度坐标
lats_f1 = f1.variables['XLAT'][0, :, :]
lons_f1 = f1.variables['XLONG'][0, :, :]
# 从 f2 文件中获取时间索引处的风速数据
u_f2 = f2.variables['U10'][idx, :, :]
v_f2 = f2.variables['V10'][idx, :, :]
# 从 f2 文件中获取经纬度坐标
lats_f2 = f2.variables['XLAT'][idx, :, :]
lons_f2 = f2.variables['XLONG'][idx, :, :]
# 从 f3 文件中获取时间索引处的风速数据
u_f3 = f3.variables['U10'][idx, :, :]
v_f3 = f3.variables['V10'][idx, :, :]
# 从 f3 文件中获取经纬度坐标
lats_f3 = f3.variables['XLAT'][idx, :, :]
lons_f3 = f3.variables['XLONG'][idx, :, :]
# 将数据转换为一维坐标和风速数据
lons_f1_flat = lons_f1.flatten()
lats_f1_flat = lats_f1.flatten()
u_f1_flat = u_f1.flatten()
v_f1_flat = v_f1.flatten()
lons_f2_flat = lons_f2.flatten()
lats_f2_flat = lats_f2.flatten()
u_f2_flat = u_f2.flatten()
v_f2_flat = v_f2.flatten()
lons_f3_flat = lons_f3.flatten()
lats_f3_flat = lats_f3.flatten()
u_f3_flat = u_f3.flatten()
v_f3_flat = v_f3.flatten()
# 执行插值到目标网格点
u_f1_interp = griddata((lons_f1_flat, lats_f1_flat), u_f1_flat, (target_lons, target_lats), method='linear')
v_f1_interp = griddata((lons_f1_flat, lats_f1_flat), v_f1_flat, (target_lons, target_lats), method='linear')
u_f2_interp = griddata((lons_f2_flat, lats_f2_flat), u_f2_flat, (target_lons, target_lats), method='linear')
v_f2_interp = griddata((lons_f2_flat, lats_f2_flat), v_f2_flat, (target_lons, target_lats), method='linear')
u_f3_interp = griddata((lons_f3_flat, lats_f3_flat), u_f3_flat, (target_lons, target_lats), method='linear')
v_f3_interp = griddata((lons_f3_flat, lats_f3_flat), v_f3_flat, (target_lons, target_lats), method='linear')
# 将插值结果重塑为网格
u_f1_interp_reshaped = u_f1_interp.reshape((len(lat_target), len(lon_target)))
v_f1_interp_reshaped = v_f1_interp.reshape((len(lat_target), len(lon_target)))
u_f2_interp_reshaped = u_f2_interp.reshape((len(lat_target), len(lon_target)))
v_f2_interp_reshaped = v_f2_interp.reshape((len(lat_target), len(lon_target)))
u_f3_interp_reshaped = u_f3_interp.reshape((len(lat_target), len(lon_target)))
v_f3_interp_reshaped = v_f3_interp.reshape((len(lat_target), len(lon_target)))
# 添加到插值结果列表中
u_f1_interp_all.append(u_f1_interp_reshaped)
v_f1_interp_all.append(v_f1_interp_reshaped)
u_f2_interp_all.append(u_f2_interp_reshaped)
v_f2_interp_all.append(v_f2_interp_reshaped)
u_f3_interp_all.append(u_f3_interp_reshaped)
v_f3_interp_all.append(v_f3_interp_reshaped)
# 保存 before_gsi_interpolate.nc 文件
before_gsi_path = output_dir + "before_gsi_interpolate.nc"
with Dataset(before_gsi_path, 'w', format='NETCDF4') as ncfile:
# 创建维度
ncfile.createDimension('time', len(time_indices))
ncfile.createDimension('lat', len(lat_target))
ncfile.createDimension('lon', len(lon_target))
# 创建变量
times = ncfile.createVariable('time', np.float64, ('time',)) # 使用float64类型存储时间索引
lats = ncfile.createVariable('lat', np.float32, ('lat',))
lons = ncfile.createVariable('lon', np.float32, ('lon',))
u_f1_interp_var = ncfile.createVariable('u', np.float32, ('time', 'lat', 'lon'))
v_f1_interp_var = ncfile.createVariable('v', np.float32, ('time', 'lat', 'lon'))
# 写入数据
times[:] = np.array(time_indices) # 直接存储时间索引
lats[:] = lat_target
lons[:] = lon_target
u_f1_interp_var[:, :, :] = np.array(u_f1_interp_all)
v_f1_interp_var[:, :, :] = np.array(v_f1_interp_all)
print(f"before_gsi_interpolate.nc 文件已保存到 {before_gsi_path}")
# 保存 after_gsi_interpolate.nc 文件
after_gsi_path = output_dir + "after_gsi_interpolate.nc"
with Dataset(after_gsi_path, 'w', format='NETCDF4') as ncfile:
# 创建维度
ncfile.createDimension('time', len(time_indices))
ncfile.createDimension('lat', len(lat_target))
ncfile.createDimension('lon', len(lon_target))
# 创建变量
times = ncfile.createVariable('time', np.float64, ('time',)) # 使用float64类型存储时间索引
lats = ncfile.createVariable('lat', np.float32, ('lat',))
lons = ncfile.createVariable('lon', np.float32, ('lon',))
u_f2_interp_var = ncfile.createVariable('u', np.float32, ('time', 'lat', 'lon'))
v_f2_interp_var = ncfile.createVariable('v', np.float32, ('time', 'lat', 'lon'))
# 写入数据
times[:] = np.array(time_indices) # 直接存储时间索引
lats[:] = lat_target
lons[:] = lon_target
u_f2_interp_var[:, :, :] = np.array(u_f2_interp_all)
v_f2_interp_var[:, :, :] = np.array(v_f2_interp_all)
print(f"after_gsi_interpolate.nc 文件已保存到 {after_gsi_path}")
# 保存 after_update_bc_interpolate.nc 文件
after_update_bc_path = output_dir + "after_update_bc_interpolate.nc"
with Dataset(after_update_bc_path, 'w', format='NETCDF4') as ncfile:
# 创建维度
ncfile.createDimension('time', len(time_indices))
ncfile.createDimension('lat', len(lat_target))
ncfile.createDimension('lon', len(lon_target))
# 创建变量
times = ncfile.createVariable('time', np.float64, ('time',)) # 使用float64类型存储时间索引
lats = ncfile.createVariable('lat', np.float32, ('lat',))
lons = ncfile.createVariable('lon', np.float32, ('lon',))
u_f3_interp_var = ncfile.createVariable('u', np.float32, ('time', 'lat', 'lon'))
v_f3_interp_var = ncfile.createVariable('v', np.float32, ('time', 'lat', 'lon'))
# 写入数据
times[:] = np.array(time_indices) # 直接存储时间索引
lats[:] = lat_target
lons[:] = lon_target
u_f3_interp_var[:, :, :] = np.array(u_f3_interp_all)
v_f3_interp_var[:, :, :] = np.array(v_f3_interp_all)
print(f"after_update_bc_interpolate.nc 文件已保存到 {after_update_bc_path}")
结果输出了3个用于对比且与seawind提取出来的数据分辨率一致的nc文件
理论上是只有两个数据需要进行对比评估,但再数据同化后多了一个updata_bc(更新边界条件)操作,为了观察更新边界条件的作用,所以放在一起看。
2. 评估
再有了相同分辨率的数据之后,就可以做评估了。
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
import os
# 打开 NetCDF 文件
f_before = Dataset("/DatadiskExt/xychen/wrfout/2024062812/before_gsi_interpolate.nc", "r")
f_after = Dataset("/DatadiskExt/xychen/wrfout/2024062812/after_gsi_interpolate.nc", "r")
f_seawind = Dataset("/DatadiskExt/xychen/seawind/20240628-0630/seawind_extracted.nc", "r")
f_new = Dataset("/DatadiskExt/xychen/wrfout/2024062812/after_update_bc_interpolate.nc", "r") # 新文件路径
# 提取 `before_gsi` 和 `after_gsi` 的 `u` 和 `v` 风速数据
before_gsi_u = f_before.variables['u'][:] # (9, 12, 16)
before_gsi_v = f_before.variables['v'][:] # (9, 12, 16)
after_gsi_u = f_after.variables['u'][:] # (9, 12, 16)
after_gsi_v = f_after.variables['v'][:] # (9, 12, 16)
# 提取 `Seawind` 的 `u_wind` 和 `v_wind` 风速数据
seawind_u = f_seawind.variables['u_wind'][:] # (7, 12, 16)
seawind_v = f_seawind.variables['v_wind'][:] # (7, 12, 16)
# 提取新文件的 `u` 和 `v` 风速数据
new_u = f_new.variables['u'][:] # 根据实际新文件变量名调整
new_v = f_new.variables['v'][:] # 根据实际新文件变量名调整
# 计算 `Seawind` 的风速
seawind_speed = np.sqrt(seawind_u**2 + seawind_v**2)
# 关闭 NetCDF 文件
f_before.close()
f_after.close()
f_seawind.close()
f_new.close()
# 提取 `before_gsi` 和 `after_gsi` 后7个时间点的数据
time_indices = range(2, 9) # 对应最后7个时间点,索引为2到8
before_gsi_u = before_gsi_u[time_indices]
before_gsi_v = before_gsi_v[time_indices]
after_gsi_u = after_gsi_u[time_indices]
after_gsi_v = after_gsi_v[time_indices]
new_u=new_u[time_indices]
new_v=new_v[time_indices]
# 计算风速
before_gsi_speed = np.sqrt(before_gsi_u**2 + before_gsi_v**2)
after_gsi_speed = np.sqrt(after_gsi_u**2 + after_gsi_v**2)
new_speed = np.sqrt(new_u**2 + new_v**2)
# 计算 Mean Bias 和 RMSE
mean_bias_before = np.mean(np.abs(before_gsi_speed - seawind_speed), axis=(1, 2))
rmse_before = np.sqrt(np.mean((before_gsi_speed - seawind_speed) ** 2, axis=(1, 2)))
mean_bias_after = np.mean(np.abs(after_gsi_speed - seawind_speed), axis=(1, 2))
rmse_after = np.sqrt(np.mean((after_gsi_speed - seawind_speed) ** 2, axis=(1, 2)))
mean_bias_new = np.mean(np.abs(new_speed - seawind_speed), axis=(1, 2))
rmse_new = np.sqrt(np.mean((new_speed - seawind_speed) ** 2, axis=(1, 2)))
# 计算 `Mean Bias` 和 `RMSE` 的平均值
mean_bias_before_avg = np.mean(mean_bias_before)
rmse_before_avg = np.mean(rmse_before)
mean_bias_after_avg = np.mean(mean_bias_after)
rmse_after_avg = np.mean(rmse_after)
mean_bias_new_avg = np.mean(mean_bias_new)
rmse_new_avg = np.mean(rmse_new)
# 计算平均相对误差
relative_error_before = np.mean(np.abs(before_gsi_speed - seawind_speed) / np.abs(seawind_speed), axis=(1, 2))
relative_error_after = np.mean(np.abs(after_gsi_speed - seawind_speed) / np.abs(seawind_speed), axis=(1, 2))
relative_error_new = np.mean(np.abs(new_speed - seawind_speed) / np.abs(seawind_speed), axis=(1, 2))
# 计算 `平均相对误差` 的平均值
relative_error_before_avg = np.mean(relative_error_before)
relative_error_after_avg = np.mean(relative_error_after)
relative_error_new_avg = np.mean(relative_error_new)
# 计算预报时间(小时)
forecast_hours = np.array([12, 18, 24, 30, 36, 42, 48]) # 横坐标时预报时间
# 创建输出路径
output_path = '/Users/xychen/pythonpainting/bias_rmse_gsi/'
if not os.path.exists(output_path):
os.makedirs(output_path) # 如果路径不存在,则创建路径
# 绘制 Mean Bias, RMSE 和 平均相对误差
plt.figure(figsize=(18, 6))
# 绘制 Mean Bias
plt.subplot(1, 3, 1)
plt.plot(forecast_hours, mean_bias_before, label=f'Before GSI (Mean Bias: {mean_bias_before_avg:.2f})', marker='o', linestyle='-', color='blue')
plt.plot(forecast_hours, mean_bias_after, label=f'After GSI (Mean Bias: {mean_bias_after_avg:.2f})', marker='o', linestyle='-', color='red')
plt.plot(forecast_hours, mean_bias_new, label=f'After Update (Mean Bias: {mean_bias_new_avg:.2f})', marker='o', linestyle='-', color='yellow') # 新文件用黄色表示
plt.xlabel('Forecast time')
plt.ylabel('Mean Bias')
plt.title('Mean Bias ')
plt.legend()
plt.grid(True)
# 绘制 RMSE
plt.subplot(1, 3, 2)
plt.plot(forecast_hours, rmse_before, label=f'Before GSI (RMSE: {rmse_before_avg:.2f})', marker='o', linestyle='-', color='blue')
plt.plot(forecast_hours, rmse_after, label=f'After GSI (RMSE: {rmse_after_avg:.2f})', marker='o', linestyle='-', color='red')
plt.plot(forecast_hours, rmse_new, label=f'After Update (RMSE: {rmse_new_avg:.2f})', marker='o', linestyle='-', color='yellow') # 新文件用黄色表示
plt.xlabel('Forecast time')
plt.ylabel('RMSE')
plt.title('RMSE ')
plt.legend()
plt.grid(True)
# 绘制 平均相对误差
plt.subplot(1, 3, 3)
plt.plot(forecast_hours, relative_error_before, label=f'Before GSI (Rel. Error: {relative_error_before_avg:.2f})', marker='o', linestyle='-', color='blue')
plt.plot(forecast_hours, relative_error_after, label=f'After GSI (Rel. Error: {relative_error_after_avg:.2f})', marker='o', linestyle='-', color='red')
plt.plot(forecast_hours, relative_error_new, label=f'After Update (Rel. Error: {relative_error_new_avg:.2f})', marker='o', linestyle='-', color='yellow') # 新文件用黄色表示
plt.xlabel('Forecast time')
plt.ylabel('Relative Error')
plt.title('Mean relative Bias')
plt.legend()
plt.grid(True)
plt.tight_layout()
# 保存图形到指定路径
plt.savefig(output_path + 'bias_rmse_relative_error_comparison.png')
plt.show()
print('图已经画好了,大王')
三张分别是每个时间点的平均误差、均方根误差和相对平均误差的变化
总结
其实结果显示同化后的误差和均方根误差都比没有同化的好,可能受到风速、插值方法、同化方案、WRF模式方案的影响,有待进一步验证。