Skip to content

GitLab

Commit fa641925 authored by Kennedy, Joseph H's avatar Kennedy, Joseph H
Browse files

Update timeseries and box plots in energy extension 

Plot characteristics now match the LIVVkit 2.1 val. paper.

Fixes #26, #25
parent f7a714b5
Hide whitespace changes
Inline Side-by-side
energy/energy/timeseries_acab.py deleted 100755 → 0
View file @ f7a714b5
import os
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from livvkit.util import elements as el
describe = """TimeSeries_SMB and BoxPlot_SMB plot."""
def make_plot(config, out_path='.'):
# ---------------- Data source in rhea, the two files save the same data ------------------------
f_cesm_glc_acab_aavg = os.path.join(config['cesm_glc_climos'],
'b.e10.BG20TRCN.f09_g16.002.cism.h.acab_aavg.nc')
f_racmo_ts_smb_aavg = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.smb.ann_tseries_aavg.nc')
# --------------------------------------------------------------
img_list = []
# read f_cesm_glc_acab_aavg
# time from 100-255, both time and acab have 156 data
ncid1 = Dataset(f_cesm_glc_acab_aavg, mode='r')
time = ncid1.variables['time'][:]
acab = ncid1.variables['acab'][:] # m/yr
ncid1.set_auto_scale(True) # put the scale-factor in the time-series plot
# acab = acab #unit transfer
acab = acab * 1000 # unit transfer to kg/m2/yr
time = time + 1751 # 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/yr
ncid2.set_auto_scale(True) # put the scale-factor in the time-series plot
# NOTE: add months to make the full year, divide my motnhs to get years, add to get 1958
time2 = (time2 + 6.5) / 12 + 1957
# print(np.max(acab))
# print(np.min(acab))
# print(np.max(time))
# print(np.min(time))
# print(np.max(smb))
# print(np.min(smb))
# print(np.max(time2))
# print(np.min(time2))
# ------- 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(-10, 250)
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_SMB.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_SMB',
' '.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=(4.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(km/m2/yr)', fontsize=14)
ax.set_ylim(0, 230)
ax.set_xticklabels('')
ax = fig.add_subplot(122)
ax.boxplot(smb, sym='gd')
ax.set_ylabel('SMB(km/m2/yr)', fontsize=14)
ax.set_ylim(0, 230)
ax.set_xticklabels('')
bx_img_path = os.path.join(out_path, 'BoxPlot_SMB.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_SMB',
' '.join(describe.split()),
img_link)
img_elem['Height'] = config['image_height']
img_list.append(img_elem)
return img_list
energy/energy/timeseries_qice.py deleted 100755 → 0
View file @ f7a714b5
import os
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from livvkit.util import elements as el
describe = """TimeSeries_QICE and BoxPlot_QICE plot."""
def make_plot(config, out_path='.'):
# ---------------- Data source in rhea, the two files save the same data ------------------------
f_cesm_lnd_ts_qice = os.path.join(config['cesm_lnd_climos'],
'b.e10.BG20TRCN.f09_g16.002.clm2.aavg.h0.yrly_ai.QICE.nc')
f_racmo_ts_smb_aavg = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.smb.ann_tseries_aavg.nc')
# --------------------------------------------------------------
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
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
# 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 sum smb: {}".format(maxsmb))
# print("Min sum smb: {}".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, 1200)
ax.set_xlim(1850, 2013)
ax.xaxis.set_ticks(np.arange(1850, 2013, 26))
ax.set_xticklabels(['1850', '1876', '1902', '1928', '1954', '1980', '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.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',
' '.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(km/m2/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(km/m2/yr)', fontsize=14)
ax.set_ylim(-50, 230)
ax.set_xticklabels('')
bx_img_path = os.path.join(out_path, 'BoxPlot_QICE.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',
' '.join(describe.split()),
img_link)
img_elem['Height'] = config['image_height']
img_list.append(img_elem)
return img_list
energy/energy/timeseries_qice_gt.py
View file @ fa641925
......@@ -57,21 +57,20 @@ def make_plot(config, out_path='.'):
# print("Min time2: {}".format(mintime2))
# ------- Time Series PLOT --------
fig = plt.figure(figsize=(8, 4))
ax = fig.add_subplot(111)
fig, ax = plt.subplots(1, 1, dpi=300, figsize=(10, 6))
ax.plot(time, acab, 'r-')
ax.plot(time2, smb, 'b-')
ax.plot(time, acab)
ax.plot(time2, smb, linestyle='--')
ax.set_ylim(0, 700)
ax.set_ylim(-50, 750)
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_ylabel('Surface mass balance ($Gt$ $a^{-1}$)', 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')
plt.tight_layout()
ts_img_path = os.path.join(out_path, 'TimeSeries_QICE_GT.png')
plt.savefig(ts_img_path)
plt.close()
......@@ -85,33 +84,27 @@ def make_plot(config, out_path='.'):
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('')
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)
plt.tight_layout()
bx_img_path = os.path.join(out_path, 'BoxPlot_QICE_GT.png')
plt.savefig(bx_img_path)
plt.close()
......
energy/energy/timeseries_trefht.py
View file @ fa641925
......@@ -51,20 +51,18 @@ def make_plot(config, out_path='.'):
# print("Max t2m: {}".format(max_t2m))
# ------- Time Series PLOT --------
fig = plt.figure(figsize=(8, 4))
ax = fig.add_subplot(111)
fig, ax = plt.subplots(1, 1, dpi=300, figsize=(10, 6))
ax.plot(time, tre, 'r-')
ax.plot(time2, t2m, 'b-')
ax.plot(time, tre)
ax.plot(time2, t2m, linestyle='--')
ax.set_ylim(-22, -4)
ax.set_ylim(-24, -14)
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('T2M (C)', fontsize=14)
ax.set_ylabel('2-meter air temperature ($^{\circ}$C)', 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_T2M.png')
plt.savefig(ts_img_path)
......@@ -83,29 +81,25 @@ def make_plot(config, out_path='.'):
# 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=(4.8, 4))
# --- the whisker shows [Q1-1.5IQR, Q3+1.5IQR] without outliners
# ax = fig.add_subplot(121)
# ax.boxplot(tre,0,'')
# ax.set_ylim([0,230])
# ax.set_ylabel('tre(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(tre, sym='gd')
ax.set_ylabel('T2M(C)', fontsize=14)
ax.set_ylim(-4, -22)
ax.set_xticklabels('')
ax = fig.add_subplot(122)
ax.boxplot(t2m, sym='gd')
ax.set_ylabel('T2m(C)', fontsize=14)
ax.set_ylim(-4, -22)
ax.set_xticklabels('')
fig, ax = plt.subplots(1, 1, dpi=300, figsize=(8, 6))
ln = ax.boxplot([tre, t2m])
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('2-meter air temperature ($^{\circ}$C)', fontsize=14)
ax.set_ylim(-24, -14)
ax.grid(False)
plt.tight_layout()
bx_img_path = os.path.join(out_path, 'BoxPlot_T2M.png')
plt.savefig(bx_img_path)
plt.close()
......
energy/energy_cesm.py
View file @ fa641925
......@@ -100,8 +100,7 @@ def run(name, config):
img_list.extend(t2m_jja.make_plot(config, out_path=img_dir))
img_list.extend(ts_trefht.make_plot(config, out_path=img_dir))
img_list.extend(ts_qice_gt.make_plot(config, out_path=img_dir))
img_list.extend(ts_qice.make_plot(config, out_path=img_dir))
img_list.extend(ts_acab.make_plot(config, out_path=img_dir))
element_list.append(el.gallery('Figures', img_list))
ref_bib = utils.bib2html(config['references'])
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment