利用seawind数据对WRF模式和GSI同化的风结果进行评估（python）

---

前言

评估是模式预报和资料同化的重要一步，这里主要介绍一下用一般的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模式方案的影响，有待进一步验证。