TEEHR Evaluation Example 3#
Hourly NWM Retrospective 3.0, CAMELS Subset (648)#
1. Get the data from S3#
For the sake of time, we prepared the individual datasets in advance and are simply copying to your 2i2c home directory. After running the cell below to copy the example_2 data.
!rm -rf ~/teehr/example-3/*
!aws s3 cp --recursive --no-sign-request s3://ciroh-rti-public-data/teehr-workshop-devcon-2024/workshop-data/example-3 ~/teehr/example-3
!tree ~/teehr/example-3/
Evaluate Model Output#
This notebook we will demonstrate how to use TEEHR to calculate metrics from a previously created joined TEEHR database containing hourly NWM3.0 Retrospective simulations and USGS observations from 1981-2022, using a range of different options for grouping and filtering. We will then create some common graphics based on the results (the same as Example 2)
In this notebook we will perform the following steps:#
- Review the contents of our joined parquet file
- Calculate metrics with different group_by options
- Calculate metrics with different filters options
- Example visualizations of TEEHR results
First setup the TEEHR class and review the contents of the joined parquet file#
from teehr.classes.duckdb_joined_parquet import DuckDBJoinedParquet
from pathlib import Path
# Define the paths to the joined parquet file and the geometry files
TEEHR_BASE = Path(Path.home(), 'teehr/example-3')
JOINED_PARQUET_FILEPATH = f"{TEEHR_BASE}/joined/configuration=nwm30_retro/variable_name=streamflow_hourly_inst/*.parquet"
GEOMETRY_FILEPATH = f"{TEEHR_BASE}/geometry/**/*.parquet"
# Initialize a teehr joined parquet class with our parquet file and geometry
joined_data = DuckDBJoinedParquet(
joined_parquet_filepath = JOINED_PARQUET_FILEPATH,
geometry_filepath = GEOMETRY_FILEPATH
)
1. Review the contents of the joined parquet files#
In practice, you may want to review the fields of data in the parquet file to plan your evaluation strategy. If the dataset is large, reading it into a dataframe may be cumbersome or even infeasible in some cases. TEEHR provides the get_joined_timeseries_schema method to quickly review the fields of the joined parquet file and the get_unique_field_values method to review the unique values contained in a specified field. The latter is particularly helpful for building dashboards for evaluation (e.g., to populate a drop down menu of possible filter or group_by values).
# Remind ourselves what fields were included
joined_data.get_joined_timeseries_schema()
# Review what configuration datasets were included
joined_data.get_unique_field_values('configuration')
# ...number of locations
len(joined_data.get_unique_field_values('primary_location_id'))
2. Calculate metrics#
%%time
gdf_all = joined_data.get_metrics(
group_by=["primary_location_id", "configuration"],
order_by=["primary_location_id", "configuration"],
include_metrics=[
'kling_gupta_efficiency_mod2',
'relative_bias',
'pearson_correlation',
'nash_sutcliffe_efficiency_normalized',
'mean_absolute_relative_error',
'primary_count'
],
include_geometry=True,
)
# view the dataframe
gdf_all
# use pandas magic to create a nice summary table of the metrics by model configuration across locations
gdf_all.groupby('configuration').describe(percentiles=[.5]).unstack(1).reset_index().rename(
columns={'level_0':'metric','level_1':'summary'}).pivot(
index=['metric','configuration'], values=0, columns='summary')
%%time
'''
Calculate metrics separately for low flows and high flows based on the
calculated field "obs_flow_category_q_mean" -> add the field to the group_by list.
'''
gdf_flowcat = joined_data.get_metrics(
group_by=["primary_location_id", "configuration", "obs_flow_category_q_mean"],
order_by=["primary_location_id", "configuration"],
include_metrics=[
'kling_gupta_efficiency_mod2',
'pearson_correlation',
'mean_absolute_relative_error',
'primary_count'
],
)
display(gdf_flowcat)
gdf_flowcat.groupby(['obs_flow_category_q_mean']).describe(percentiles=[.5]).unstack().reset_index().rename(
columns={'level_0':'metric','level_1':'summary'}).pivot(
index=['metric','obs_flow_category_q_mean'], values=0, columns='summary')
%%time
'''
Now add the location characteristics you want included in the metrics table
(for output tables and visualization)
To include location-specific attributes in the metrics table, those attributes
must be added to the group_by list. If grouping across locations (.e.g., all locations
within an RFC region), you should only add attributes that are already aggregated by that
same region (TEEHR does not check for this). An example of including location characteristic
attributes is included below.
'''
# list the attributes that are location characteristics that you want to include
# in the metrics results tables
# in the metrics results tables
include_location_characteristics = [
'aridity_none',
'runoff_ratio_none',
'baseflow_index_none',
'stream_order_none',
'q_mean_cms',
'slope_fdc_none',
'frac_urban_none',
'frac_snow_none',
'forest_frac_none',
'ecoregion_L2_none',
'river_forecast_center_none',
]
df_atts = joined_data.get_metrics(
group_by=["primary_location_id", "configuration"] + include_location_characteristics,
order_by=["primary_location_id", "configuration"],
include_metrics=[
'kling_gupta_efficiency_mod2',
'pearson_correlation',
'mean_absolute_relative_error',
'relative_bias',
'primary_count'
],
include_geometry=False,
)
# view the dataframe
display(df_atts)
# summarize just the median results across locations by attribute (river forecast center)
df_atts_summary = df_atts.groupby(['configuration','river_forecast_center_none'])\
.describe(percentiles=[.5]).unstack().unstack().reset_index()\
.rename(columns={'level_0':'metric','level_1':'summary'})
df_atts_summary[df_atts_summary['summary'].isin(['50%'])].pivot(
index=['river_forecast_center_none','configuration'],values=0, columns=['metric','summary'])
%%time
import geoviews as gv
import holoviews as hv
import colorcet as cc
hv.extension('bokeh', logo=False)
gv.extension('bokeh', logo=False)
basemap = hv.element.tiles.CartoLight()
gdf_filters = joined_data.get_metrics(
group_by=["primary_location_id", "configuration", "stream_order_none"],
order_by=["primary_location_id", "configuration"],
include_metrics=[
'kling_gupta_efficiency_mod2',
'relative_bias',
'pearson_correlation',
'nash_sutcliffe_efficiency_normalized',
'mean_absolute_relative_error',
'primary_count'
],
filters = [
{
"column": "stream_order_none",
"operator": "in",
"value": ['1','2','3','4']
#"value": ['5','6','7','8']
},
# {
# "column": "month",
# "operator": "in",
# "value": ['5','6','7','8','9']
# },
# {
# "column": "river_forecast_center_none",
# "operator": "=",
# "value": "SERFC"
# },
],
include_geometry=True,
)
#display(gdf_filters.head())
# make a quick map of locations - see how it changes as you make different filter selections
basemap * gv.Points(gdf_filters, vdims=['kling_gupta_efficiency_mod2','configuration']).select(
configuration='nwm30_retro').opts(
color='kling_gupta_efficiency_mod2',
height=400, width=600, size=7,
cmap=cc.rainbow[::-1], colorbar=True, clim=(0,1))
4. More visualizations#
# set up color and abbrevation settings to use across multiple plots
metric_abbrev=dict(
kling_gupta_efficiency_mod2 = "KGE''",
mean_absolute_relative_error = "MAE",
pearson_correlation = "Corr",
relative_bias = "Rel.Bias",
nash_sutcliffe_efficiency_normalized = "NNSE",
)
cmap_lin = cc.rainbow[::-1]
cmap_div = cc.CET_D1A[::-1]
metric_colors=dict(
kling_gupta_efficiency_mod2 = {'cmap': cmap_lin, 'clim': (0,1)},
relative_bias = {'cmap': cmap_div, 'clim': (-1,1)},
pearson_correlation = {'cmap': cmap_lin, 'clim': (0,1)},
nash_sutcliffe_efficiency_normalized = {'cmap': cmap_lin, 'clim': (0,1)},
mean_absolute_relative_error = {'cmap': cmap_lin, 'clim': (0,2)},
)
metrics = list(metric_colors.keys())
configs = ['nwm30_retro']
4a. Side by side metric maps#
First we will create side-by-side maps of the first query results above (all locations and configurations, no filters), showing metric values at each location, where dots are colored by metric value and sized by sample size. See how the comparison changes for each metric.
# map_metric = 'kling_gupta_efficiency_mod2'
# map_metric = 'pearson_correlation'
# map_metric = 'nash_sutcliffe_efficiency_normalized'
# map_metric = 'mean_absolute_relative_error'
map_metric = 'relative_bias'
# factor to size dots based on sample size
size_factor = 15/max(gdf_filters[('primary_count')])
polys = gv.Points(
gdf_all,
vdims = metrics + ['primary_location_id','configuration','primary_count'],
label = 'metric value (color), sample size (size)',
).opts(
height = 400,
width = 600,
line_color = 'gray',
colorbar = True,
size = hv.dim('primary_count') * 15/max(gdf_filters[('primary_count')]),
tools = ['hover'],
xaxis = 'bare',
yaxis = 'bare',
show_legend = True
)
maps = []
config = configs[0]
for map_metric in ['kling_gupta_efficiency_mod2','relative_bias']:
maps.append(basemap * polys.select(configuration=config).opts(
title=f"{config} | {metric_abbrev[map_metric]}",
color = map_metric,
clim = metric_colors[map_metric]['clim'],
cmap = metric_colors[map_metric]['cmap']
)
)
maps[0] + maps[1]
4b. Dataframe table and bar chart side by side#
Next we will summarize results across locations by creating a summary table with pandas (as we did above) and juxtapose it with a bar chart using holoviews and panel.
# Display dataframes and simple plots side by side using Panel
import panel as pn
gdf_summary = gdf_all.groupby('configuration').describe(percentiles=[.5]).unstack(1).reset_index().rename(
columns={'level_0':'metric','level_1':'summary'}).pivot(
index=['metric','configuration'], values=0, columns='summary')
gdf_bars = gdf_summary.drop('primary_count', axis=0)['50%'].reset_index().replace({'metric':metric_abbrev})
bars = hv.Bars(gdf_bars, kdims=['metric', 'configuration']).opts(
xrotation=90, height=400, width=300, ylabel='median',xlabel='')
pn.Row(pn.pane.DataFrame(gdf_summary, width=800), bars)
4c. Box-whisker plots of results by metric and model#
Next we’ll create box-whisker plots to see the distribution of metrics across locations for each metric and configuration.
# remove geometry so holoviews knows this is not a map.
df = gdf_all.drop('geometry', axis=1)
opts = dict(
show_legend=False,
width=100,
cmap='Set1',
xrotation=45,
labelled=[]
)
boxplots = []
for metric in metrics:
boxplots.append(
hv.BoxWhisker(df, 'configuration', metric, label=metric_abbrev[metric]).opts(
**opts,
box_fill_color=hv.dim('configuration')
)
)
hv.Layout(boxplots).cols(len(metrics))
4d. Histograms by metric and model#
Every good scientist loves a histogram. The below example creates a layout of histograms by configuration and metric, which gives us a more complete understanding of the metric distributions.
import hvplot.pandas
histograms =[]
for metric in metrics:
histograms.append(
df[df['configuration']==config].hvplot.hist(
y=metric,
ylim=(0,200),
bin_range=metric_colors[metric]['clim'],
xlabel=metric_abbrev[metric],
).opts(height = 200, width=250, title = config)
)
hv.Layout(histograms).cols(len(metrics))
4e. CDFs overlays by metric#
Every good scientist loves a CDF even more. The below example creates a layout of histograms by configuration and metric, which gives us a more complete understanding of the metric distributions. We include metrics here with (mostly) the same range (0,1) and ‘good’ value (1).
We kept this graphic for consistiency, but it is not that interesting with 1 model scenario.
import numpy as np
layout = []
for metric in [
'kling_gupta_efficiency_mod2',
'pearson_correlation',
'nash_sutcliffe_efficiency_normalized',
]:
xlim = metric_colors[metric]['clim']
xlabel = metric_abbrev[metric]
cdfs = hv.Curve([])
for config in ['nwm30_retro']:
data = df[df['configuration']==config]
data[xlabel] = np.sort(data[metric])
n = len(data[xlabel])
data['y'] = 1. * np.arange(n) / (n - 1)
cdfs = cdfs * hv.Curve(data, xlabel, 'y', label=config)
layout.append(
cdfs.opts(
width = 300,
legend_position='top_left',
xlim=xlim,
xlabel=xlabel,
title=metric_abbrev[metric],
shared_axes=False,
)
)
hv.Layout(layout).cols(5)
4f. Bar charts by attribute#
In one of the queries above, we demonstrate how to add attributes to the resulting dataframe for summary and visualization purposes. In that example we generated a summary table to RFC region. The below example uses those result to build bar charts of the median performance metric across locations within each RFC region.
df_bars = df_atts_summary.set_index('metric').drop('primary_count', axis=0).reset_index().set_index('summary').loc['50%']
df_bars = df_bars.replace({'metric': metric_abbrev}) \
.rename(columns={'river_forecast_center_none':'rfc',0:'median'}) \
.reset_index().drop('summary', axis=1)
df_bars.loc[df_bars['metric'] == 'MAE', 'median'] = 1 - df_bars.loc[df_bars['metric'] == 'MAE', 'median']
df_bars = df_bars.replace('MAE','1-MAE')
bars = hv.Bars(df_bars, kdims=['metric','configuration','rfc'], vdims=['median']).opts(
xrotation=90, height=300, width=200, ylabel='median',xlabel='')
layout = []
for rfc in df_bars['rfc'].unique():
layout.append(bars.select(rfc=rfc).opts(title=rfc))
hv.Layout(layout).cols(6)
4g Scatter plots by attribute#
Scatter plots of location metric values and location characteristics can provide insight about the relationship between the two - i.e., does model performance have a clear relationship with any of the characteristics?
# As examples, let's create scatter plots of KGE with each of the numeric attributes
import pandas as pd
import numpy as np
import geoviews as gv
import holoviews as hv
import colorcet as cc
hv.extension('bokeh', logo=False)
gv.extension('bokeh', logo=False)
basemap = hv.element.tiles.CartoLight()
location_chars = [
'aridity',
'runoff_ratio',
'baseflow_index',
'stream_order',
'q_mean_cms',
'slope_fdc',
'frac_urban',
'frac_snow',
'forest_frac'
]
df_atts.columns = df_atts.columns.str.replace('_none', '')
df_atts[location_chars] = df_atts[location_chars].apply(pd.to_numeric)
df_atts['config_num'] = np.where(df_atts['configuration']=='nwm30_retro',1,2)
metrics = [
'kling_gupta_efficiency_mod2',
'pearson_correlation',
'mean_absolute_relative_error',
'relative_bias',
]
from bokeh.models import FixedTicker
scatter_layout = []
for char in location_chars:
scatter_layout.append(
hv.Scatter(
df_atts,
kdims=[char],
vdims=['kling_gupta_efficiency_mod2', 'relative_bias', 'primary_location_id','config_num'],
label="nwm3.0"
).opts(
width = 400, height = 300,
#color = 'relative_bias',
#color = 'config_num',
#cmap = ['#377EB8', '#E41A1C'],
#colorbar = True,
clim=(0.5,2.5),
ylabel = "KGE''",
tools=['hover'],
ylim=(-1,1),
size=4,
alpha=0.8,
show_legend = True,
# colorbar_opts={
# 'ticker': FixedTicker(ticks=[1,2]),
# 'major_label_overrides': {
# 1: 'nwm30_retro',
# },
# 'major_label_text_align': 'left',
# },
))
hv.Layout(scatter_layout).opts(show_legends = True).cols(3)