timeseries_qice_gt.py 4.03 KB
Newer Older
1
2
3
4
5
6
7
8
import os
import numpy as np
import matplotlib.pyplot as plt

from netCDF4 import Dataset

from livvkit.util import elements as el

9
describe = """TimeSeries_QICE_GT and BoxPlot_QICE_GT plot."""
10
11


12
def make_plot(config, out_path='.'):
13
    # ---------------- Data source in rhea, the two files save the same data ------------------------
14
    f_cesm_lnd_ts_qice = os.path.join(config['cesm_lnd_climos'],
15
                                      'b.e10.BG20TRCN.f09_g16.002.clm2.aavg.h0.yrly_ai.QICE.nc')
16
    f_racmo_ts_smb_aavg = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.smb.ann_tseries_aavg.nc')
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
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
    # --------------------------------------------------------------

    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 --------
60
    fig, ax = plt.subplots(1, 1, dpi=300, figsize=(10, 6))
61

62
63
    ax.plot(time, acab)
    ax.plot(time2, smb, linestyle='--')
64

65
    ax.set_ylim(-50, 750)
66
67
68
    ax.set_xlim(1851, 2013)
    ax.xaxis.set_ticks(np.arange(1851, 2013, 26))
    ax.set_xticklabels(['1851', '1877', '1903', '1929', '1955', '1981', '2006'])
69
    ax.set_ylabel('Surface mass balance ($Gt$ $a^{-1}$)', fontsize=14)
70
71
    ax.set_xlabel('time', fontsize=14)

72
73

    plt.tight_layout()
74
75
76
77
78
79
80
81
82
83
84
85
86
    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',
                        ' '.join(describe.split()),
                        img_link)
    if config:
        img_elem['Height'] = config['image_height']
    img_list.append(img_elem)

87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
    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'])
    ax.set_ylabel('Surface mass balance ($Gt$ $a^{-1}$)', fontsize=14)
    ax.set_ylim(-50, 750)
    ax.grid(False)
106

107
    plt.tight_layout()
108
109
110
111
112
113
114
115
116
    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',
                        ' '.join(describe.split()),
                        img_link)
117
    img_elem['Height'] = config['image_height']
118
119
120
    img_list.append(img_elem)

    return img_list