[METEO FRANCE DATA] add gev visualization in comparisons_visualization

experiment/meteo_france_data/stations_data/comparison_analysis.py

@@ -246,7 +246,8 @@ class ComparisonAnalysis(object):
         df_study = pd.concat([df, df_coord, df_obs], axis=1)
         return df_study
 
-    def load_main_df_merged_intersection_clean(self):
+    @cached_property
+    def df_merged_intersection_clean(self):
         df_stations = self.load_main_df_station_intersection_clean()
         df_study = self.load_main_df_study_intersection_clean()
         diff = set(df_study.columns).symmetric_difference(set(df_stations.columns))

experiment/meteo_france_data/stations_data/main_station_comparison.py

@@ -5,22 +5,30 @@ from experiment.meteo_france_data.stations_data.visualization.comparisons_visual
     ComparisonsVisualization
 
 
-def visualize_full_comparison():
+def visualize_all_stations():
     vizu = ComparisonsVisualization(altitudes=ALL_ALTITUDES_WITH_20_STATIONS_AT_LEAST, margin=150)
     vizu.visualize_maximum()
 
 
-def visualize_fast_comparison():
-    vizu = ComparisonsVisualization(altitudes=[900])
-    vizu.visualize_maximum()
+def visualize_non_nan_station():
+    vizu = ComparisonsVisualization(altitudes=ALL_ALTITUDES_WITH_20_STATIONS_AT_LEAST,
+                                    keep_only_station_without_nan_values=True,
+                                    normalize_observations=True)
+    # vizu.visualize_maximum()
+    vizu.visualize_gev()
 
 
 def example():
-    vizu = ComparisonsVisualization(altitudes=[900])
-    vizu._visualize_maximum(vizu.comparisons[0], 'Beaufortain', show=True)
+    # this is a really good example for the maxima at least
+    vizu = ComparisonsVisualization(altitudes=[900], normalize_observations=False)
+    vizu._visualize_ax_main(vizu.plot_maxima, vizu.comparisons[0], 'Beaufortain', show=True)
+
+    # vizu = ComparisonsVisualization(altitudes=[900], normalize_observations=False, keep_only_station_without_nan_values=True)
+    # vizu._visualize_ax_main(vizu.plot_gev, vizu.comparisons[0], 'Beaufortain', show=True)
 
 if __name__ == '__main__':
     # visualize_fast_comparison()
-    visualize_full_comparison()
+    # visualize_all_stations()
+    visualize_non_nan_station()
     # example()
 

experiment/meteo_france_data/stations_data/visualization/comparisons_visualization/comparisons_visualization.py

+from itertools import chain
 from typing import Dict, List
 
 import math
 import matplotlib.pyplot as plt
+import numpy as np
 import pandas as pd
 
 from experiment.meteo_france_data.scm_models_data.visualization.study_visualization.study_visualizer import \
     VisualizationParameters
 from experiment.meteo_france_data.stations_data.comparison_analysis import ComparisonAnalysis, MASSIF_COLUMN_NAME, \
     REANALYSE_STR, ALTITUDE_COLUMN_NAME
-from itertools import chain
+from extreme_estimator.extreme_models.result_from_fit import ResultFromIsmev
+from extreme_estimator.extreme_models.utils import r, safe_run_r_estimator, ro
 
 
 class ComparisonsVisualization(VisualizationParameters):
 
-    def __init__(self, altitudes=None, keep_only_station_without_nan_values=False, margin=150):
+    def __init__(self, altitudes=None, keep_only_station_without_nan_values=False, margin=150,
+                 normalize_observations=False):
         self.keep_only_station_without_nan_values = keep_only_station_without_nan_values
         if self.keep_only_station_without_nan_values:
             self.nb_columns = 5
@@ -24,7 +28,7 @@ class ComparisonsVisualization(VisualizationParameters):
         self.altitude_to_comparison = {}  # type: Dict[int, ComparisonAnalysis]
         for altitude in altitudes:
             comparison = ComparisonAnalysis(altitude=altitude,
-                                            normalize_observations=False,
+                                            normalize_observations=normalize_observations,
                                             one_station_per_massif=False,
                                             transformation_class=None,
                                             margin=margin,
@@ -43,7 +47,7 @@ class ComparisonsVisualization(VisualizationParameters):
     def massifs(self):
         return sorted(set(chain(*[c.intersection_massif_names for c in self.comparisons])))
 
-    def _visualize_main(self, visualization_ax_function, title=''):
+    def _visualize_main(self, plot_function, title=''):
         nb_rows = math.ceil(self.nb_plot / self.nb_columns)
         fig, axes = plt.subplots(nb_rows, self.nb_columns, figsize=self.figsize)
         fig.subplots_adjust(hspace=self.subplot_space, wspace=self.subplot_space)
@@ -52,39 +56,54 @@ class ComparisonsVisualization(VisualizationParameters):
         ax_idx = 0
         for massif in self.massifs:
             for c in [c for c in self.comparisons if massif in c.intersection_massif_names]:
-                visualization_ax_function(c, massif, axes[ax_idx])
+                self._visualize_ax_main(plot_function, c, massif, axes[ax_idx])
                 ax_idx += 1
         plt.suptitle(title)
         plt.show()
 
-    def visualize_maximum(self):
-        return self._visualize_main(self._visualize_maximum, 'Recent trend of Annual maxima of snowfall')
-
-    def _visualize_maximum(self, comparison: ComparisonAnalysis, massif, ax=None, show=False):
+    def _visualize_ax_main(self, plot_function, comparison: ComparisonAnalysis, massif, ax=None, show=False):
         if ax is None:
             _, ax = plt.subplots(1, 1, figsize=self.figsize)
 
-        df = comparison.load_main_df_merged_intersection_clean()
+        df = comparison.df_merged_intersection_clean.copy()
         ind = df[MASSIF_COLUMN_NAME] == massif
         df.drop([MASSIF_COLUMN_NAME], axis=1, inplace=True)
         assert sum(ind) > 0
-        df = df.loc[ind] # type: pd.DataFrame
+        df = df.loc[ind]  # type: pd.DataFrame
         colors_station = ['r', 'tab:orange', 'tab:purple', 'm', 'k']
         for color, (i, s) in zip(colors_station, df.iterrows()):
             label = i
             label += ' ({}m)'.format(s[ALTITUDE_COLUMN_NAME])
             s_values = s.iloc[3:].to_dict()
             plot_color = color if REANALYSE_STR not in label else 'g'
-            ax.plot(list(s_values.keys()), list(s_values.values()), label=label, color=plot_color)
-            ax.legend(prop={'size': 5})
+            plot_function(ax, s_values, label, plot_color)
             ax.set_title('{} at {}'.format(massif, comparison.altitude))
+            ax.legend(prop={'size': 5})
 
         if show:
-
             plt.show()
 
+    def visualize_maximum(self):
+        return self._visualize_main(self.plot_maxima, 'Recent trend of Annual maxima of snowfall')
+
+    def plot_maxima(self, ax, s_values, label, plot_color):
+        ax.plot(list(s_values.keys()), list(s_values.values()), label=label, color=plot_color)
+
     def visualize_gev(self):
-        return self._visualize_main(self._visualize_gev)
+        return self._visualize_main(self.plot_gev)
+
+    def plot_gev(self, ax, s_values, label, plot_color):
+        # todo should I normalize here ?
+        # fit gev
+        data = list(s_values.values())
+        res = safe_run_r_estimator(function=r('gev.fit'), xdat=ro.FloatVector(data),
+                                   use_start=True)
+        res = ResultFromIsmev(res, {})
+        gev_params = res.stationary_gev_params
+
+        lim = 1.5 * max(data)
+        x = np.linspace(0, lim, 1000)
+        y = gev_params.density(x)
+        # display the gev distribution that was obtained
+        ax.plot(x, y, label=label, color=plot_color)
 
-    def _visualize_gev(self):
-        pass

extreme_estimator/extreme_models/result_from_fit.py

@@ -11,6 +11,7 @@ from spatio_temporal_dataset.coordinates.abstract_coordinates import AbstractCoo
 def convertFloatVector_to_float(f):
     return np.array(f)[0]
 
+
 class ResultFromFit(object):
 
     def __init__(self, result_from_fit: robjects.ListVector) -> None:
@@ -74,6 +75,12 @@ class ResultFromIsmev(ResultFromFit):
                 i += 1
         return coef_dict
 
+    @property
+    def stationary_gev_params(self) -> GevParams:
+        params = {k.split('Coeff1')[0]: v for k, v in self.margin_coef_dict.items()
+                  if 'Coeff1' in k and 'temp' not in k}
+        return GevParams.from_dict(params)
+
     @property
     def all_parameters(self):
         return self.margin_coef_dict
@@ -91,7 +98,6 @@ class ResultFromIsmev(ResultFromFit):
         return convertFloatVector_to_float(self.name_to_value['conv']) == 0
 
 
-
 class ResultFromSpatialExtreme(ResultFromFit):
     """
     Handler from any result with the result of a fit functions from the package Spatial Extreme

extreme_estimator/margin_fits/gev/gev_params.py

@@ -21,9 +21,22 @@ class GevParams(ExtremeParams):
         else:
             return r.qgev(p, self.location, self.scale, self.shape)[0]
 
+    def density(self, x):
+        if self.has_undefined_parameters:
+            return np.nan
+        else:
+            res = r.dgev(x, self.location, self.scale, self.shape)
+            if isinstance(x, float):
+                return res[0]
+            else:
+                return np.array(res)
+
     @property
     def param_values(self):
         if self.has_undefined_parameters:
             return [np.nan for _ in range(3)]
         else:
-            return [self.location, self.scale, self.shape]
\ No newline at end of file
+            return [self.location, self.scale, self.shape]
+
+    def __str__(self):
+        return self.to_dict().__str__()
\ No newline at end of file