Local frequency analysis #20

sebastienlanglois · 2023-08-28T05:13:02Z

Pull Request Checklist:

This PR addresses an already opened issue (for bug fixes / features)
- This PR fixes #xyz
(If applicable) Documentation has been added / updated (for bug fixes / features).
(If applicable) Tests have been added.
HISTORY.rst has been updated (with summary of main changes).
- Link to issue (:issue:number) and pull request (:pull:number) has been added.

What kind of change does this PR introduce?

Local frequency analysis module
Notebook example

Does this PR introduce a breaking change?

No

Other information:

A few things are missing in order to complete this pull request :

local.py should be simplified and generalised. For instance, we should only assume 2 dimensions as inputs (time, id). The rest should be carried forward if desired but shouldn't be required.
Once local.py is completed, complete the notebook example
Add notebook to docs
Add tests

…ass-Data-in-frequency_analysis.local

…s-Data-in-frequency_analysis.local

github-actions · 2023-08-28T05:13:15Z

Welcome, new contributor!

It appears that this is your first Pull Request. To give credit where it's due, we ask that you add your information to the AUTHORS.rst and .zenodo.json.:

The relevant author information has been added to AUTHORS.rst and .zenodo.json.

Please make sure you've read our contributing guide. We look forward to reviewing your Pull Request shortly ✨

RondeauG · 2023-08-30T16:08:43Z

Je vois plusieurs améliorations potentielles qui pourraient être faites à local.py, par exemple en utilisant des fonctions qui existent déjà tel quel dans xclim, plutôt que de les réécrire à même la classe Data.

Pour que mes changements soient le plus explicite possible et qu'on puisse facilement les commenter, je partirais une nouvelle branche à partir de cette PR. Qu'en pensez-vous ?

sebastienlanglois · 2023-08-31T16:15:03Z

@RondeauG Oui je suis d'accord avec ton plan de match. On est bien conscient qu'il y a des optimisations à faire autant au niveau du code python qu'au niveau de la réutilisation des librairies existantes. N'hésites pas à contacter @TC-FF pour toute question sur les fonctionnalités du code.

J'ai commencé à monter un notebook (non complété) dans ce pull request pour voir plus en détails le genre d'opérations qu'on voudrait faire et éventuellement, pour ajouter à la doc

RondeauG · 2023-08-31T16:19:05Z

Alright, excellent. Je vais tenter de faire le maximums de changements d'ici la réunion de septembre.

RondeauG

Still working on my PR, but I think that I spotted an error in the code.

RondeauG · 2023-09-01T16:20:32Z

xhydro/frequency_analysis/local.py

+                .sum()
+                .where(ds.get_bool_over_tolerence(tolerence, "Volumes"), drop=True)
+            )
+            self.rm_season("Volumes")
+            # Transform tp hm³
+            # TODO add start and end and clear other attributes
+            grouped_ds = (
+                grouped_ds
+                * xr.apply_ufunc(
+                    conversion_factor_to_hm3,
+                    grouped_ds["timestep"],
+                    input_core_dims=[[]],
+                    vectorize=True,
+                )
+                * (dates[1] - dates[0])


As it is written right now, I think that you basically sum twice. The first time is at Line 508 (and would result in biased results if there are NaNs, so .mean() is more robust), and the second time is at Line 522 when you multiply the results by the number of days.

Suggested change

.sum()

.where(ds.get_bool_over_tolerence(tolerence, "Volumes"), drop=True)

)

self.rm_season("Volumes")

# Transform tp hm³

# TODO add start and end and clear other attributes

grouped_ds = (

grouped_ds

* xr.apply_ufunc(

conversion_factor_to_hm3,

grouped_ds["timestep"],

input_core_dims=[[]],

vectorize=True,

)

* (dates[1] - dates[0])

.mean()

.where(ds.get_bool_over_tolerence(tolerence, "Volumes"), drop=True)

)

self.rm_season("Volumes")

# Transform tp hm³

# TODO add start and end and clear other attributes

grouped_ds = (

grouped_ds

* xr.apply_ufunc(

conversion_factor_to_hm3,

grouped_ds["timestep"],

input_core_dims=[[]],

vectorize=True,

)

* (dates[1] - dates[0] + 1)

Let's say I have: ds2.streamflow.isel(time=slice(0,5)).values
array([[4.56, 4.25, 3.99, 3.77, 3.62]], dtype=float32)

Daily volumes: ds2.streamflow.isel(time=slice(0,5)).values * 60*60*24
array([[393984., 367200., 344736., 325728., 312768.]], dtype=float32)

Sum: ds2.streamflow.isel(time=slice(0,5)).values * 60*60*24).sum()
1744416.0 m³, or 1.744416 hm³. This should be the right answer, no?

However, ds_fa_vol = ds_fa.calculate_volume(tolerence=0.15, dates=[1, 5])
6.97766418, which is 1.744416 * (5 - 1)

@TC-FF Qu'en penses-tu ?

@RondeauG Il faut vraiment mettre en place des tests unitaires ! Mais en effet, c'est un bug. Les résultats des pointes ont étés (un peu) validés, mais pas du tout ceux des volumes.
Pour la correction, je pense que ce serait mieux de laisser sum() plutôt que de faire mean() * (dates[1] - dates[0] + 1).

Est-ce qu'il n'y a pas un danger avec la combinaison de .sum() et tolerance ? Si tu as 15% de données manquantes pour une année donnée, la somme risque d'être pas mal biaisée.

EDIT: Cela dit, facile d'ajouter les deux options.

Pour la tolérance, ça me semble ok, comme c'est par année, si on dépasse, l'année n'est pas inclue. Mais en effet si l'on est sous la tolérance, est-ce qu'on veut pondérer selon le nb de valeurs manquante ? Dans un cas de datas fixes, je dirais oui, dans un cas de dates variable, ¸a me semble moins nécessaire.

RondeauG · 2023-09-01T18:14:06Z

Je viens de terminer de regarder local.py plus en détails. Voici mon review plus "high level":

J'ai pu reproduire exactement les mêmes résultats en utilisant Data/Local et une version brouillon basées sur les fonctionnalités xclim que je compte vous proposer, ce qui est une excellente nouvelle. La seule exception est pour les volumes (voir mon commentaire plus tôt aujourd'hui).
On va pouvoir se débarrasser d'essentiellement toute la classe Data (et de l'équivalent qui se trouve dans Local), sauf peut-être pour quelques goodies. Les fonctionnalités qui s'y trouvent peuvent être remplacées par des indicateurs custom dans xclim.

# Get the maximum streamflow between DOY 1 and 120, tolerance of 15% missing values
qmax_winter = xc.core.indicator.Indicator.from_dict(
                                 data={"base": "stats",
                                            "input": {"da": "streamflow"},
                                            "parameters":
                                                            {"op": "max",
                                                             "indexer": {"doy_bounds": [1, 120]}},
                                            "missing": "pct",
                                            "missing_options": {"tolerance": 0.15}
                                            },
                                  identifier="qmax_winter",
                                  module="fa",
                              )

L'équivalent est possible aussi pour les volumes.

Pour Local, l'idée est la bonne. Tel que spécifié dans la documentation, il faut y aller un peu plus à la main si on veut autant les paramètres que les périodes de retour, et pour traiter adéquatement les données manquantes. Je crois que ça vaudrait éventuellement une PR dans xclim, mais entre temps on peut rester avec ce que vous avez. Je vais quand même probablement proposer des amélioration à la structure de la classe elle-même.

sebastienlanglois · 2023-09-06T13:44:32Z

Je viens de terminer de regarder local.py plus en détails. Voici mon review plus "high level":

J'ai pu reproduire exactement les mêmes résultats en utilisant Data/Local et une version brouillon basées sur les fonctionnalités xclim que je compte vous proposer, ce qui est une excellente nouvelle. La seule exception est pour les volumes (voir mon commentaire plus tôt aujourd'hui).

On va pouvoir se débarrasser d'essentiellement toute la classe Data (et de l'équivalent qui se trouve dans Local), sauf peut-être pour quelques goodies. Les fonctionnalités qui s'y trouvent peuvent être remplacées par des indicateurs custom dans xclim.
# Get the maximum streamflow between DOY 1 and 120, tolerance of 15% missing values
qmax_winter = xc.core.indicator.Indicator.from_dict(
                                 data={"base": "stats",
                                            "input": {"da": "streamflow"},
                                            "parameters":
                                                            {"op": "max",
                                                             "indexer": {"doy_bounds": [1, 120]}},
                                            "missing": "pct",
                                            "missing_options": {"tolerance": 0.15}
                                            },
                                  identifier="qmax_winter",
                                  module="fa",
                              )
L'équivalent est possible aussi pour les volumes.

Pour Local, l'idée est la bonne. Tel que spécifié dans la documentation, il faut y aller un peu plus à la main si on veut autant les paramètres que les périodes de retour, et pour traiter adéquatement les données manquantes. Je crois que ça vaudrait éventuellement une PR dans xclim, mais entre temps on peut rester avec ce que vous avez. Je vais quand même probablement proposer des amélioration à la structure de la classe elle-même.

Oui, ça me semble une bonne idée de ne plus avoir Data. J'avais un inconfort avec puisqu'il aurait fallu encore la redéfinir pour les analyses fréquentielles régionales et que ça pouvait devenir mélangeant.

Pour ce qui est des fonctionnalités de Local, on voulait couvrir : les périodes de retour, les critères d'ajustement (AIC, BIC, AICC, etc.), les quantiles, les méthodes d'ajustement des paramètres (max vraisemblance, L-moments, etc.), etc. @TC-FF avait aussi commencé une fonction pour visualiser le fit avec hvplot qui pourrait être ajoutée. Idéalement, on pourrait aussi ajouter une dimension aléatoire. Par exemple, on peut vouloir faire une analyse fréquentielle sur le vent, auquel cas on aurait une dimension de direction (NNE, SSE, etc.). Idéalement, dans le résultat final, on aurait nos résultats pour chaque direction dans ce cas. Aussi, il y avait l'argument calculated dont l'objectif était de passer directement les maximas à Local, car parfois ils sont précalculés.

…cy_analysis.local # Conflicts: # environment.yml # setup.py # xhydro/__init__.py

for more information, see https://pre-commit.ci

RondeauG · 2023-10-12T17:29:06Z

This PR was replaced by #27.

TC-FF and others added 24 commits May 10, 2023 16:09

replaced github username by organisation's name

f97c84a

changed layout in documentation

7a40b36

first functions

f29380f

added select_catchments function

b9dff83

Merge branch 'master' of https://github.com/hydrologie/xhydro into Cl…

7465d49

…ass-Data-in-frequency_analysis.local

custom_group_by

d09c450

custom_group_by comment in english

fd58a71

seasons related functions

46f640d

get_bool_over_tolerence function added

9abd34a

added functions to get maximum and represent them

e187afa

added imports and function to calculate volume

5057003

changed to comply with flake8 standards

c910dc2

ran black formatter

ce72f2f

Merge branch 'main' of https://github.com/hydrologie/xhydro into Clas…

c09a79a

…s-Data-in-frequency_analysis.local

added class local

62c0869

Merge branch 'main' of https://github.com/hydrologie/xhydro into Clas…

db77d2c

…s-Data-in-frequency_analysis.local

Added utils functions to go with local and data classes

60c55d1

updated init to run first Notebook

40e74ad

made changes to adrees new dataset dimensions

02b94bb

NB's first comit. Code still need upgrades for end

db600c5

change paths to general format

c123558

add xdatasets to retrieve data in example nb

14368b5

first draft notebook for local frequency analysis

a3f6ced

make frequency_analysis a module

dce5076

sebastienlanglois added the enhancement New feature or request label Aug 28, 2023

sebastienlanglois marked this pull request as draft August 28, 2023 05:16

remove first paragraph

a7bb04c

RondeauG reviewed Sep 1, 2023

View reviewed changes

RondeauG and others added 2 commits October 2, 2023 16:44

Merge remote-tracking branch 'origin/main' into Class-Data-in-frequen…

4ac2f97

…cy_analysis.local # Conflicts: # environment.yml # setup.py # xhydro/__init__.py

Merge branch 'main' into Class-Data-in-frequency_analysis.local

7b4fb8a

Zeitsperre force-pushed the Class-Data-in-frequency_analysis.local branch from 67e875e to 7b4fb8a Compare October 3, 2023 15:03

exclude jupyter notebooks

b479682

Zeitsperre force-pushed the Class-Data-in-frequency_analysis.local branch from 11ee11a to b479682 Compare October 3, 2023 15:05

[pre-commit.ci] auto fixes from pre-commit.com hooks

5b18b64

for more information, see https://pre-commit.ci

RondeauG mentioned this pull request Oct 4, 2023

major refactoring of local.py #27

Merged

5 tasks

RondeauG added the wontfix This will not be worked on label Oct 12, 2023

RondeauG closed this Oct 12, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Local frequency analysis #20

Local frequency analysis #20

sebastienlanglois commented Aug 28, 2023

github-actions bot commented Aug 28, 2023

RondeauG commented Aug 30, 2023 •

edited

sebastienlanglois commented Aug 31, 2023 •

edited

RondeauG commented Aug 31, 2023

RondeauG left a comment

RondeauG Sep 1, 2023 •

edited

RondeauG Sep 1, 2023

sebastienlanglois Sep 6, 2023

TC-FF Sep 8, 2023

RondeauG Sep 8, 2023 •

edited

TC-FF Sep 25, 2023

RondeauG commented Sep 1, 2023 •

edited

sebastienlanglois commented Sep 6, 2023

RondeauG commented Oct 12, 2023

Local frequency analysis #20

Local frequency analysis #20

Conversation

sebastienlanglois commented Aug 28, 2023

Pull Request Checklist:

What kind of change does this PR introduce?

Does this PR introduce a breaking change?

Other information:

github-actions bot commented Aug 28, 2023

RondeauG commented Aug 30, 2023 • edited

sebastienlanglois commented Aug 31, 2023 • edited

RondeauG commented Aug 31, 2023

RondeauG left a comment

Choose a reason for hiding this comment

RondeauG Sep 1, 2023 • edited

Choose a reason for hiding this comment

RondeauG Sep 1, 2023

Choose a reason for hiding this comment

sebastienlanglois Sep 6, 2023

Choose a reason for hiding this comment

TC-FF Sep 8, 2023

Choose a reason for hiding this comment

RondeauG Sep 8, 2023 • edited

Choose a reason for hiding this comment

TC-FF Sep 25, 2023

Choose a reason for hiding this comment

RondeauG commented Sep 1, 2023 • edited

sebastienlanglois commented Sep 6, 2023

RondeauG commented Oct 12, 2023

RondeauG commented Aug 30, 2023 •

edited

sebastienlanglois commented Aug 31, 2023 •

edited

RondeauG Sep 1, 2023 •

edited

RondeauG Sep 8, 2023 •

edited

RondeauG commented Sep 1, 2023 •

edited