利用grads描述文件中的pdef做插值的实用方法，转模式的同学看过来~~

东风急流

在数值模拟结束后，常会用模拟结果与实况做对比，有时候两张图比较相似，难以直观的体现出区别。这时候可以通过模拟结果与实况做差值的方法画图分析，这样可以更加具体的分析出模拟的差异（哪里大哪里小），但是fnl资料和模式输出结果在投影方式、范围上均存在一定差异，需要插值到同一分辨率网格上才能进行差值分析，现在提供一种仅仅利用grads中pdef的简单方法给大家参考下，不完善的地方还请大家帮着改进。
1、首先要理解描述文件中pdef的含义
这是我们模式结束后利用ARWpost生成的kz3.ctl中一段：
pdef  300 300 lcc  34.500  100.500  150.500  150.500  60.00000  30.00000  100.50000  18000.000  18000.000
xdef 1080 linear   56.71197   0.08108108
ydef  646 linear    7.71816   0.08108108
因为wrf中投影通常是lambert，而grads需要的是等经纬网格的数据，所以要进行投影转换，pedf告诉的原始数据的排列方式，各项含义如下
300 300 ：300×300格点
lcc:             Lambert投影
34.5，100.5 ：模式的参考点位置；
150.500 150.500 ：该点对应的网格点坐标为x方向第150.5个格点，y方向第150.5个格点；
30 60: Lambert投影中的标准纬度，
110.5：标准经度；
180000 网格距 18KM

而xdef 1080 linear 56.71197 0.08108108
ydef 646 linear 7.71816 0.08108108
则告诉我们要插值到的等经纬网格的信息，
现在比如我们要将模拟数据在10~55N，60~140E范围内与1度*1度的fnl资料做差值，在可以将上两段改为：
ydef 46 linear 10 1
xdef 81 linear 60 1

2、提取数据，将数据插值到需要的格点上
利用fwrite来提取模拟结果中的某个变量，脚本很简单，假设提取500hp高度场
'reinit'
'open E:\sk\fnl4\kz3.ctl'
'set gxout fwrite'
'set fwrite E:\sk\fnl4\wrfHGT500.dat'
'set z 11'
'set t 5'
'define a=(PH+PHB)/9.8'
'd a'
'disable fwrite'
;
这里要注意的是，我们写的数据是81*46个，非300*300个，这也是这个方法的核心，否则后面的ctl会写错！

3、写新数据的ctl文件
wrfHGT500.ctl:
dset E:\sk\fnl4\wrfHGT500.dat
title my data
undef -999000000
ydef 46 linear 10 1
xdef 81 linear 60 1
zdef 1 levels 500
tdef 1 linear 00Z08aug2015 6hr
vars 1
HGT500 1 99
endvars

4、最后是同fnl资料做差值
从fnl资料中提取出相同范围的500hp高度场方法同上，不赘述，但注意两个ctl文件要高度保持一致，否则还会出错
其描述文件为：
sk500HGT.ctl
dset e:/sk/fnl4/skHGT500.dat
undef 9.999E+20
title e:/sk/fnl4/fnl_20150809_00_00.grib2
ydef 46 linear 10 1
xdef 81 linear 60 1
tdef 1 linear 00Z08aug2015 6hr
zdef 1 levels 500
vars 1
HGT500   1 99  
ENDVARS

插值脚本主要语句为：
'open e:/sk/fnl4/sk500HGT.ctl ’
'open e:/sk/fnl4/wrf500HGT.ctl ’
'set gxout shaded’
'd HGT500.2-HGT500'
'cbarn'
贴一个我画的图，优点丑，哈哈~~

                               
登录/注册后可看大图

                               
登录/注册后可看大图

根据兰溪大神的提示，确实是我的方法太繁琐啦，抱歉给大家一些误导，简单方法总结如下：
1、修改模式结果的描述文件kz3.ctl，写成同fnl资料相同的分辨率
xdef 1080 linear 56.71197 0.08108108
ydef 646 linear 7.71816 0.08108108
现在比如我们要将模拟数据在10~55N，60~140E范围内与1度*1度的fnl资料做差值，在可以将上两段改为：
ydef 46 linear 10 1
xdef 81 linear 60 1

2、第二步是跳过提取数据步骤，直接画图，在grads中选定做差值的变量
脚本如下：
'open e:\sk\fnl4\kz3.ctl'
'open e:\sk\fnl4\fnl_20150809_00_00.ctl'
'c'
'set map 1 1 4'
'set mpdset cnworld'
'set lat 10 55'
'set lon 60 140'
'set t 5'
'set z 13'
'set xlint 10'
'set ylint 10'
'set xlopts 1 8 0.2'
'set ylopts 1 8 0.2'
'set cterp on'
'set csmooth on'
'set grads off'
'set grid off'
'set gxout shaded'
'd (PH+PHB)/9.8-HGTprs.2'
'cbarn'
'printim E:\sk\fnl4\500.png white'
;

东风急流

有同学说我这种提取wrf中数据的方法有问题，想请兰溪大神帮看看是否有问题呢？

river

就脚本上看应该是没有问题的

东风急流

river 发表于 2016-1-27 21:55
就脚本上看应该是没有问题的

恩恩，谢谢~

兰溪之水

东风急流 发表于 2016-1-27 20:39
有同学说我这种提取wrf中数据的方法有问题，想请兰溪大神帮看看是否有问题呢？

粗略看了下，应该没问题，不过并不需要这么麻烦，
只要把xdef和ydef改成与FNL的一致，然后同时打开两个文件即可做差值~

东风急流

谢谢，按你的方法已经弄好了，我把帖子修改下。

东风急流

兰溪之水 发表于 2016-1-28 08:11
粗略看了下，应该没问题，不过并不需要这么麻烦，
只要把xdef和ydef改成与FNL的一致，然后同时打开两个 ...

谢谢，按你的方法已经弄好啦，我把帖子修改一下。

香草拿铁

学习了
多谢分享

兰溪之水

东风急流 发表于 2016-1-28 10:39
谢谢，按你的方法已经弄好啦，我把帖子修改一下。

话说你的FNL数据是区域的？不是全球的？

心如止水的饭团

解决了我之前对于PDEF解释的困惑，但是我有个困惑，看到楼主的ctl,是不是也把pdef  300 300 lcc  34.500  100.500  150.500  150.500  60.00000  30.00000  100.50000  18000.000  18000.000
这行数据删除了，没有了pdef也能画图吗？这样是能读网格经纬度画没有pdef的图吗？