diff --git a/climatedata_api/charts.py b/climatedata_api/charts.py index 90b38f5..c233900 100644 --- a/climatedata_api/charts.py +++ b/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() diff --git a/climatedata_api/download.py b/climatedata_api/download.py index a75d909..0480df7 100644 --- a/climatedata_api/download.py +++ b/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']] diff --git a/climatedata_api/utils.py b/climatedata_api/utils.py index 0953f55..7e679d0 100644 --- a/climatedata_api/utils.py +++ b/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 diff --git a/default_settings.py b/default_settings.py index b4a551d..34f9d00 100644 --- a/default_settings.py +++ b/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']