PDEF to nc files

[zwpzwp01]

Professor Qian,

I can read the ctl file (with PDEF of lccr) based on xgrads.

Descriptor File as follows,

DSET ^surf_HPB_m001_198101_rain_subset.dat
UNDEF -9.9999e+20
OPTIONS BIG_ENDIAN
PDEF 191 155 LCCR 35 135 97 77 30 60 135 20000 20000
XDEF 321 LINEAR 105.000000 0.2
YDEF 205 LINEAR 15.000000 0.2
ZDEF 1 LINEAR 1.000000 1
TDEF 744 LINEAR 01:00Z1JAN1981 60MN
VARS 1
rain 0 61 , 1, 0 precipitation [mm/h]
ENDVARS

One problem is 'LCCR' in this descriptor file is capital, which is different from 'lccr' used in IO.py of xgrads.

This can be solved eailsy by changing 'LCCR' to 'lccr'.

Another thing is that the 'rain 0 61' in the descriptor file gives the following error in python:

Exception: invalid storage 61, only "99" or "-1,20" are supported

Thank you for your help,

Zhao

datasets.zip

[miniufo]

Hi 小赵，你在哪里得到的数据？GrADS手册里面没看到61的描述，你看看改成99能不能打开并正确画图？我会看看把大小写关键词支持给加上。

[zwpzwp01]

钱老师，您好, 是日本的d4PDF气候模型数据，可以公开下载，但是需要注册账户

https://www.miroc-gcm.jp/d4PDF/index_en.html

我改为99后，能顺利转变为nc文件,但是转变后的数据的grid type = generic.

下面是数据的ctl文件（已经修改LCCR 为lccr, 61到99）

DSET ^surf_HPB_m001_198101_rain_subset.dat
UNDEF -9.9999e+20
OPTIONS BIG_ENDIAN
PDEF 14 11 lccr 35 135 -13 2 30 60 135 20000 20000
XDEF 321 LINEAR 105.000000 0.2
YDEF 205 LINEAR 15.000000 0.2
ZDEF 1 LINEAR 1.000000 1
TDEF 744 LINEAR 01:00Z1JAN1981 60MN
VARS 1
rain 0 99 , 1, 0 precipitation [mm/h]
ENDVARS

然后转换之后的nc数据在cdo文件中的信息为：

$ cdo griddes surf_HPB_m001_198101_rain_subset.nc

gridtype = generic
gridsize = 154
xsize = 14
ysize = 11
xname = x
yname = y
xfirst = 0
xinc = 20000
yfirst = 0
yinc = 20000

我想主要原因是在于

PDEF 191 155 Lccr 35 135 97 77 30 60 135 20000 20000
XDEF 321 LINEAR 105.000000 0.2
YDEF 205 LINEAR 15.000000 0.2
PDEF中的最后两个参数是20km（20000m），而下面的XDEF和YDEF是0.2°

所以我可能需要自定义一个lonlat的网格，但是对于61改变为99，还是不是很清楚。

如果上文有错误的话，还请老师指正！

我去联系一下这个数据的维护者，有结果的话，与您分享！

[miniufo]

GrADS在PDEF存在的时候做两个事情，一个是插值到lat/lon网格，另一个就是用插值后的数据画图。xgrads只有读取数据的功能，插值功能提供在utils.py里面。你可以试试插值以后用xarray写成nc再用cdo转。

[zwpzwp01]

钱老师，按照您说的，只需加一步转换即可，然后转为nc文件。
整体转换过来没有问题，但是，不知道为何，右下角的数据少一块，左上角的数据多一块。（请看附件）

。

[miniufo]

能看看你转换的代码吗？我用grads和你的数据画了个图如下：

理论上应该和grads一样就差不多了。但结果和你给的差很远

[zwpzwp01]

对了老师，我用的是一个很小体积的数据，进行验证的，请查看新数据（附件）
ctl文件（已经修改LCCR 为lccr, 61到99）
from xgrads import open_CtlDataset, CtlDescriptor,interp_to_latlon

ds = open_CtlDataset(r'D:\pythonwork\Japan_zhao\d4PDF_RCM\HPB\m001\surf_HPB_m001_198101_rain_subset.ctl')

obtained_ctl = CtlDescriptor(encoding="UTF-8", file='D:\pythonwork\Japan_zhao\d4PDF_RCM\HPB\m001\surf_HPB_m001_198101_rain_subset.ctl')

data = interp_to_latlon(ds, obtained_ctl)

data.to_netcdf(r'D:\pythonwork\Japan_zhao\d4PDF_RCM\HPB\m001\11.nc')

New_data.zip

[miniufo]

我用你的New_data画了第一个时刻的图，GrADS也有一个缺口，我觉得这个和插值还有画图的细节有关，应该影响不大。插值在边界上或多或少会存在不太一致的处理办法（软件内部实现），比如xarray用了线性插值，GrADS可能会用其他方法，导致哪个格点开始是缺测，有可能不完全一样。

另外python已经可以画图了，就没有必要用cdo处理再画图了，python画图可以参考这个notebook，无需插值，直接可以画在特定投影底图上。我觉得更方便。如果需要其他计算那可能还是插值一下会方便写。

[zwpzwp01]

谢谢老师的建议，我的结果与您的一模一样，就还剩参数61的问题，目前我还没有得到回复，不过根据出图效果，应该感觉不大。

[miniufo]

最近utils里面的get_coordinates_from_PDEF，interp_to_latlon，get_data_projection有比较重要的bug修复了，具体参考#37。不过你的数据不是WRF的，可能影响不大。如果没有问题的话我就关闭这个issue了。

[zwpzwp01]

谢谢老师的建议，现在没有问题了。

[ray2anya]

谢谢老师的建议，我的结果与您的一模一样，就还剩参数61的问题，目前我还没有得到回复，不过根据出图效果，应该感觉不大。

不好意思，请问一下参数61的问题有解决吗，我也是需要用d4pdf的数据，现在还是会有61的问题

[miniufo]

可以试试把61改成99？CTL手册没有详细说明61的意义，所以不知道如何处理。如果改成99也可以使用，那应该不是很重要的参数。

[ray2anya]

可以读取了，谢谢老师的回复！