气候数据订正

免责声明

首先声明本人不是气候方向，对于气候模式完全不了解，对于订正技术完全不了解。

之所以写这篇文章纯粹是有个需求用到，在自己的初步研究下做出来了。

对于结果不做任何保证。

数据介绍

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()

经过上述步骤相信我们得到了不同变量不同模式不同情景的各类文件，只需要将文件路径传递给这个函数即可。

再次声明，本人既不了解订正过程，更没有对订正结果没有做过验证。

以上仅仅记录本人一上午的工作记录。