免责声明
首先声明本人不是气候方向,对于气候模式完全不了解,对于订正技术完全不了解。
之所以写这篇文章纯粹是有个需求用到,在自己的初步研究下做出来了。
对于结果不做任何保证。
数据介绍
54511站每天一个数据,包括多种要素,如气温、降水等
不同气候模式不同情境的模拟,时间为2015-2050
1979-2014实况数据
1979-2014气候模式模拟的历史数据
Quantile Mapping
使用该方法对气候模式预测数据进行订正,参考这篇文章,尤其要感谢作者不仅实现了该方法,还进行了开源,在github上搜索pycat即可找到,用户手册作为重要参考。
其实来说我不是很了解算法实现的细节,经过咨询相关领域专家,大体上是使用实况的概率分布去订正气候模式预报的结果。
程序编写
开源的pycat不能直接使用,其在用户手册中提供的示例以及数据均具有一定的特殊性。因为时间关系,至今我尚未搞明白在netcdf文件下写入投影信息时如何做到的。不过经过修改源码不影响对数据进行订正。
程序使用iris来读取netcdf数据,iris是英国气象局开源气象处理软件,我早期使用过,但后来发现了xarray就不再使用了...
在pycat中包含了插值过程、匹配组合历史和气候模式模拟历史数据等过程。
由于本人的数据是站点数据,因此需要先写成netcdf文件,在写入的过程中需要注意几个重要的点:
首先保证所有数据没有缺失数据,如果有请自动剔除或使用pandas的interpolate函数进行处理;
netcdf文件中coord要增加一些属性,如axis,其中在x轴要命名为大写X。
hist_nc.x.attrs.update({
'axis':'X',
'units':'degrees_east',
'standard_name':'longitude',
'long_name':'x of longitude'
})
hist_nc.y.attrs.update({
'axis':'Y',
'units':'degrees_north',
'standard_name':'latitude',
'long_name':'x of latitude'
})
hist_nc.time.attrs.update({
'axis':'T',
'standard_name':'time',
# 'units':'hours since 1970-01-01 00:00:00',
# 'calendar':'standard'
})
稍作解释,为什么对于time变量没有写入untis和calendar,因为xarray在写入时会自动匹配当前的时间数据格式。
- 由于pycat的作者每次仅对一个变量进行订正,且考虑变量的概率分布不同还可能采用不同的方法,因此在写入数据时要增加变量的各类属性。
if iname=='tt':
standard_name = 'air_temperature'
elif iname == 'mm':
standard_name = 'precipitation_amount'
else:
standard_name = 'air_temperature'</pre>
tt为温度、mm为降水,因为我这里还订正了其他变量所以其他变量的standard_name均设置为与温度一致。standard_name会影响使用的订正算法,我们来看pycat的具体代码如下:
所对应的具体方法如下:
- air_temperature
:meth:absolute_sdm
using normal distribution- precipitation_amount, surface_downwelling_shortwave_flux_in_air
:meth:relative_sdm
using gamma distribution
- 变量的单位要按照实际情况进行设置。例如温度是
celsius
,降水meters
,其他自行查阅。
修改源码
刚才讨论了,pycat示例中涉及投影情况,因此气候模拟的数据包含了3个点的数据,在进行订正前进行插值,而我这里仅仅包含一个点,且不会在netcdf中写入投影信息让iris直接能够获取,因此需要修改涉及投影这一块。首先在上一步写入netcdf中增加x和y的coord,分别设置为116.47和39.8。
修改pycat使之单独对我们这个需求应用等经纬度投影,对pycat的io模块代码进行调整如下:
这样如果在netcdf文件中的投影信息时None,则会自动设置为等经纬度投影。
修改提方法中插值部分,如果我们的数据仅仅包含1个点,则不需要插值处理。
当且仅当模式数据点不是1个点时进行插值处理。
运行订正
from pycat.io import Dataset
from pycat.esd import ScaledDistributionMapping
def esd_func(obsdir, moddir, scedir, work_dir):
#work-dir即结果输出目录,pycat默认在临时文件夹
obs = Dataset(os.path.dirname(obsdir), os.path.basename(obsdir))
mod = Dataset(os.path.dirname(moddir), os.path.basename(moddir))
sce = Dataset(os.path.dirname(scedir), os.path.basename(scedir))
qm = ScaledDistributionMapping(obs, mod, sce, work_dir=work_dir)
qm.correct()
经过上述步骤相信我们得到了不同变量
不同模式
不同情景
的各类文件,只需要将文件路径传递给这个函数即可。
再次声明,本人既不了解订正过程,更没有对订正结果没有做过验证。
以上仅仅记录本人一上午的工作记录。