提取区域CMIP6/isimip3b数据，并作模式不确定性区间图

利用python从isimip3b数据中提取所需的变量，并简单作图。这里以提取黄河源区2015-2100年的prsn变量（包含SSP245及SSP585）为例。

准备数据

下载数据

• 下载isimip3b数据。地址：https://data.isimip.org/search/
• 数据样式：
  
• 准备研究区域的shape文件。
  

定义文件名和路径

1. 导入需要的包

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr
import geopandas as geo
import salem

import warnings
plt.rc('font',family='Arial',size=15) 
warnings.filterwarnings("ignore")
%matplotlib inline

2. 定义文件名和文件路径

path_cmip = '/isimip3b/'
path_shape = '/outlet_WAT_wgs84/'
path_save = '/yellow/'

ssp_name = ['ssp245','ssp585']
gcm_name = ['ESM4','CM6A','ESM1','ESM2','UKESM1']
gcm_name_full = ['GFDL_ESM4','IPSL_CM6A_LR','MPI_ESM1_2_HR','MRI_ESM2_0','UKESM1_0_LL']
file_name = ['gfdl-esm4_r1i1p1f1','ipsl-cm6a-lr_r1i1p1f1','mpi-esm1-2-hr_r1i1p1f1','mri-esm2-0_r1i1p1f1','ukesm1-0-ll_r1i1p1f2']
year_name = ['2015_2020','2021_2030','2031_2040','2041_2050','2051_2060','2061_2070','2071_2080','2081_2090','2091_2100']

3. 打开示例文件，查看数据

xr.open_dataset(path_cmip+f'/{gcm_name_full[1]}/{ssp_name[0]}/global/{file_name[1]}_w5e5_{ssp_name[0]}_prsn_global_daily_{year_name[0]}.nc')

可以得到该文件的基本信息，包括经纬度和时间信息，以及最重要的变量名信息 prsn及单位 kg m-2 s-1

通过下载和查看数据可以知道，我们需要提取的prsn变量一共有2个SSP，5个不同模式（其中SSP245下GFDL_ESM4模式缺失），每个模式由9个文件拼接而成。由此可以确定思路：

1. 通过SSP创建第一层循环
2. 通过GCM创建第二层循环，并跳过SSP245下的GFDL_ESM4模式
3. 通过被分割的文件年份创建第三层循环
4. 在三层循环下，打开某个SSP，某个GCM，某个年份下的文件。利用shape文件将数据裁剪成自己需要的区域。最终需要的是流域平均的prsn，因此需要做lat和lon维的平均。最后换算单位为mm
5. 拼接数据并退出循环，保存数据。

提取目标变量数据

ishape = geo.read_file(path_shape+'yellow_river.shp')   # 读取流域shape文件

for i in range(len(ssp_name)):

    basin_data = pd.DataFrame()

    for j in range(len(gcm_name)):        
        
        # SSP245下，GFDL_ESM4 数据缺失，由此跳过
        if i == 0 and j == 0:
            
            continue

        gcm_data = np.array([])

        for k in range(len(year_name)):
            
            # 读取数据
            cmip = xr.open_dataset(path_cmip+f'/{gcm_name_full[j]}/{ssp_name[i]}/global/{file_name[j]}_w5e5_{ssp_name[i]}_prsn_global_daily_{year_name[k]}.nc')
            
            # 裁剪数据并转换单位
            cmip_b = cmip.salem.roi(shape=ishape).mean(dim=['lat','lon']).prsn.values*86400
            
            # 拼接数据
            gcm_data = np.concatenate((gcm_data, cmip_b))

        # 添加拼接好的数据
        basin_data['prsn-'+gcm_name[j]] = gcm_data
    
    # 为数据添加日期索引
    basin_data.index = pd.date_range('20150101','21001231',freq='D')

    # 保存数据
    basin_data.to_csv(path_save+f'/data/yellow_{ssp_name[i]}_prsn.csv')

代码稍加更改可以直接保存年数据，但是为了方便后续可能得研究，这里建议保留成日或者月数据，保留更多细节。

由此得到黄河源区的每日降雪量，重采样为年数据查看：

画图

划分不确定性

定义一个函数，划分不确定性，由于不同SSP包含的GCM数量不一样，这里采用*args参数来定义函数

def cal_uncertain(*args):
    
    date = pd.date_range('20150101', '21001231', freq='Y')
    data_zip = [data.values for data in args]

    p75 = []
    p25 = []
    pmean = []

    for i in range(len(date)):
        
        data_temp = [data[i] for data in data_zip]
        
        p75.append(np.percentile(data_temp, 75))
        p25.append(np.percentile(data_temp, 25))
        pmean.append(np.mean(data_temp))

    return pd.Series(p75, index=date), pd.Series(p25, index=date), pd.Series(pmean, index=date)

得到不同SSP的GCM不确定性范围

for i in range(len(ssp_name)):
    
    data = pd.read_csv(path_save+f'/data/yellow_{ssp_name[i]}_prsn.csv',index_col=0)
    
    data.index = pd.date_range('20150101', '21001231', freq='D')
    
    data = data.resample('Y').sum()
        
    locals()[ssp_name[i]+'_p75'], locals()[ssp_name[i]+'_p25'], locals()[ssp_name[i]+'_p50'] = cal_uncertain(data)

出图

1. 先简单画一张

date_y =  pd.date_range('20150101','21001231',freq='Y')
fig,ax = plt.subplots(figsize=(8,4), dpi=300)

for i in range(len(ssp_name)):

    ax.fill_between(date_y,locals()[ssp_name[i]+'_p25'],locals()[ssp_name[i]+'_p75'],)
    
    ax.plot(date_y,locals()[ssp_name[i]+'_p50'])

得到：

2. 添加亿点细节

date_y =  pd.date_range('20150101','21001231',freq='Y')
colors = ['#099dd2','#d23e09']
fig,ax = plt.subplots(figsize=(10,4), dpi=300)

for i in range(len(ssp_name)):

    ax.fill_between(date_y,
                    locals()[ssp_name[i]+'_p25'],locals()[ssp_name[i]+'_p75'],\
                    alpha=0.3,color=colors[i])
    
    ax.plot(date_y,locals()[ssp_name[i]+'_p50'],label=ssp_name[i],alpha=1,color=colors[i])

ax.set_title('Yellow')
ax.legend(loc = 'upper left',ncol=2,frameon=False,fontsize=13)
ax.set_xlabel('year')
ax.set_ylabel('Prsn (mm)')
plt.savefig(path_save+'/fig/yellow_prsn.png',format='png',bbox_inches = 'tight');  
# 这里的分号是为了出图的时候没有其他输出干扰，
# 比如<matplotlib.collections.PolyCollection at 0x7f9c594aa810>

最终成图：