评估方法

#!/usr/bin/env python
from __future__ import print_function 
from netCDF4 import Dataset 
import numpy as np 
from scipy import stats 
import sys, getopt 
import argparse 
from WRF import getvar, interplevel 
  
parser = argparse.ArgumentParser(description='xxxxx.') 
#parser.add_argument('--key', '-k') 
parser.add_argument('file1') 
parser.add_argument('file2') 
  
args = parser.parse_args() 
  
try: 
    f1 = Dataset(args.file1,'r') 
    f2 = Dataset(args.file2,'r') 
except: 
    print("file open failed") 
    exit(2) 
# 计算的要素 
elements = ('T2', 'RAINNC') 
indexnumber = len(f1.variables[elements[0]]) 
  
for index in range(indexnumber): 
    print ('%-24s %16s %16s %16s %16s' %('Var', 'Pearson', 'ME', 'MAE', 'RMSE')) 
    for var in elements: 
        var_split = var.split('-') 
        if len(var_split) > 1: 
            if var_split[0]=='wspd_wdir': 
                var0_f1 = getvar(f1, var_split[0])[0,:] 
                var0_f2 = getvar(f2, var_split[0])[0,:] 
  
            else: 
                var0_f1 = getvar(f1, var_split[0]) 
                var0_f2 = getvar(f2, var_split[0]) 
  
            var1_f1 = getvar(f1, var_split[1]) 
            f1_para = interplevel(var0_f1, var1_f1, int(var_split[2]))[index] 
            var1_f2 = getvar(f2, var_split[1]) 
            f2_para = interplevel(var0_f2, var1_f2, int(var_split[2]))[index] 
        else: 
            f1_para = f1.variables[var][index] 
            f2_para = f2.variables[var][index] 
        # 计算空间相关系数SCC 
        f1_mean = np.mean(f1_para) 
        f2_mean = np.mean(f2_para) 
  
        divisor = np.sum((f1_para - f1_mean) * (f2_para - f2_mean)) 
        f1_dividend = np.sqrt(np.sum((f1_para - f1_mean) ** 2.)) 
        f2_dividend = np.sqrt(np.sum((f2_para - f2_mean) ** 2.)) 
        
  
        if f1_dividend * f2_dividend == 0: 
            pearson = 'nan' 
        else: 
            pearson = divisor / (f1_dividend * f2_dividend) 
        # 计算平均值误差       
        mean = np.mean(f1_para - f2_para) 
        # 计算平均绝对误差 
        mae = np.mean(abs(f1_para - f2_para)) 
        # 计算均方根误差        
        rmse = np.sqrt(np.mean((f1_para - f2_para)**2)) 
         
        print ('%-24s %10.6f %10.6f %10.6f %10.6f'  %(var, pearson, mean, mae, rmse))

WRF应用自带比较工具。WRF编译之后生成在external/io_netcdf目录，可以计算两个结果文件各要素的均方根（RMS）和均方根误差（RMSE），输出不一致的要素。如两个结果文件中全部RMSE值为0，则无输出。如下图所示：

md5sum是一种计算机程序，它通过对文件内容进行哈希编码生成MD5校验码，用于计算与校验RFC 1321所描述的128位MD5哈希值。MD5值可以作为一个文件的数字指纹使用，使用md5sum命令可以查看不同平台上同一算例、同一时刻的结果文件。MD5值相同，则说明结果文件完全一致。

NCL（NCAR Command Language）是由美国大气研究中心（NCAR）推出的一款用于科学数据计算的免费可视化软件。它有着非常强大的文件输入和输出功能，可读写netCDF-3、netCDF-4 classic、HDF4、binary、ASCII等数据。WRF输出的ncl文件可使用wrf的ncl插件进行对比出图。

NCL对比出图脚本：

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" 
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" 
 
begin 
f=addfile("kunpeng/wrfout_d01_2020-07-08_00:00:00","r") 
f1=addfile("x86/wrfout_d01_2020-07-08_00:00:00","r") 
 
res = True                           ; Set up some basic plot resources 
res@MainTitle = "METGRID FILES" 
res@Footer = False 
res@cnFillOn=True 
res@gsnSpreadColors = True 
opts = res 
pltres=True 
mpres=True 
 
wks = gsn_open_wks("pdf","wrf") 
 
gsn_define_colormap(wks,"NCV_bright") 
 
ter = wrf_user_getvar(f,"RAINNC",0) 
contour = wrf_contour(f,wks,ter,res) 
plot = wrf_map_overlays(f,wks,(/contour/),pltres,mpres) 
ter1 = wrf_user_getvar(f1,"RAINNC",0) 
contour = wrf_contour(f1,wks,ter1,res) 
plot = wrf_map_overlays(f1,wks,(/contour/),pltres,mpres) 
print(ter-ter1) 
 
contour = wrf_contour(f,wks,ter-ter1,res) 
plot = wrf_map_overlays(f,wks,(/contour/),pltres,mpres) 
end

图1 NCL对比出图实例

GrADS（Grid Analysis and Display System）是当今气象界广泛使用的一种数据处理和显示软件系统。该软件系统通过其集成环境对气象数据进行读取、加工、图形显示和打印输出。

GrADS输出的grads格式文件可使用grads脚本进行对比出图。

#!/bin/sh 
 
level=500 
 
source ./env-grads.sh 
 
for v in {'h',} 
do 
value=$v 
file=$value'_at_'$level'.png' 
cp draw-all.gs.bak draw.gs 
sed -i "s/FILE/$file/g" draw.gs 
sed -i "s/LEVEL/$level/g" draw.gs 
sed -i "s/VALUE/$value/g" draw.gs 
grads -bpcx draw.gs 
done 

draw.gs脚本：

'reinit' 
 
infile='../post.ctl_2018110112_192' 
outdir='./out' 
outfile=outdir'/w_at_925.png' 
level=925 
value=w 
 
'open 'infile 
 
'query time' 
res=subwrd(result,3) 
 
'clear' 
'set lev 'level 
'set gxout shade2b' 
'd 'value 
'cbarn' 
'draw title 'value' at 'level'   'res 
'set grid off' 
'set grads off' 
'printim 'outfile' PNG white x800 y800'

图2 grads对比出图实例