In [ ]:
%pylab notebook
from __future__ import print_function
from charistools.hypsometry import Hypsometry
from charistools.timeSeries import TimeSeries
import glob
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import pandas as pd
In [ ]:
drainageIDs = ['AM_Vakhsh_at_Komsomolabad',
'IN_Hunza_at_DainyorBridge']
#'GA_Karnali_at_Benighat',
#'GA_Narayani_at_Devghat']
#'GA_SaptaKosi_at_Chatara']
drainageIDs
In [ ]:
%cd /Users/brodzik/projects/CHARIS/derived_hypsometries/MODSCAG_GF_v09_fromFile_MERRA_less_ET
%ls
In [ ]:
def plot_modeled_vs_measured(drainageID, year):
# Get the 3 melt Hyps files for this drainageID and year
meltTypes = ["exposed_glacier_ice_melt_by_elev",
"snow_on_ice_melt_by_elev",
"snow_on_land_melt_by_elev"]
for i, meltType in enumerate(meltTypes):
files = glob.glob(
"%s/%s*%4d*%s*.txt" % (
drainageID, drainageID, year, meltType))
if i == 0:
EGI_melt_file = files[0]
elif i == 1:
SOI_melt_file = files[0]
else:
SOL_melt_file = files[0]
# Read the melt hyps data
EGImelt = Hypsometry(EGI_melt_file)
SOImelt = Hypsometry(SOI_melt_file)
SOLmelt = Hypsometry(SOL_melt_file)
# Aggregate the modeled data by doy
modelHyps = Hypsometry()
for hyps in [EGImelt, SOImelt, SOLmelt]:
modelHyps.data = modelHyps.data.add(hyps.data, fill_value=0)
total_melt = modelHyps.data_by_doy().sum()
print("total melt = %.2f" % total_melt)
# Plot melt by elevation, and daily modeled melt
fig, ax = plt.subplots(3,1, figsize=(8,8))
modelHyps.imshow(ax[0], title="Modeled melt for %s (%4d)" % (drainageID, year))
model_by_doy = modelHyps.data_by_doy()
model_by_doy.plot(ax=ax[1], label="Modeled")
ax[1].set_title("Modeled melt for %s (%4d)" % (drainageID, year))
ax[1].set_ylabel("Volume ($km^3$)")
divider1 = make_axes_locatable(ax[1])
cax1 = divider1.append_axes("right", size="7%", pad=0.08)
cax1.axis('off')
# To compare to discharge by month, group the modeled melt by month
# Start a DataFrame to hold the various pieces
# Coerce the date indices to be the 15th of each month
melt_by_month = model_by_doy.groupby(
[pd.TimeGrouper('M')]).sum().to_frame(name='Modeled')
yyyymm_index = pd.to_datetime(
melt_by_month.index.map(lambda x: x.strftime('%Y-%m-15')))
df = pd.DataFrame(
data=melt_by_month.values, index=yyyymm_index, columns=['Modeled'])
df['SOLmelt'] = SOLmelt.data_by_doy().groupby([pd.TimeGrouper('M')]).sum().values
df['SOImelt'] = SOImelt.data_by_doy().groupby([pd.TimeGrouper('M')]).sum().values
df['EGImelt'] = EGImelt.data_by_doy().groupby([pd.TimeGrouper('M')]).sum().values
# Get the monthly discharge data for this basin
streamflowDir = '/Users/brodzik/projects/CHARIS/streamflow'
files = glob.glob("%s/%s.month_runoff.dat" % (streamflowDir, drainageID))
runoffFile = files[0]
print("Streamflow file: %s" % runoffFile)
runoffSeries = TimeSeries(filename=runoffFile)
this_runoff = runoffSeries.data["%4d-01-01" % year:"%4d-12-01" % year]
this_runoff[this_runoff.loc[:] < 0] = 0.
# Reset the index to 15th of each month
yyyymm_index = pd.to_datetime(
this_runoff.index.map(lambda x: x.strftime('%Y-%m-15')))
this_runoff = pd.DataFrame(
data=this_runoff.values, index=yyyymm_index, columns=['runoff'])
df["Measured"] = this_runoff["runoff"]
# I can't believe I have to reset the index like this to get the
# xaxis formatting to work properly
df.index = df.index.format(formatter=lambda x: x.strftime('%b'))
meltcolor = "blue"
runoffcolor= "gray"
df[["Measured", "Modeled"]].plot.bar(
ax=ax[2],
title="Modeled melt vs. Discharge for %s (%4d)" % (drainageID, year),
edgecolor=(0.9, 0.9, 0.9),
color=[runoffcolor, meltcolor])
ax[2].set_ylabel("Volume ($km^3$)")
divider2 = make_axes_locatable(ax[2])
cax2 = divider2.append_axes("right", size="7%", pad=0.08)
cax2.axis('off')
#ax[2].legend(loc=2)
# Save the plot
fig.tight_layout()
outfile = "%s/%s.%4d.modeled_vs_discharge.png" % (
"/Users/brodzik/ipython_notebooks/charis/",
drainageID, year)
fig.savefig(outfile)
In [ ]:
years = np.arange(14) + 2001
years
In [ ]:
for year in years:
for drainageID in drainageIDs:
print("%s: %d" % (drainageID, year))
plot_modeled_vs_measured(drainageID, year)
In [ ]:
runoffFile = "/Users/brodzik/projects/CHARIS/streamflow/IN_Hunza_at_DainyorBridge.month_runoff.dat"
runoffSeries = TimeSeries(filename=runoffFile)
this_runoff = runoffSeries.data["%4d-01-01" % year:"%4d-12-01" % year]
this_runoff
In [ ]:
def plot_mm_modeled_vs_measured(drainageID):
years = np.arange(14) + 2001
fig, ax = plt.subplots(14,1, figsize=(10,25))
for y, year in enumerate(years):
# Get the 3 melt Hyps files for this drainageID and year
meltTypes = ["exposed_glacier_ice_melt_by_elev",
"snow_on_ice_melt_by_elev",
"snow_on_land_melt_by_elev"]
for i, meltType in enumerate(meltTypes):
files = glob.glob(
"%s/%s*%4d*%s*.txt" % (
drainageID, drainageID, year, meltType))
if i == 0:
EGI_melt_file = files[0]
elif i == 1:
SOI_melt_file = files[0]
else:
SOL_melt_file = files[0]
# Read the melt hyps data
EGImelt = Hypsometry(EGI_melt_file)
SOImelt = Hypsometry(SOI_melt_file)
SOLmelt = Hypsometry(SOL_melt_file)
# Aggregate the modeled data by doy
modelHyps = Hypsometry()
for hyps in [EGImelt, SOImelt, SOLmelt]:
modelHyps.data = modelHyps.data.add(hyps.data, fill_value=0)
total_melt = modelHyps.data_by_doy().sum()
print("total melt = %.2f" % total_melt)
# To compare to discharge by month, group the modeled melt by month
# Start a DataFrame to hold the various pieces
# Coerce the date indices to be the 15th of each month
model_by_doy = modelHyps.data_by_doy()
melt_by_month = model_by_doy.groupby(
[pd.TimeGrouper('M')]).sum().to_frame(name='Modeled')
yyyymm_index = pd.to_datetime(
melt_by_month.index.map(lambda x: x.strftime('%Y-%m-15')))
df = pd.DataFrame(
data=melt_by_month.values, index=yyyymm_index, columns=['Modeled'])
# Get the monthly discharge data for this basin
streamflowDir = '/Users/brodzik/projects/CHARIS/streamflow'
files = glob.glob("%s/%s.month_runoff.dat" % (streamflowDir, drainageID))
runoffFile = files[0]
print("Streamflow file: %s" % runoffFile)
runoffSeries = TimeSeries(filename=runoffFile)
this_runoff = runoffSeries.data["%4d-01-01" % year:"%4d-12-01" % year]
this_runoff[this_runoff.loc[:] < 0] = 0.
# Reset the index to 15th of each month
yyyymm_index = pd.to_datetime(
this_runoff.index.map(lambda x: x.strftime('%Y-%m-15')))
this_runoff = pd.DataFrame(
data=this_runoff.values, index=yyyymm_index, columns=['runoff'])
df["Measured"] = this_runoff["runoff"]
# I can't believe I have to reset the index like this to get the
# xaxis formatting to work properly
df.index = df.index.format(formatter=lambda x: x.strftime('%b'))
meltcolor = "blue"
runoffcolor= "gray"
df[["Measured", "Modeled"]].plot.bar(
ax=ax[y],
title="Modeled melt vs. Discharge for %s (%4d)" % (drainageID, year),
edgecolor=(0.9, 0.9, 0.9),
color=[runoffcolor, meltcolor])
ax[y].set_ylabel("Volume ($km^3$)")
divider2 = make_axes_locatable(ax[y])
cax2 = divider2.append_axes("right", size="7%", pad=0.08)
cax2.axis('off')
# Save the plot
fig.tight_layout()
outfile = "%s/%s.mm_modeled_vs_discharge.png" % (
"/Users/brodzik/ipython_notebooks/charis/",
drainageID)
fig.savefig(outfile)
In [ ]:
for drainageID in drainageIDs:
plot_mm_modeled_vs_measured(drainageID)
In [ ]:
def plot_mm_avg_modeled_vs_measured(drainageID):
nyears = 14
#nyears = 5
years = np.arange(nyears) + 2001
# fig, ax = plt.subplots(14,1, figsize=(10,25))
for y, year in enumerate(years):
# Get the 3 melt Hyps files for this drainageID and year
meltTypes = ["exposed_glacier_ice_melt_by_elev",
"snow_on_ice_melt_by_elev",
"snow_on_land_melt_by_elev"]
for i, meltType in enumerate(meltTypes):
files = glob.glob(
"%s/%s*%4d*%s*.txt" % (
drainageID, drainageID, year, meltType))
if i == 0:
EGI_melt_file = files[0]
elif i == 1:
SOI_melt_file = files[0]
else:
SOL_melt_file = files[0]
# Read the melt hyps data
EGImelt = Hypsometry(EGI_melt_file)
SOImelt = Hypsometry(SOI_melt_file)
SOLmelt = Hypsometry(SOL_melt_file)
# Aggregate the modeled data by doy and then by month
# Should give 12 values for each of 3 files.
EGImelt_by_doy = EGImelt.data_by_doy()
EGImelt_by_month = EGImelt_by_doy.groupby(
[pd.TimeGrouper('M')]).sum().to_frame(name="EGIMelt")
SOImelt_by_doy = SOImelt.data_by_doy()
SOImelt_by_month = SOImelt_by_doy.groupby(
[pd.TimeGrouper('M')]).sum().to_frame(name="SOIMelt")
SOLmelt_by_doy = SOLmelt.data_by_doy()
SOLmelt_by_month = SOLmelt_by_doy.groupby(
[pd.TimeGrouper('M')]).sum().to_frame(name="SOLMelt")
if y > 0:
EGImelt_all = pd.concat([EGImelt_all, EGImelt_by_month])
SOImelt_all = pd.concat([SOImelt_all, SOImelt_by_month])
SOLmelt_all = pd.concat([SOLmelt_all, SOLmelt_by_month])
else:
EGImelt_all = EGImelt_by_month.copy()
SOImelt_all = SOImelt_by_month.copy()
SOLmelt_all = SOLmelt_by_month.copy()
# Figure out monthly median values
EGImelt_median = EGImelt_all.groupby(EGImelt_all.index.month).median()
SOImelt_median = SOImelt_all.groupby(SOImelt_all.index.month).median()
SOLmelt_median = SOLmelt_all.groupby(SOLmelt_all.index.month).median()
# Get the monthly discharge data for this basin
streamflowDir = '/Users/brodzik/projects/CHARIS/streamflow'
files = glob.glob("%s/%s.month_runoff.dat" % (streamflowDir, drainageID))
runoffFile = files[0]
print("runoff file=%s" % runoffFile)
runoffSeries = TimeSeries(filename=runoffFile)
runoff = runoffSeries.data["%4d-01-01" % years[0]:"%4d-12-01" % years[-1]]
runoff.replace(to_replace=[-9999.99], value=[np.nan], inplace=True)
runoff_median = runoff.groupby(runoff.index.month).median()
df = pd.concat([EGImelt_median['EGIMelt'],
SOImelt_median['SOIMelt'],
SOLmelt_median['SOLMelt'],
runoff_median['runoff']],
axis=1, keys=['EGIMelt', 'SOIMelt', 'SOLMelt', 'runoff'])
print(df)
fig, ax = plt.subplots(1,1, figsize=(10,4))
# Change the index from integers to monthnames for a nicer plot
df = df.reset_index()
df.rename(columns={'index': 'Month'}, inplace=True)
monthname = {1: 'Jan', 2: 'Feb', 3: 'Mar', 4: 'Apr', 5: 'May', 6: 'Jun',
7: 'Jul', 8: 'Aug', 9: 'Sep', 10: 'Oct', 11: 'Nov', 12: 'Dec'}
df['Month'] = df['Month'].apply(lambda x: monthname[x])
df.set_index(['Month'], inplace=True)
EGIColor = 'darkorchid'
SOIColor = 'royalblue'
SOLColor = 'seagreen'
runoffColor = (0.8, 0.8, 0.8)
edgeColor = (0.9, 0.9, 0.9)
right_bar_list = ['EGIMelt', 'SOIMelt', 'SOLMelt']
right_bar_colors = [ EGIColor, SOIColor, SOLColor ]
df[["runoff"]].plot(
ax=ax,
kind="bar",
edgecolor=edgeColor,
color=[runoffColor])
#monthly_annotation = '%s = %.2f $km^3$\nrunoff = %.2f $km^3$' % (
# right_bar_title, right_bar_sum, df["runoff"].sum())
df[right_bar_list].plot(
ax=ax,
kind="bar",
stacked=True,
position=0.,
edgecolor=edgeColor,
color=right_bar_colors)
ax.set_title("%s Median Modeled Melt vs. Discharge (%4d-%4d)" % (
drainageID, years[0], years[-1]))
ax.set_ylabel('Volume (' + r'$km^3$' + ')')
#ax[2].xaxis.set_major_formatter(plt.NullFormatter())
for container in ax.containers:
plt.setp(container, width=0.25)
#ax[2].set_ylim([0, 1.1 * ylim])
#ax[2].text(ax[2].get_xlim()[0] + (0.03 * (ax[2].get_xlim()[1] - ax[2].get_xlim()[0])),
# ax[2].get_ylim()[1] * 0.7,
# monthly_annotation,
# style='italic',
# bbox={'facecolor':'gray', 'alpha':0.1, 'pad':10})
handles, labels = ax.get_legend_handles_labels()
ax.legend(reversed(handles), reversed(labels))
# Annual labels:
annual = df.sum()
annual_annotation = 'Annual:\n%s = %.2f $km^3$\n%s = %.2f $km^3$\n%s = %.2f $km^3$\n%s = %.2f $km^3$' % (
"SOL", annual["SOLMelt"],
"SOI", annual["SOIMelt"],
"EGI", annual["EGIMelt"],
"Runoff", annual["runoff"])
ax.text(ax.get_xlim()[0] + (0.03 * (ax.get_xlim()[1] - ax.get_xlim()[0])),
ax.get_ylim()[1] * 0.55,
annual_annotation,
style='italic',
bbox={'facecolor':'gray', 'alpha':0.1, 'pad':10})
# Save the plot
fig.tight_layout()
outfile = "%s/%s.median_modeled_vs_discharge.png" % (
"/Users/brodzik/ipython_notebooks/charis/",
drainageID)
fig.savefig(outfile)
In [ ]:
for drainageID in drainageIDs:
plot_mm_avg_modeled_vs_measured(drainageID)