timeseries_qice_gt.py 5.42 KB
Newer Older
1
2
# coding=utf-8

3
4
5
6
7
8
9
10
import os
import numpy as np
import matplotlib.pyplot as plt

from netCDF4 import Dataset

from livvkit.util import elements as el

11
12
13
14
15
16
17
18
19
20
21
22
23
24
describe_ts = """
Time series of the ({model_color}) CESM surface mass balance (Gt a^-1) for the 
entire {model_start:.0f}--{model_end:.0f} simulation compared to the ({obs_color}) 
{obs_start:.0f}--{obs_end:.0f} RACMO 2.3 surface mass balance.   
"""

describe_bx = """
Box plot of the CESM surface mass balance (Gt a^-1) for the entire 
{model_start:.0f}--{model_end:.0f} simulation compared to the {obs_start:.0f}--{obs_end:.0f} 
RACMO 2.3 surface mass balance. Here, the colored lines represent the median 
values, the boxes are drawn from the first quartile to the third quartile 
(covering the interqartile range, or IQR), and the whiskers cover 1.5x the IQR.
Diamonds outside the whiskers represent suspected outliers. 
"""
25
26


27
def make_plot(config, out_path='.'):
28
    # ---------------- Data source in rhea, the two files save the same data ------------------------
29
    f_cesm_lnd_ts_qice = os.path.join(config['cesm_lnd_climos'],
30
                                      'b.e10.BG20TRCN.f09_g16.002.clm2.aavg.h0.yrly_ai.QICE.nc')
31
    f_racmo_ts_smb_aavg = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.smb.ann_tseries_aavg.nc')
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
    # --------------------------------------------------------------

    img_list = []

    # read f_cesm_lnd_ts_qice
    # time from 100-255, both time and acab have 156 data
    ncid1 = Dataset(f_cesm_lnd_ts_qice, mode='r')
    time = ncid1.variables['time'][:]
    acab = ncid1.variables['QICE'][:]  # m/yr

    # unit transfer
    acab = acab * 3.156e7  # convert from mm/s to kg/m^2/yr
    acab = acab * 1.7  # Convert to Gt/yr, for validation, compare to V13
    time = (time - time[1]) / 365 + 1851  # time from 1851-2006

    # read f_racmo_ts_smb_aavg
    ncid2 = Dataset(f_racmo_ts_smb_aavg, mode='r')
    time2 = ncid2.variables['time'][:]
    smb = ncid2.variables['smb'][:]  # mmWE = kg/m2
    smb = smb * 1.7  # unit Transfer to Gt from kg/m2

    # NOTE: add months to make the full year, divide my motnhs to get years, add to get 1958
    time2 = (time2 + 6.5) / 12 + 1957

    # maxacab = (np.max(acab))
    # minacab = (np.min(acab))
    # maxtime = (np.max(time))
    # mintime = (np.min(time))
    # maxsmb = (np.max(smb))
    # minsmb = (np.min(smb))
    # maxtime2 = (np.max(time2))
    # mintime2 = (np.min(time2))

    # print("Max QICE: {}".format(maxacab))
    # print("Min QICE: {}".format(minacab))
    # print("Max time: {}".format(maxtime))
    # print("Min time: {}".format(mintime))
    # print("Max smb Gt/yr: {}".format(maxsmb))
    # print("Min smb Gt/yr: {}".format(minsmb))
    # print("Max time2: {}".format(maxtime2))
    # print("Min time2: {}".format(mintime2))

    # ------- Time Series PLOT --------
75
    fig, ax = plt.subplots(1, 1, dpi=300, figsize=(10, 6))
76

77
78
    ax.plot(time, acab)
    ax.plot(time2, smb, linestyle='--')
79

80
    ax.set_ylim(-50, 750)
81
82
83
    ax.set_xlim(1851, 2013)
    ax.xaxis.set_ticks(np.arange(1851, 2013, 26))
    ax.set_xticklabels(['1851', '1877', '1903', '1929', '1955', '1981', '2006'])
Kennedy, Joseph H's avatar
Kennedy, Joseph H committed
84
    ax.set_ylabel('Surface mass balance (Gt a$^{-1}$)', fontsize=14)
85
86
    ax.set_xlabel('time', fontsize=14)

87
    plt.tight_layout()
88
89
90
91
92
93
94
    ts_img_path = os.path.join(out_path, 'TimeSeries_QICE_GT.png')
    plt.savefig(ts_img_path)
    plt.close()

    img_link = os.path.join(os.path.basename(out_path),
                            os.path.basename(ts_img_path))
    img_elem = el.image('TimeSeries_QICE_GT',
95
96
97
98
99
100
101
                        ' '.join(describe_ts.format(model_color='blue',
                                                    model_start=time[0],
                                                    model_end=time[-1],
                                                    obs_color='orange',
                                                    obs_start=time2[0],
                                                    obs_end=time2[-1],
                                                    ).split()),
102
103
104
105
106
                        img_link)
    if config:
        img_elem['Height'] = config['image_height']
    img_list.append(img_elem)

107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
    fig, ax = plt.subplots(1, 1, dpi=300, figsize=(8, 6))

    ln = ax.boxplot([acab, smb])

    [box.set_color('k') for box in ln['boxes']]
    [box.set_color('k') for box in ln['whiskers']]

    ln['medians'][0].set_color('C0')
    ln['fliers'][0].set_marker('d')
    ln['fliers'][0].set_markeredgecolor('C0')

    ln['medians'][1].set_color('C1')
    ln['fliers'][1].set_marker('d')
    ln['fliers'][1].set_markeredgecolor('C1')

    ax.set_xticklabels(['CESM', 'RACMO 2.3'])
Kennedy, Joseph H's avatar
Kennedy, Joseph H committed
123
    ax.set_ylabel('Surface mass balance (Gt a$^{-1}$)', fontsize=14)
124
125
    ax.set_ylim(-50, 750)
    ax.grid(False)
126

127
    plt.tight_layout()
128
129
130
131
132
133
134
    bx_img_path = os.path.join(out_path, 'BoxPlot_QICE_GT.png')
    plt.savefig(bx_img_path)
    plt.close()

    img_link = os.path.join(os.path.basename(out_path),
                            os.path.basename(bx_img_path))
    img_elem = el.image('BoxPlot_QICE_GT',
135
136
137
138
139
                        ' '.join(describe_bx.format(model_start=time[0],
                                                    model_end=time[-1],
                                                    obs_start=time2[0],
                                                    obs_end=time2[-1],
                                                    ).split()),
140
                        img_link)
141
    img_elem['Height'] = config['image_height']
142
143
144
    img_list.append(img_elem)

    return img_list