Skip to content
This repository has been archived by the owner on Jun 2, 2023. It is now read-only.

Edits to spatial train/val/test, additional performance metrics #211

Merged
merged 24 commits into from
Apr 26, 2023
Merged

Edits to spatial train/val/test, additional performance metrics #211

Show file tree
Hide file tree
Changes from 9 commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
de51eea
add kge for log values, and kge for upper and lower 10th of dataset
jds485 Mar 7, 2023
e30e3cc
add args to filter the data by an earliest and latest time index
jds485 Mar 8, 2023
13ee6cb
add arg to make a strict spatial partition based on the provided val …
jds485 Mar 8, 2023
e014068
removing training sites from validation data
jds485 Mar 8, 2023
f5d132f
ensuring that this spatial data filter will work if training and vali…
jds485 Mar 8, 2023
85c6917
ensure that validation sites are not in the test set
jds485 Mar 8, 2023
37fc891
add options for biweekly, monthly, and yearly timeseries summaries, a…
jds485 Mar 10, 2023
fe69f86
add metrics, edit comments
jds485 Mar 14, 2023
81dca56
add count of observations to the metrics files
jds485 Mar 14, 2023
89ec644
change groups arg to several args that describe how to group_temporal…
jds485 Mar 21, 2023
a4930cd
revert back to 10 as minimum number of observations
jds485 Mar 29, 2023
a7e35b5
change sum to mean, and change name of function arg to time_aggregati…
jds485 Mar 29, 2023
c60adc5
update comment
jds485 Mar 29, 2023
5535ba4
change variable name from sum to mean
jds485 Apr 18, 2023
4ec8543
update parameter description
jds485 Apr 18, 2023
03a8ae0
check that trn, val, and tst partitions have unique time-site ids
jds485 Apr 18, 2023
1f1d6f7
check partitions with a nan filtering
jds485 Apr 18, 2023
029a933
make the handling of site selection consistent in evaluate and prepro…
jds485 Apr 18, 2023
71d4bfd
remove sites from partitions only if sites are not provided for those…
jds485 Apr 18, 2023
4875ffc
fix line indentation
jds485 Apr 18, 2023
2b59d01
allow overriding pretraining partition check
jds485 Apr 18, 2023
7d630f8
Merge branch 'main' into jds-edits
jds485 Apr 20, 2023
6c2f18d
add small offset for log metrics
jds485 Apr 26, 2023
c2e08ad
handle case when any of trn, val, tst are not used
jds485 Apr 26, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
238 changes: 204 additions & 34 deletions river_dl/evaluate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ def rmse_logged(y_true, y_pred):

def nse_logged(y_true, y_pred):
"""
compute the rmse of the logged data
compute the nse of the logged data
:param y_true: [array-like] observed y_dataset values
:param y_pred: [array-like] predicted y_dataset values
:return: [float] the nse of the logged data
Expand All @@ -77,16 +77,32 @@ def nse_logged(y_true, y_pred):

def kge_eval(y_true, y_pred):
y_true, y_pred = filter_nan_preds(y_true, y_pred)
r, _ = pearsonr(y_pred, y_true)
mean_true = np.mean(y_true)
mean_pred = np.mean(y_pred)
std_true = np.std(y_true)
std_pred = np.std(y_pred)
r_component = np.square(r - 1)
std_component = np.square((std_pred / std_true) - 1)
bias_component = np.square((mean_pred / mean_true) - 1)
return 1 - np.sqrt(r_component + std_component + bias_component)
#Need to have > 1 observation to compute correlation.
#This could be < 2 due to percentile filtering
if len(y_true) > 1:
r, _ = pearsonr(y_pred, y_true)
mean_true = np.mean(y_true)
mean_pred = np.mean(y_pred)
std_true = np.std(y_true)
std_pred = np.std(y_pred)
r_component = np.square(r - 1)
std_component = np.square((std_pred / std_true) - 1)
bias_component = np.square((mean_pred / mean_true) - 1)
result = 1 - np.sqrt(r_component + std_component + bias_component)
else:
result = np.nan
return result

def kge_logged(y_true, y_pred):
"""
compute the kge of the logged data
:param y_true: [array-like] observed y_dataset values
:param y_pred: [array-like] predicted y_dataset values
:return: [float] the nse of the logged data
"""
y_true, y_pred = filter_nan_preds(y_true, y_pred)
y_true, y_pred = filter_negative_preds(y_true, y_pred)
return kge_eval(np.log(y_true), np.log(y_pred))

def filter_by_percentile(y_true, y_pred, percentile, less_than=True):
"""
Expand Down Expand Up @@ -136,7 +152,7 @@ def calc_metrics(df):
pred = df["pred"].values
obs, pred = filter_nan_preds(obs, pred)

if len(obs) > 10:
jds485 marked this conversation as resolved.
Show resolved Hide resolved
if len(obs) > 20:
metrics = {
"rmse": rmse_eval(obs, pred),
"nse": nse_eval(obs, pred),
Expand All @@ -162,12 +178,11 @@ def calc_metrics(df):
),
"nse_logged": nse_logged(obs, pred),
"kge": kge_eval(obs, pred),
"rmse_logged": rmse_logged(obs, pred),
"nse_top10": percentile_metric(obs, pred, nse_eval, 90, less_than=False),
"nse_bot10": percentile_metric(obs, pred, nse_eval, 10, less_than=True),
"nse_logged": nse_logged(obs, pred),
"kge_logged": kge_logged(obs, pred),
jds485 marked this conversation as resolved.
Show resolved Hide resolved
"kge_top10": percentile_metric(obs, pred, kge_eval, 90, less_than=False),
"kge_bot10": percentile_metric(obs, pred, kge_eval, 10, less_than=True),
"n_obs": len(obs)
}

else:
metrics = {
"rmse": np.nan,
Expand All @@ -182,10 +197,10 @@ def calc_metrics(df):
"nse_bot10": np.nan,
"nse_logged": np.nan,
"kge": np.nan,
"rmse_logged": np.nan,
"nse_top10": np.nan,
"nse_bot10": np.nan,
"nse_logged": np.nan,
"kge_logged": np.nan,
"kge_top10": np.nan,
"kge_bot10": np.nan,
"n_obs": len(obs)
}
return pd.Series(metrics)

Expand Down Expand Up @@ -224,27 +239,79 @@ def partition_metrics(
:param outfile: [str] file where the metrics should be written
jds485 marked this conversation as resolved.
Show resolved Hide resolved
:param val_sites: [list] sites to exclude from training and test metrics
:param test_sites: [list] sites to exclude from validation and training metrics
:param train_sites: [list] sites to exclude from test metrics
:param train_sites: [list] sites to exclude from validation and test metrics
:return: [pd dataframe] the condensed metrics
"""
var_data = fmt_preds_obs(preds, obs_file, spatial_idx_name,
time_idx_name)
var_metrics_list = []

for data_var, data in var_data.items():
#multiindex df
data_multiind = data.copy(deep=True)
jds485 marked this conversation as resolved.
Show resolved Hide resolved
data.reset_index(inplace=True)
# mask out validation and test sites from trn partition
if val_sites and partition == 'trn':
data = data[~data[spatial_idx_name].isin(val_sites)]
if test_sites and partition == 'trn':
data = data[~data[spatial_idx_name].isin(test_sites)]
# mask out test sites from val partition
if test_sites and partition=='val':
data = data[~data[spatial_idx_name].isin(test_sites)]
if train_sites and partition=='tst':
data = data[~data[spatial_idx_name].isin(train_sites)]
if val_sites and partition=='tst':
data = data[~data[spatial_idx_name].isin(val_sites)]
if train_sites and partition == 'trn':
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that there is an omission in the original code and that the following code should be added to correct that:

if train_sites and partition=='val':
    data = data[~data[spatial_idx_name].isin(train_sites)]

If I'm understanding them correctly, I think the suggested edits change the functionality of this section. As I read the original code, it appears that train_sites (as well as test_sites and val_sites) were sites that only appeared in that partition, but they weren't necessarily the only sites in that partition. In the revised code, it appears that if train_sites is specified, it will use only those sites in the training evaluations (and remove those sites from the test and validation partition).

If the intention is to change the functionality of train_sites, etc, then it's probably good to have a broader conversation. I am not using that option in my projects currently, but I don't know if others are.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As I read the original code, it appears that train_sites (as well as test_sites and val_sites) were sites that only appeared in that partition, but they weren't necessarily the only sites in that partition.

Yes, I agree. I wasn't sure if this was intentional. I could add another explicit_spatial_split parameter here to allow for the previous functionality when False and this new functionality when True. I'll hold off on that before receiving feedback from others

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds good.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh wow. Yeah. This was definitely a bug! 😮 Thanks for catching this.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@janetrbarclay: in a comment below, Jeff suggests we assume that sites within the train/val/test are the only sites in those partitions. That's also what I would expect. Do you know of anyone who is relying on the previous method wherein sites that are not within train_sites/val_sites/test_sites could be in the train/val/test partitions?

# simply use the train sites when specified.
data = data[data[spatial_idx_name].isin(train_sites)]
data_multiind = data_multiind.loc[data_multiind
.index
.get_level_values(level=spatial_idx_name)
.isin(train_sites)]
else:
#check if validation or testing sites are specified
if val_sites and partition == 'trn':
data = data[~data[spatial_idx_name].isin(val_sites)]
data_multiind = data_multiind.loc[~data_multiind
.index
.get_level_values(level=spatial_idx_name)
.isin(val_sites)]
if test_sites and partition == 'trn':
data = data[~data[spatial_idx_name].isin(test_sites)]
data_multiind = data_multiind.loc[~data_multiind
.index
.get_level_values(level=spatial_idx_name)
.isin(test_sites)]
jds485 marked this conversation as resolved.
Show resolved Hide resolved
# mask out training and test sites from val partition
if val_sites and partition == 'val':
data = data[data[spatial_idx_name].isin(val_sites)]
data_multiind = data_multiind.loc[data_multiind
.index
.get_level_values(level=spatial_idx_name)
.isin(val_sites)]
else:
if test_sites and partition=='val':
data = data[~data[spatial_idx_name].isin(test_sites)]
data_multiind = data_multiind.loc[~data_multiind
.index
.get_level_values(level=spatial_idx_name)
.isin(test_sites)]
if train_sites and partition=='val':
data = data[~data[spatial_idx_name].isin(train_sites)]
data_multiind = data_multiind.loc[~data_multiind
.index
.get_level_values(level=spatial_idx_name)
.isin(train_sites)]
# mask out training and validation sites from val partition
if test_sites and partition == 'tst':
data = data[data[spatial_idx_name].isin(test_sites)]
data_multiind = data_multiind.loc[data_multiind
.index
.get_level_values(level=spatial_idx_name)
.isin(test_sites)]
else:
if train_sites and partition=='tst':
data = data[~data[spatial_idx_name].isin(train_sites)]
data_multiind = data_multiind.loc[~data_multiind
.index
.get_level_values(level=spatial_idx_name)
.isin(train_sites)]
if val_sites and partition=='tst':
data = data[~data[spatial_idx_name].isin(val_sites)]
data_multiind = data_multiind.loc[~data_multiind
.index
.get_level_values(level=spatial_idx_name)
.isin(val_sites)]

if not group:
metrics = calc_metrics(data)
Expand All @@ -268,6 +335,99 @@ def partition_metrics(
.apply(calc_metrics)
.reset_index()
)
elif group == "year":
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is getting to be quite verbose. I suggest we split the group argument into two, maybe (group_spatially, and group_temporally). The group_spatially would be just a boolean. The group_temporally would be a str.

then the function could be something like:

if not group_spatially and not group_temporally:
    metrics = calc_metrics(data)
    # need to convert to dataframe and transpose so it looks like the
    # others
    metrics = pd.DataFrame(metrics).T
elif group_spatially and not group_temporally:
    metrics = data.groupby(spatial_idx_name).apply(calc_metrics).reset_index()
elif not group_spatially and group_temporally:
    metrics = data.groupby(pd.Grouper(index=time_idx_name, freq=group_temporally)).apply(calc_metrics).reset_index()
elif group_spatially and group_temporally:
    metrics = data.groupby([pd.Grouper(index=time_idx_name, freq=group_temporally),
                            pd.Grouper(index=spatial_idx_name)]
                            ).apply(calc_metrics).reset_index()

I think that should work. We'd just have to document how the group_temporally argument needs to work.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Definitely a bigger change, but I think it would be worth trying. It would also require propagating the change all the way up including any Snakefiles that are using this function.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would also resolve #212

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, that is much cleaner. I think one more argument would be needed for how to do temporal aggregation. I think what you've programmed would compute metrics for native timestep only (daily). I used a sum of the daily data to get biweekly, monthly, and yearly timesteps and compute metrics for those. Let me try out this edit because it will make the code much cleaner

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I addressed this change in the most recent commit. I needed 4 group args to replicate the previous functionality. I think it's more generic now to using different timesteps of aggregation, but it's more difficult to understand how to apply the 4 args to compute the desired metrics. For reference, here is the function that I used in the snakemake file to assign the 4 new args to this function to compute different metrics for different groups. See the description of the 4 args in the function docstring.

#Order in the list is: 
#group_spatially (bool), group_temporally (False or timestep to use), sum_aggregation (bool), site_based (bool)
def get_grp_arg(wildcards):
	if wildcards.metric_type == 'overall':
		return [False, False, False, False]
	elif wildcards.metric_type == 'month':
		return [False, 'M', False, False]
	elif wildcards.metric_type == 'reach':
		return [True, False, False, False]
	elif wildcards.metric_type == 'month_reach':
		return [True, 'M', False, False]
	elif wildcards.metric_type == 'monthly_site_based':
		return [False, 'M', True, True]
	elif wildcards.metric_type == 'monthly_all_sites':
		return [False, 'M', True, False]
	elif wildcards.metric_type == 'monthly_reach':
		return [True, 'M', True, False]

metrics = (
data.groupby(
data[time_idx_name].dt.year)
.apply(calc_metrics)
.reset_index()
)
elif group == ["seg_id_nat", "year"]:
metrics = (
data.groupby(
[data[time_idx_name].dt.year,
spatial_idx_name])
.apply(calc_metrics)
.reset_index()
)
elif group == "biweekly":
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the addition of biweekly and yearly options for the metrics is great.

If I'm reading this correctly, "biweekly", "monthly", and "yearly" all also use "seg_id_nat". For consistency with the other grouped metrics, it seems good to have that included in the group list. (so group = ['seg_id_nat','biweekly'])

(and as an aside, I'm noticing we should remove the hardcoded reference to seg_id_nat and replace it with spatial_idx_name. I think it's just the 3 references in this section. Would you want to fix that in this PR since you're already editing this section?)

Also, without running (which I haven't done) I'm not sure how monthly and yearly are different from ['seg_id_nat','month'] and ['seg_id_nat','year'] since they are both grouping on the same things.

Copy link
Member Author

@jds485 jds485 Mar 14, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure how monthly and yearly are different from ['seg_id_nat','month'] and ['seg_id_nat','year'] since they are both grouping on the same things

The biweekly, monthly and yearly options are resampling the daily timeseries to those time steps by taking the sum of the data within those time periods (only for the days with observations). I'm not sure that sum is the best option and am open to other suggestions.

The resulting performance metrics are computed over all reaches, not by reach as with the ['seg_id_nat','time'] options, so I can add the group option that reports these metrics by reach

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with removing the "seg_id_nat" comparison to be more generic, but that will affect snakemake workflows. For example, the workflow examples all have a function that defines group using seg_id_nat. Might be better to address this problem in a separate issue

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The resulting performance metrics are computed over all reaches, not by reach as with the ['seg_id_nat','time'] options, so I can add the group option that reports these metrics by reach

I'm not super familiar with the pandas grouper (so maybe that's the source of my confusion), but both monthly and yearly use 2 pandas groupers, 1 on time_idx_name and one on spatial_idx_name, right? So are you summing by reach and then calculating metrics across all the reaches?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So are you summing by reach and then calculating metrics across all the reaches?

yes, that's right

Copy link
Member Author

@jds485 jds485 Mar 14, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can add the group option that reports these metrics by reach

Edit: I added metrics for reach-biweekly, reach-monthly and reach-yearly timeseries.

We could also have reach-biweekly-month (summarize the biweekly timeseries by month), reach-biweekly-year, and reach-monthly-year. reach-biweekly-time would require an additional function to define a biweekly index for Python datetime objects.

#filter the data to remove nans before computing the sum
#so that the same days are being summed in the month.
data_calc = (data_multiind.dropna()
.groupby(
[pd.Grouper(level=time_idx_name, freq='2W'),
pd.Grouper(level=spatial_idx_name)])
.sum()
)
metrics = calc_metrics(data_calc)
metrics = pd.DataFrame(metrics).T
elif group == ["seg_id_nat", "biweekly"]:
#filter the data to remove nans before computing the sum
#so that the same days are being summed in the year.
data_calc = (data_multiind.dropna()
.groupby(
[pd.Grouper(level=time_idx_name, freq='2W'),
pd.Grouper(level=spatial_idx_name)])
.sum()
)
data_calc = data_calc.reset_index()
metrics = (data_calc
.groupby(spatial_idx_name)
.apply(calc_metrics)
.reset_index()
)
elif group == "monthly":
#filter the data to remove nans before computing the sum
#so that the same days are being summed in the month.
data_calc = (data_multiind.dropna()
.groupby(
[pd.Grouper(level=time_idx_name, freq='M'),
pd.Grouper(level=spatial_idx_name)])
.sum()
)
metrics = calc_metrics(data_calc)
metrics = pd.DataFrame(metrics).T
elif group == ["seg_id_nat", "monthly"]:
#filter the data to remove nans before computing the sum
#so that the same days are being summed in the year.
data_calc = (data_multiind.dropna()
.groupby(
[pd.Grouper(level=time_idx_name, freq='M'),
pd.Grouper(level=spatial_idx_name)])
.sum()
)
data_calc = data_calc.reset_index()
metrics = (data_calc
.groupby(spatial_idx_name)
.apply(calc_metrics)
.reset_index()
)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure why you would need to break it into two steps. Did you try this (below) and it didn't work?

Suggested change
elif group == ["seg_id_nat", "monthly"]:
#filter the data to remove nans before computing the sum
#so that the same days are being summed in the year.
data_calc = (data_multiind.dropna()
.groupby(
[pd.Grouper(level=time_idx_name, freq='M'),
pd.Grouper(level=spatial_idx_name)])
.sum()
)
data_calc = data_calc.reset_index()
metrics = (data_calc
.groupby(spatial_idx_name)
.apply(calc_metrics)
.reset_index()
)
elif group == ["seg_id_nat", "monthly"]:
#filter the data to remove nans before computing the sum
#so that the same days are being summed in the year.
data_calc = (data_multiind.dropna()
.groupby(
[pd.Grouper(level=time_idx_name, freq='M'),
pd.Grouper(level=spatial_idx_name)])
.apply(calc_metrics)
.reset_index()
)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No additional . functions were allowed after .sum(), so I broke it into several lines.

Also, noting that the original and suggested code compute metrics for different groups. The original code creates a monthly timeseries for each reach and computes reach-based metrics using the monthly timeseries. The suggested code uses data in the native timestep (e.g., daily) to compute metrics for each month-reach. Both of these are options in the latest commit

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess my question was, why do you need to do .sum() in the first place?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I used .sum() to convert the daily timeseries to a biweekly, monthly, and yearly timeseries for each reach. I think this is analogous to the overall metrics computed for the native timestep, which is also indexed by date-reach. Performance metrics for these timeseries indicate how well the model is able to predict at those timesteps.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay. A little slow, but I think I got it now. So it's a two-step process because first you are aggregating by the time-space dimensions (e.g., site-biweekly) and then you are aggregating by the space (e.g., site). So you end up with just one number per site, but it's looking at bigger windows of time at a time (e.g., two weeks instead of every day).

Am I getting that right?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If that's right, why did you choose sum instead of mean? It seems like mean would be better.

For example, if I do a rmse on the biweekly sums, I get

site_id
01466500      8.862246
01472104     23.094643
01472119     23.529100
014721254    20.738631
014721259    16.497109
01473499     17.143101
01473500     17.544453
01473675     17.654814
01473780     15.881052
01474500     12.225832
01480617     35.625263
01480870     36.393360
01481000     24.799420
01481500     26.861900

which is way high.

If I do it on the means:

site_id
01466500     0.633631
01472104     1.651097
01472119     1.682208
014721254    1.481826
014721259    1.179091
01473499     1.225882
01473500     1.254535
01473675     1.262302
01473780     1.136314
01474500     0.874386
01480617     2.546138
01480870     2.601188
01481000     1.772298
01481500     1.919950

A lot lower.

Is that what you'd expect?

Copy link
Member Author

@jds485 jds485 Mar 29, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So it's a two-step process because first you are aggregating by the time-space dimensions (e.g., site-biweekly) and then you are aggregating by the space (e.g., site). So you end up with just one number per site, but it's looking at bigger windows of time at a time

Yes, that's right.

why did you choose sum instead of mean?

Good point - the daily values are daily averages, not the total flow in a day. For some reason I was thinking of total flow, salinity, etc. I think the sum metrics will also be inflated compared to the mean. Here's an example for biweekly metrics for specific conductivity. Log RMSE is the same because log(sum_pred/n) - log(sum_obs/n) = log(sum_pred) - log(sum_obs)

metrics_sum
Out[126]: 
rmse               10407.172144
nse                    0.946307
rmse_top10         28507.780134
rmse_bot10          1128.871582
rmse_logged            0.162555
mean_bias          -3456.110397
mean_bias_top10   -11188.212460
mean_bias_bot10     -618.612363
nse_top10              0.793622
nse_bot10              0.803559
nse_logged             0.963944
kge                    0.899966
kge_logged             0.977389
kge_top10              0.881643
kge_bot10              0.875451
n_obs                966.000000
dtype: float64

metrics_mean
Out[127]: 
rmse                47.075655
nse                  0.649618
rmse_top10         121.434744
rmse_bot10          21.972157
rmse_logged          0.162555
mean_bias          -21.274822
mean_bias_top10    -88.179638
mean_bias_bot10    -15.051474
nse_top10           -1.496829
nse_bot10            0.258352
nse_logged           0.698174
kge                  0.756805
kge_logged           0.887974
kge_top10            0.220643
kge_bot10            0.732286
n_obs              966.000000
dtype: float64

elif group == "yearly":
#filter the data to remove nans before computing the sum
#so that the same days are being summed in the year.
data_calc = (data_multiind.dropna()
.groupby(
[pd.Grouper(level=time_idx_name, freq='Y'),
pd.Grouper(level=spatial_idx_name)])
.sum()
)
metrics = calc_metrics(data_calc)
metrics = pd.DataFrame(metrics).T
elif group == ["seg_id_nat", "yearly"]:
#filter the data to remove nans before computing the sum
#so that the same days are being summed in the year.
data_calc = (data_multiind.dropna()
.groupby(
[pd.Grouper(level=time_idx_name, freq='Y'),
pd.Grouper(level=spatial_idx_name)])
.sum()
)
data_calc = data_calc.reset_index()
metrics = (data_calc
.groupby(spatial_idx_name)
.apply(calc_metrics)
.reset_index()
)
else:
raise ValueError("group value not valid")

Expand Down Expand Up @@ -321,7 +481,17 @@ def combined_metrics(
:param group: [str or list] which group the metrics should be computed for.
Currently only supports 'seg_id_nat' (segment-wise metrics), 'month'
(month-wise metrics), ['seg_id_nat', 'month'] (metrics broken out by segment
and month), and None (everything is left together)
and month), 'year' (year-wise metrics), ['seg_id_nat', 'year']
(metrics broken out by segment and year), 'biweekly' (metrics for
biweekly timeseries, aggregated by summing data in the original timestep)
'monthly' (metrics for monthly timeseries, aggregated by summing data
in the original timestep), 'yearly' (metrics for yearly timeseries,
aggregated by summing data in the original timestep),
["seg_id_nat", "biweekly"] (metrics for biweekly timeseries broken out
by segment), ["seg_id_nat", "monthly"] (metrics for monthly timeseries broken out
by segment), ["seg_id_nat", "yearly"] (metrics for yearly timeseries
broken out by segment), and None (metrics in the original timestep computed
across all space)
:param id_dict: [dict] dictionary of id_dict where dict keys are the id
names and dict values are the id values. These are added as columns to the
metrics information
Expand Down Expand Up @@ -356,7 +526,7 @@ def combined_metrics(
group=group,
val_sites = val_sites,
test_sites = test_sites,
train_sites=train_sites)
train_sites = train_sites)
df_all.extend([metrics])

df_all = pd.concat(df_all, axis=0)
Expand Down
Loading