timeseries_qice_gt.py 4.43 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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
    # --------------------------------------------------------------

    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 --------
    fig = plt.figure(figsize=(8, 4))
    ax = fig.add_subplot(111)

    ax.plot(time, acab, 'r-')
    ax.plot(time2, smb, 'b-')

    ax.set_ylim(0, 700)
    ax.set_xlim(1851, 2013)
    ax.xaxis.set_ticks(np.arange(1851, 2013, 26))
    ax.set_xticklabels(['1851', '1877', '1903', '1929', '1955', '1981', '2006'])
    ax.set_ylabel('SMB (kg/m2/yr)', fontsize=14)
    ax.set_xlabel('time', fontsize=14)
    # ax.text(-12,5,'T-2m, JJA',fontsize=18)
    # ax.legend([s1,s2],['>20%','>99%'],loc='lower right')

    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)

    # ------ Box plot of statistics ---------
    # If I have 3 sets of data, put them a matrix with 3 columns,
    # then boxplot(matrix) in one figure data=[a1,a2,a3]

    fig = plt.figure(figsize=(8, 4))

    # --- the whisker shows [Q1-1.5IQR, Q3+1.5IQR] without outliners
    # ax = fig.add_subplot(121)
    # ax.boxplot(acab,0,'')
    # ax.set_ylim([0,230])
    # ax.set_ylabel('acab(kg/m2/yr)',fontsize=14)
    # ax.set_xticklabels('')

    # --- the whisker shows [Q1-1.5IQR, Q3+1.5IQR] with outliners
    # --- where the Q1 is 25% percentile and Q3 is 75% percential, IQR=Q3-Q1
    # --- the outlier is data point located outside the whiskers.
    ax = fig.add_subplot(121)
    ax.boxplot(acab, sym='gd')
    ax.set_ylabel('SMB(Gt/yr)', fontsize=14)
    ax.set_ylim(-50, 230)
    ax.set_xticklabels('')
    ax = fig.add_subplot(122)
    ax.boxplot(smb, sym='gd')
    ax.set_ylabel('SMB(Gt/yr)', fontsize=14)
    ax.set_ylim(-50, 230)
    ax.set_xticklabels('')

    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)
124
    img_elem['Height'] = config['image_height']
125
126
127
    img_list.append(img_elem)

    return img_list