CLIM-1389: update for return periods

climatedata_api/charts.py

@@ -44,7 +44,8 @@ def _format_slices_to_highcharts_series(observations_location_slice, bccaq_locat
     chart_series['modeled_historical_range'] = convert_time_series_dataset_to_list(
         xr.merge([bccaq_location_slice_historical[f'{scenarios[0]}_{var}_p10'],
                   bccaq_location_slice_historical[f'{scenarios[0]}_{var}_p90']]), decimals)
-    # we return values in historical for all scenarios after HISTORICAL_DATE_LIMIT
+
+    # For projection data, filter to only include values after the historical period (HISTORICAL_DATE_LIMIT)
     bccaq_location_slice = bccaq_location_slice.where(
         bccaq_location_slice.time >= np.datetime64(app.config['HISTORICAL_DATE_LIMIT_AFTER'][dataset_name]), drop=True)
 
@@ -80,6 +81,55 @@ def _format_slices_to_highcharts_series(observations_location_slice, bccaq_locat
     return chart_series
 
 
+def _format_slices_to_highcharts_series_return_periods(observations_location_slice, delta_30y_slice, var, decimals, dataset_name):
+    """
+    Format observations, bccaq and delta_30y slices to dictionary of series ready for highcharts
+    """
+    scenarios = app.config['SCENARIOS'][dataset_name]
+    delta_naming = app.config['DELTA_NAMING'][dataset_name]
+    # Offset return period by 15 years (half of 30 year period),
+    # to have data points displayed at the middle of the 30y interval on the charts
+    year_offset = 15
+
+    chart_series = {
+        'observations': convert_time_series_dataset_to_list(observations_location_slice, decimals, year_offset),
+        '30y_observations': convert_time_series_dataset_to_dict(observations_location_slice, decimals, year_offset)
+    }
+
+    for scenario in scenarios:
+        # Skip the scenario if it's not available for this variable
+        if f'{scenario}_{var}_p50' not in delta_30y_slice:
+            continue
+
+        # Add chart data points
+        chart_series[scenario + '_median'] = convert_time_series_dataset_to_list(
+            delta_30y_slice[f'{scenario}_{var}_p50'], decimals, year_offset)
+        chart_series[scenario + '_range'] = convert_time_series_dataset_to_list(
+            xr.merge([delta_30y_slice[f'{scenario}_{var}_p10'],
+                      delta_30y_slice[f'{scenario}_{var}_p90']]), decimals, year_offset)
+
+        # Add 30-year changes
+        chart_series[f"delta7100_{scenario}_median"] = convert_time_series_dataset_to_dict(
+            delta_30y_slice[f'{scenario}_{var}_{delta_naming}_p50'], decimals, year_offset)
+        chart_series[f"delta7100_{scenario}_range"] = convert_time_series_dataset_to_dict(
+            xr.merge([delta_30y_slice[f'{scenario}_{var}_{delta_naming}_p10'],
+                      delta_30y_slice[f'{scenario}_{var}_{delta_naming}_p90']]), decimals, year_offset)
+
+    # Add 30-year values
+    for scenario in scenarios:
+        # Skip the scenario if it's not available for this variable
+        if f'{scenario}_{var}_p50' not in delta_30y_slice:
+            continue
+        chart_series[f"30y_{scenario}_median"] = convert_time_series_dataset_to_dict(
+            delta_30y_slice[f'{scenario}_{var}_p50'],
+            decimals, year_offset)
+        chart_series[f"30y_{scenario}_range"] = convert_time_series_dataset_to_dict(
+            xr.merge([delta_30y_slice[f'{scenario}_{var}_p10'],
+                      delta_30y_slice[f'{scenario}_{var}_p90']]), decimals, year_offset)
+
+    return chart_series
+
+
 def generate_charts(var, lat, lon, month='ann'):
     """
     Rewrite of get_values, generating a JSON ready for highcharts
@@ -127,16 +177,23 @@ def generate_charts(var, lat, lon, month='ann'):
         observations_location_slice = None
         observations_dataset = None
 
-    bccaq_dataset = open_dataset(dataset_name, 'allyears', var, msys, monthpath)
-    bccaq_location_slice = bccaq_dataset.sel(lon=loni, lat=lati, method='nearest').drop(['lat', 'lon']).dropna('time')
-
-    delta_30y_dataset = bccaq_dataset = open_dataset(dataset_name, '30ygraph', var, msys, monthpath)
+    delta_30y_dataset = open_dataset(dataset_name, '30ygraph', var, msys, monthpath)
     delta_30y_slice = delta_30y_dataset.sel(lon=loni, lat=lati, method='nearest').drop(
         [i for i in delta_30y_dataset.coords if i != 'time']).dropna('time')
 
-    chart_series = _format_slices_to_highcharts_series(observations_location_slice, bccaq_location_slice, delta_30y_slice,
-                                                       var, decimals, dataset_name)
-    bccaq_dataset.close()
+    if var.startswith('rl'):  # return period variables
+        if dataset_name != 'CMIP6':
+            return f"Bad request : return period variables only use the CMIP6 dataset, and has no {dataset_name} data available.\n", 400
+        chart_series = _format_slices_to_highcharts_series_return_periods(
+            observations_location_slice, delta_30y_slice, var, decimals, dataset_name)
+    else:
+        bccaq_dataset = open_dataset(dataset_name, 'allyears', var, msys, monthpath)
+        bccaq_location_slice = bccaq_dataset.sel(lon=loni, lat=lati, method='nearest').drop(['lat', 'lon']).dropna('time')
+
+        chart_series = _format_slices_to_highcharts_series(observations_location_slice, bccaq_location_slice, delta_30y_slice,
+                                                           var, decimals, dataset_name)
+        bccaq_dataset.close()
+
     if observations_dataset:
         observations_dataset.close()
     delta_30y_dataset.close()
@@ -284,26 +341,33 @@ def generate_regional_charts(partition, index, var, month='ann'):
         observations_dataset = None
         observations_location_slice = None
 
-    bccaq_dataset = open_dataset(dataset_name, 'allyears', var, msys, partition=partition)
-    bccaq_location_slice = bccaq_dataset.sel(geom=indexi).drop(
-        [i for i in bccaq_dataset.coords if i != 'time']).dropna('time')
-
     delta_30y_dataset = open_dataset(dataset_name, '30ygraph', var, msys, partition=partition)
     delta_30y_slice = delta_30y_dataset.sel(geom=indexi).drop(
         [i for i in delta_30y_dataset.coords if i != 'time']).dropna('time')
 
-    # we filter the appropriate month/season from the MS or QS-DEC file
-    if msys in ["MS", "QS-DEC"]:
-        bccaq_location_slice = bccaq_location_slice.sel(time=(bccaq_location_slice.time.dt.month == monthnumber))
-        if observations_location_slice:
-            observations_location_slice = observations_location_slice.sel(
-                time=(observations_location_slice.time.dt.month == monthnumber))
-        delta_30y_slice = delta_30y_slice.sel(
-            time=(delta_30y_slice.time.dt.month == monthnumber))
-
-    chart_series = _format_slices_to_highcharts_series(observations_location_slice, bccaq_location_slice, delta_30y_slice,
-                                                       var, decimals, dataset_name)
-    bccaq_dataset.close()
+    if var.startswith('rl'):  # return period variables
+        if dataset_name != 'CMIP6':
+            return f"Bad request : return period variables only use the CMIP6 dataset, and has no {dataset_name} data available.\n", 400
+        chart_series = _format_slices_to_highcharts_series_return_periods(
+            observations_location_slice, delta_30y_slice, var, decimals, dataset_name)
+    else:
+        bccaq_dataset = open_dataset(dataset_name, 'allyears', var, msys, partition=partition)
+        bccaq_location_slice = bccaq_dataset.sel(geom=indexi).drop(
+            [i for i in bccaq_dataset.coords if i != 'time']).dropna('time')
+
+        # we filter the appropriate month/season from the MS or QS-DEC file
+        if msys in ["MS", "QS-DEC"]:
+            bccaq_location_slice = bccaq_location_slice.sel(time=(bccaq_location_slice.time.dt.month == monthnumber))
+            if observations_location_slice:
+                observations_location_slice = observations_location_slice.sel(
+                    time=(observations_location_slice.time.dt.month == monthnumber))
+            delta_30y_slice = delta_30y_slice.sel(
+                time=(delta_30y_slice.time.dt.month == monthnumber))
+
+        chart_series = _format_slices_to_highcharts_series(observations_location_slice, bccaq_location_slice, delta_30y_slice,
+                                                           var, decimals, dataset_name)
+        bccaq_dataset.close()
+
     if observations_dataset:
         observations_dataset.close()
     delta_30y_dataset.close()

climatedata_api/download.py

@@ -338,6 +338,14 @@ def download():
             monthnumber = app.config['MONTH_NUMBER_LUT'][month]
             datasets[0] = datasets[0].sel(time=(datasets[0].time.dt.month == monthnumber))
         limit = app.config['SPEI_DATE_LIMIT']
+    elif var.startswith('rl'):  # return period variables
+        if dataset_name != 'CMIP6':
+            return f"Bad request : `{var}` variable only uses the CMIP6 dataset, and has no {dataset_name} data available.\n", 400
+        if dataset_type != "30ygraph":
+            return f"Bad request : `{var}` variable only has 30-year averages, use the `30ygraph` dataset_type parameter.\n", 400
+        if month != "ann":
+            return f"Bad request : `{var}` variable only supports the annual frequency, use the `ann` month parameter.\n", 400
+        datasets = [open_dataset(dataset_name, dataset_type, var, freq, monthpath)]
     else:
         if month == 'all':
             datasets = [open_dataset(dataset_name, dataset_type, var, freq, m) for m in app.config['ALLMONTHS']]

climatedata_api/utils.py

@@ -57,7 +57,7 @@ def open_dataset_by_path(path):
         raise FileNotFoundError(f"Dataset file not found: {path}")
 
 
-def convert_time_series_dataset_to_list(dataset, decimals):
+def convert_time_series_dataset_to_list(dataset, decimals, year_offset: int=0):
     """
     Converts xarray dataset to a list.
     We assume that the coordinates are timestamps, which are converted to milliseconds since 1970-01-01 (integer)
@@ -67,11 +67,11 @@ def convert_time_series_dataset_to_list(dataset, decimals):
     def _convert(float_list):
         return float_list if decimals > 0 else list(map(int, float_list))
 
-    return [[int(a[0].timestamp() * 1000)] + _convert(a[1:]) for a in
-            dataset.to_dataframe().astype('float64').round(decimals).reset_index().values.tolist()]
+    df_list = dataset.to_dataframe().astype('float64').round(decimals).reset_index().values.tolist()
+    return [[int(a[0].replace(year=a[0].year + year_offset).timestamp() * 1000)] + _convert(a[1:]) for a in df_list]
 
 
-def convert_time_series_dataset_to_dict(dataset, decimals):
+def convert_time_series_dataset_to_dict(dataset, decimals, year_offset: int=0):
     """
     Converts xarray dataset to a dict.
     We assume that the coordinates are timestamps, which are converted to milliseconds since 1970-01-01 (integer)
@@ -81,8 +81,8 @@ def convert_time_series_dataset_to_dict(dataset, decimals):
     def _convert(float_list):
         return float_list if decimals > 0 else list(map(int, float_list))
 
-    return {int(a[0].timestamp() * 1000): _convert(a[1:]) for a in
-            dataset.to_dataframe().astype('float64').round(decimals).reset_index().values.tolist()}
+    df_list = dataset.to_dataframe().astype('float64').round(decimals).reset_index().values.tolist()
+    return {int(a[0].replace(year=a[0].year + year_offset).timestamp() * 1000): _convert(a[1:]) for a in df_list}
 
 
 SAFE_CHARACTERS = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_-"  # Safe for URL

default_settings.py

@@ -125,6 +125,7 @@
              'HXmax30',
              'HXmax35',
              'HXmax40',
+             # Snow/rain fall
              'first_snowfall',
              'last_snowfall',
              'snowfall_season_length',
@@ -134,6 +135,25 @@
              'sn10mm',
              'ratot',
              'rax1day',
+             # Return periods
+             'rl10pr',
+             'rl10tasmax',
+             'rl10tasmin',
+             'rl20pr',
+             'rl20tasmax',
+             'rl20tasmin',
+             'rl25pr',
+             'rl25tasmax',
+             'rl25tasmin',
+             'rl30pr',
+             'rl30tasmax',
+             'rl30tasmin',
+             'rl50pr',
+             'rl50tasmax',
+             'rl50tasmin',
+             'rl5pr',
+             'rl5tasmax',
+             'rl5tasmin',
 ]
 
 SPEI_VARIABLES = ['spei_3m', 'spei_12m']