diff --git a/pepsico/app_calc.py b/pepsico/app_calc.py index ff8d8323..a52b860f 100644 --- a/pepsico/app_calc.py +++ b/pepsico/app_calc.py @@ -1,4 +1,6 @@ import xarray as xr +import numpy as np +from scipy.stats import norm #This is what we should need for the app @@ -19,7 +21,9 @@ #-- then maybe some other things to beautify map tbd -def read_data(scenario, model, variable, region, unit_convert=False): +def read_data( + scenario, model, variable, region, time_res="monthly", unit_convert=False +): if region == "US-CA": xslice = slice(-154, -45) yslice = slice(60, 15) @@ -34,10 +38,14 @@ def read_data(scenario, model, variable, region, unit_convert=False): yslice = slice(28, 2) data = xr.open_zarr( f'/Data/data24/ISIMIP3b/InputData/climate/atmosphere/bias-adjusted/global' - f'/monthly/{scenario}/{model}/zarr/{variable}' + f'/{time_res}/{scenario}/{model}/zarr/{variable}' )[variable].sel(X=xslice, Y=yslice) if unit_convert : data = unit_conversion(data) + #Turns out that some models are centered on noon + if time_res == "daily" : + if data["T"][0].dt.hour == 12 : + data = data.assign_coords({"T" : data["T"] - np.timedelta64(12, "h")}) return data @@ -92,6 +100,80 @@ def seasonal_data(monthly_data, start_month, end_month, start_year=None, end_yea ) +def seasonal_wwc( + labelled_season_data, variable, frost_threshold, wet_threshold +): + # Boolean variables need the additional where to return NaNs from False/0 to Nans + # and the sum parameters for entirely NaNs seasons to remain NaNs and not turn to + # 0. This is necessary because histo and projections scenarios don't cover the + # full time period and it's not creating fake results as there are no missing + # data + if variable == "frost_days": + data_ds = ( + (labelled_season_data <= frost_threshold) + .where(~np.isnan(labelled_season_data)) + .groupby(labelled_season_data["seasons_starts"]) + .sum(skipna=True, min_count=1) + ) + wwc_units = "days" + if variable == "dry_days": + data_ds = ( + (labelled_season_data <= wet_threshold) + .where(~np.isnan(labelled_season_data)) + .groupby(labelled_season_data["seasons_starts"]) + .sum(skipna=True, min_count=1) + ) + wwc_units = "days" + if variable in ["mean_Tmax", "mean_Tmin"]: + data_ds = labelled_season_data.groupby( + labelled_season_data["seasons_starts"] + ).mean() + wwc_units = "˚C" + # It takes several 10s of minutes to get empirical quantiles maps + if variable in ["Tmax_90", "Tmin_10"]: + quantile = 0.1 if variable == "Tmin_10" else 0.9 + data_ds = labelled_season_data.groupby( + labelled_season_data["seasons_starts"] + ) + data_std = data_ds.std() + data_ds = data_ds.mean() + for var in data_ds.data_vars: + data_ds[var].data = xr.apply_ufunc( + norm.ppf, quantile, kwargs={"loc": data_ds[var], "scale": data_std[var]}, + ) + wwc_units = "˚C" + # This option is also commented out as it didn't work in its present form. Didn't + # really expect that it would though. + if variable == "frost_season_length": + data_ds = (labelled_season_data <= frost_threshold).groupby( + labelled_season_data["seasons_starts"] + ) + data_ds = ( + data_ds[-1:0].idxmax() + - data_ds.idxmax() + + np.timedelta64(1, "D") + ) + wwc_units = "days" + if variable == "wet_days": + data_ds = ( + (labelled_season_data > wet_threshold) + .where(~np.isnan(labelled_season_data)) + .groupby(labelled_season_data["seasons_starts"]) + .sum(skipna=True, min_count=1) + ) + wwc_units = "days" + # This is all a bit tedious but I didn't figure out another way to keep + # seasons_ends and renaming time dim T + # Can revisit later if this code has a future + data_ds = data_ds.rename({"seasons_starts" : "T"}) + seasons_ends = labelled_season_data["seasons_ends"].rename({"group": "T"}) + data_ds = data_ds.drop_vars(["seasons_ends", "group"]).assign_coords({"seasons_ends" : seasons_ends}) + # + for var in data_ds.data_vars: + data_ds[var].attrs["units"] = wwc_units + return data_ds + + def unit_conversion(variable): #if precipitation variable, change from kg to mm per day if variable.name == 'pr': @@ -109,6 +191,156 @@ def unit_conversion(variable): # This is in enacts calc.py + + +def daily_tobegroupedby_season( + daily_data, start_day, start_month, end_day, end_month, time_dim="T" +): + """Group daily data by season. + + Returns dataset ready to be grouped by with the daily data where all days not in season of interest are dropped. + + season_starts: + An array where the non-dropped days are indexed by the first day of their season. + -- to use to groupby + seasons_ends: + An array with the dates of the end of the seasons. + Can then apply groupby on daily_data against seasons_starts, and preserving seasons_ends for the record. + + If starting day-month is 29-Feb, uses 1-Mar. + + If ending day-month is 29-Feb, uses 1-Mar and uses < rather than <= + That means that the last day included in the season will be 29-Feb in leap years and 28-Feb otherwise. + + Parameters + ----------- + daily_data : DataArray + Daily data to be grouped. + start_day : int + Day of the start date of the season. + start_month : int + Month of the start date of the season. + end_day : int + Day of the end date of the season. + end_month : int + Day of the end date of the season. + time_dim : str, optional + Time coordinate in `daily_data` (default `time_dim`="T"). + Returns + ------- + daily_tobegroupedby_season : Dataset + Daily data grouped by season using season start date. Dataset includes grouped data + and `season_starts`, `season_ends` as output variables. + See Also + -------- + Notes + ----- + Examples + -------- + """ + # Deal with leap year cases + if start_day == 29 and start_month == 2: + start_day = 1 + start_month = 3 + # Find seasons edges + start_edges = sel_day_and_month(daily_data[time_dim], start_day, start_month) + if end_day == 29 and end_month == 2: + end_edges = sel_day_and_month(daily_data[time_dim], 1 , 3, offset=-1) + else: + end_edges = sel_day_and_month(daily_data[time_dim], end_day, end_month) + # Drop dates outside very first and very last edges + # -- this ensures we get complete seasons with regards to edges, later on + daily_data = daily_data.sel(**{time_dim: slice(start_edges[0], end_edges[-1])}) + start_edges = start_edges.sel(**{time_dim: slice(start_edges[0], end_edges[-1])}) + end_edges = end_edges.sel( + **{time_dim: slice(start_edges[0], end_edges[-1])} + ).assign_coords(**{time_dim: start_edges[time_dim]}) + # Drops daily data not in seasons of interest + days_in_season = ( + daily_data[time_dim] >= start_edges.rename({time_dim: "group"}) + ) & (daily_data[time_dim] <= end_edges.rename({time_dim: "group"})) + days_in_season = days_in_season.sum(dim="group") + daily_data = daily_data.where(days_in_season == 1, drop=True) + # Creates seasons_starts that will be used for grouping + # and seasons_ends that is one of the outputs + seasons_groups = (daily_data[time_dim].dt.day == start_day) & ( + daily_data[time_dim].dt.month == start_month + ) + seasons_groups = seasons_groups.cumsum() - 1 + seasons_starts = ( + start_edges.rename({time_dim: "toto"})[seasons_groups] + .drop_vars("toto") + .rename("seasons_starts") + ) + seasons_ends = end_edges.rename({time_dim: "group"}).rename("seasons_ends") + # I actually changed this from enacts by making both seasons_starts/ends cords + # rather than vars. This seems better following a comment Aaron made on xcdat and + # because we local data case provides an xr.Dataset of vars to groupby and we + # don't want seasons_starts/ends to be groupedby + daily_tobegroupedby_season = daily_data.assign_coords( + {"seasons_starts" : seasons_starts, "seasons_ends" : seasons_ends} + ) + return daily_tobegroupedby_season + + +def sel_day_and_month(daily_dim, day, month, offset=0): + """Return a subset of `daily_dim` daily time dimension of corresponding + `day`/`month` + `offset` day(s) for all years. + + The returned time dimension can then be used to select daily DataArrays. + Offset is convenient to get days prior to a 1st of March, + that are not identifiable by a common `day` (28, 29). + + Parameters + ---------- + daily_dim : DataArray[datetime64[ns]] + A daily time dimension. + day : int + day of the `month`. + month : int + month of the year. + offset : int, optional + number of days to add to `day`/`month` to offset the selection + (the default is 0, which implies no offset). + + Returns + ------- + DataArray[datetime64[ns]] + a subset of `daily_dim` with all and only `day`-`month` points, offset by `offset` days. + + See Also + -------- + + Examples + -------- + >>> t = pd.date_range(start="2000-05-01", end="20002-04-30", freq="1D") + >>> values = numpy.arrange(t.size) + >>> toto = xarray.DataArray(numpy.arrange(61), dims=["T"], coords={"T": t}) + >>> sel_day_and_month(toto["T"], 6, 5) + + array(['2000-05-06T00:00:00.000000000', '2001-05-06T00:00:00.000000000',] + dtype='datetime64[ns]') + Coordinates: + * T (T) datetime64[ns] 2000-05-06 2001-05-06 + + With an offset of -1 day + + >>> t = pd.date_range(start="2000-01-01", end="20002-01-30", freq="1D") + >>> toto = xarray.DataArray(numpy.arrange(t.size), dims=["T"], coords={"T": t}) + >>> sel_day_and_month(toto["T"], 1, 3, -1) + + array(['2000-02-29T00:00:00.000000000', '2001-02-28T00:00:00.000000000',] + dtype='datetime64[ns]') + Coordinates: + * T (T) datetime64[ns] 2000-02-29 2001-02-28 + """ + return daily_dim.where( + lambda x: ((x - np.timedelta64(offset, "D")).dt.day == day) + & ((x - np.timedelta64(offset, "D")).dt.month == month), + drop=True + ) + + def strftimeb2int(strftimeb): """Convert month values to integers (1-12) from strings. diff --git a/pepsico/config-defaults.yaml b/pepsico/config-defaults.yaml index 4df5fed6..6baa20da 100644 --- a/pepsico/config-defaults.yaml +++ b/pepsico/config-defaults.yaml @@ -92,3 +92,7 @@ maprooms: # App set up - title: Projections core_path: projections + + proj_wwc: + - title: WWC Projections + core_path: proj_wwc diff --git a/pepsico/config-dev-sample.yaml b/pepsico/config-dev-sample.yaml index f71905db..c361213a 100644 --- a/pepsico/config-dev-sample.yaml +++ b/pepsico/config-dev-sample.yaml @@ -15,4 +15,8 @@ maprooms: # App set up - title: Projections core_path: projections + + proj_wwc: + - title: WWC Projections + core_path: proj_wwc \ No newline at end of file diff --git a/pepsico/proj_wwc/__init__.py b/pepsico/proj_wwc/__init__.py new file mode 100644 index 00000000..0379666a --- /dev/null +++ b/pepsico/proj_wwc/__init__.py @@ -0,0 +1 @@ +from .maproom import register diff --git a/pepsico/proj_wwc/layout.py b/pepsico/proj_wwc/layout.py new file mode 100644 index 00000000..c069bcaf --- /dev/null +++ b/pepsico/proj_wwc/layout.py @@ -0,0 +1,216 @@ +from dash import html +from fieldsets import Block, Select, PickPoint, Month, Number +import layout_utilities as lou + + +from globals_ import GLOBAL_CONFIG + +IRI_BLUE = "rgb(25,57,138)" +IRI_GRAY = "rgb(113,112,116)" +LIGHT_GRAY = "#eeeeee" + + +def app_layout(): + + return lou.app_1( + + lou.navbar("CCA - WWC", + Block("Region", + Select( + id="region", + options=["Thailand", "US-CA"], + labels=["Thailand", "United States and Canada"], + init=1, + ), + ), + PickPoint(width="8em"), + Block("Submit", + Block("Scenario", Select( + id="scenario", + options=["picontrol", "ssp126", "ssp370", "ssp585"], + init=1, + )), + Block("Model", Select(id="model", options=[ + "Multi-Model-Average", "GFDL-ESM4", "IPSL-CM6A-LR", + "MPI-ESM1-2-HR", "MRI-ESM2-0", "UKESM1-0-LL", + ])), + Block("WWC statistics", Select( + id="variable", + options=[ + # "warm_nights", + # "rain_events", + "mean_Tmax", + "mean_Tmin", + "dry_days", + "wet_days", + # "heatwave_duration", + #"frost_season_length", + "frost_days", + "Tmax_90", + "Tmin_10", + # "longest_dry_spell", + # "longest_wet_spell", + # "wet_day_persistence", + # "dry_day_persistence", + # "dry_spells_mean_length", + # "dry_spells_median_length", + ], + labels=[ + # "Warm Nights", + # "Count of Rain Events", + "Mean Max Temperature", + "Mean Min Temperature", + "Count of Dry Days", + "Count of Wet Days", + # "Mean Heatwaves Duration", + #"Frost Season Length", + "Count of Frost Days", + "Max Temperature 90th %-ile", + "Min Temperature 10th %-ile", + # "Longest Dry Spell", + # "Longest Wet Spell", + # "Mean Wet Day Persistence", + # "Mean Dry Day Persistence", + # "Mean Dry Spells Length", + # "Median Dry Spells Length", + ], + init=0, + )), + Block("Definitions", + "Frost <=", + Number( + id="frost", + default=-2, + min=-99, + max=0, + width="5em", + debounce=False, + ), + "˚C; ", + "Dry/Wet day <= / >", + Number( + id="wet", + default=0, + min=0, + max=999, + width="5em", + debounce=False, + ), + "mm", + ), + Block("Season", + Number( + id="start_day", + default=1, + min=1, + max=31, + width="5em", + debounce=False, + ), + Month(id="start_month", default="Jan"), + "-", + Number( + id="end_day", + default=31, + min=1, + max=31, + width="5em", + debounce=False, + ), + Month(id="end_month", default="Mar"), + ), + Block("Projected Years", + Number( + id="start_year", + default=2015, + min=2015, + max=2099, + width="5em", + debounce=False, + ), + "-", + Number( + id="end_year", + default=2019, + min=2015, + max=2099, + width="5em", + debounce=False, + ), + ), + Block("Reference Years", + Number( + id="start_year_ref", + default=1981, + min=1951, + max=2014, + width="5em", + debounce=False, + ), + "-", + Number( + id="end_year_ref", + default=2010, + min=1951, + max=2014, + width="5em", + debounce=False, + ), + ), + button_id="submit_controls", + ), + ), + + lou.description( + "Weather-within-Climate Change Analysis", + """ + This Maproom displays seasonal projected change of key + weather-within-climate agronomic and climatic variables with respect to + historical records. + """, + html.P( + """ + Use the controls in the top banner to choose other variables, models, + scenarios, seasons, projected years and reference to compare with. + """ + ), + html.P( + """ + Click the map (or enter coordinates) to show historical seasonal time + series for this variable of this model, followed by a plume of + possible projected scenarios. + """ + ), + html.P( + """ + Change is expressed as the difference between average over projected + years and average over reference historical years (in the variables + units). + """ + ), + html.H4(["Definitions"]), + html.P([html.B("Mean Max/Min Temperature (mean_Tmax/min):"), """ + Average maximum/minimum temperature in season + """]), + html.P([html.B("Count of Dry/Wet Days (dry/wet_days):"), """ + Number of dry/wet days in the season. A dry/wet day is defined as + lesser or equal / greather than a user-defined threshold. + """]), + html.P([html.B("Count of Frost Days (frost_days):"), """ + Number of frost days in the season. A frost day is defined as + lesser or equal than a user-defined threshold. + """]), + html.P([ + html.B("Max/Min Temperature 90/10th %-ile (Tmax/min_90/10):"),""" + Maximum/Minimum temperature 90/10th percentile in the season. + Obtained through parametric Normal distributions. + """ + ]), + ), + + lou.map(GLOBAL_CONFIG["zoom"]), + + lou.local_single_tabbed( + "Local History and Projections", download_button=True + ), + ) diff --git a/pepsico/proj_wwc/maproom.py b/pepsico/proj_wwc/maproom.py new file mode 100644 index 00000000..2597851b --- /dev/null +++ b/pepsico/proj_wwc/maproom.py @@ -0,0 +1,636 @@ +import dash +import dash_bootstrap_components as dbc +from dash.dependencies import Output, Input, State +import pingrid +from pingrid import CMAPS +from . import layout +import plotly.graph_objects as pgo +import xarray as xr +import pandas as pd +from dateutil.relativedelta import * +from globals_ import FLASK, GLOBAL_CONFIG +import app_calc as ac +import maproom_utilities as mru +import numpy as np + + +def register(FLASK, config): + PFX = f"{GLOBAL_CONFIG['url_path_prefix']}/{config['core_path']}" + TILE_PFX = f"{PFX}/tile" + + # App + + APP = dash.Dash( + __name__, + server=FLASK, + external_stylesheets=[ + dbc.themes.BOOTSTRAP, + ], + url_base_pathname=f"{PFX}/", + meta_tags=[ + {"name": "description", "content": "WWC-CCA"}, + {"name": "viewport", "content": "width=device-width, initial-scale=1.0"}, + ], + ) + APP.title = "WWC-CCA" + + APP.layout = layout.app_layout() + + + @APP.callback( + Output("lat_input", "min"), + Output("lat_input", "max"), + Output("lat_input_tooltip", "children"), + Output("lng_input", "min"), + Output("lng_input", "max"), + Output("lng_input_tooltip", "children"), + Output("map", "center"), + Output("map", "zoom"), + Input("region", "value"), + Input("location", "pathname"), + ) + def initialize(region, path): + scenario = "ssp126" + model = "GFDL-ESM4" + variable = "pr" + data = ac.read_data(scenario, model, variable, region, time_res="daily") + zoom = {"US-CA": 4, "Thailand": 5} + return mru.initialize_map(data) + (zoom[region],) + + + @APP.callback( + Output("loc_marker", "position"), + Output("lat_input", "value"), + Output("lng_input", "value"), + Input("submit_lat_lng","n_clicks"), + Input("map", "click_lat_lng"), + Input("region", "value"), + State("lat_input", "value"), + State("lng_input", "value"), + ) + def pick_location(n_clicks, click_lat_lng, region, latitude, longitude): + # Reading + scenario = "ssp126" + model = "GFDL-ESM4" + variable = "pr" + data = ac.read_data(scenario, model, variable, region, time_res="daily") + initialization_cases = ["region"] + return mru.picked_location( + data, initialization_cases, click_lat_lng, latitude, longitude + ) + + + def select_var(variable): + if variable in [ + "killing_frost", + "mean_Tmin", + "Tmin_10", + "frost_season_length", + "frost_days", + ]: + data_var = "tasmin" + if variable in [ + # "warm_nights", + "mean_Tmax", + "Tmax_90", + # "heatwave_duration", + ]: + data_var = "tasmax" + if variable in [ + # "rain_events", + "dry_days", + "wet_days", + # "longest_dry_spell", + # "longest_wet_spell", + "wet_day_persistence", + "dry_day_persistence", + "dry_spells_mean_length", + "dry_spells_median_length", + ]: + data_var = "pr" + return data_var + + + def local_data( + lat, lng, region, + model, variable, + start_day, start_month, end_day, end_month, + frost_threshold, wet_threshold, + ): + model = [model] if model != "Multi-Model-Average" else [ + "GFDL-ESM4", "IPSL-CM6A-LR", "MPI-ESM1-2-HR","MRI-ESM2-0", "UKESM1-0-LL" + ] + data_var = select_var(variable) + data_ds = xr.concat([xr.Dataset({ + "histo" : ac.read_data( + "historical", m, data_var, region, + time_res="daily", unit_convert=True, + ), + "picontrol" : ac.read_data( + "picontrol", m, data_var, region, + time_res="daily", unit_convert=True, + ), + "ssp126" : ac.read_data( + "ssp126", m, data_var, region, + time_res="daily", unit_convert=True, + ), + "ssp370" : ac.read_data( + "ssp370", m, data_var, region, + time_res="daily", unit_convert=True, + ), + "ssp585" : ac.read_data( + "ssp585", m, data_var, region, + time_res="daily", unit_convert=True, + ), + }) for m in model], "M").assign_coords({"M": [m for m in model]}) + error_msg = None + missing_ds = xr.Dataset() + if any([var is None for var in data_ds.data_vars.values()]): + #This is not supposed to happen: + #would mean something happened to that data + data_ds = missing_ds + error_msg="Data missing for this model or variable" + try: + data_ds = pingrid.sel_snap(data_ds, lat, lng) + except KeyError: + data_ds = missing_ds + error_msg="Grid box out of data domain" + if error_msg == None : + data_ds = ac.seasonal_wwc( + ac.daily_tobegroupedby_season( + data_ds, start_day, start_month, end_day, end_month, + ), + variable, frost_threshold, wet_threshold, + ) + return data_ds, error_msg + + + def plot_ts(ts, name, color, start_format, units): + return pgo.Scatter( + x=ts["T"].dt.strftime(mru.STD_TIME_FORMAT), + y=ts.values, + customdata=ts["seasons_ends"].dt.strftime("%d %B %Y"), + hovertemplate=("%{x|"+start_format+"}%{customdata}: %{y:.2f} " + units), + name=name, + line=pgo.scatter.Line(color=color), + connectgaps=False, + ) + + + def add_period_shape( + graph, data, start_year, end_year, fill_color, line_color, annotation + ): + return graph.add_vrect( + x0=data["T"].where( + lambda x : (x.dt.year == int(start_year)), drop=True + ).dt.strftime(mru.STD_TIME_FORMAT).values[0], + #it's hard to believe this is how it is done + x1=( + pd.to_datetime(data["seasons_ends"].where( + lambda x : (x.dt.year == int(end_year)), drop=True + ).dt.strftime(mru.STD_TIME_FORMAT).values[0] + ) + relativedelta(months=+1)).strftime(mru.STD_TIME_FORMAT), + fillcolor=fill_color, opacity=0.2, + line_color=line_color, line_width=3, + layer="below", + annotation_text=annotation, annotation_position="top left", + #editable=True, #a reminder it might be the way to interact + ) + + + @APP.callback( + Output("btn_csv", "disabled"), + Input("lat_input", "value"), + Input("lng_input", "value"), + Input("lat_input", "min"), + Input("lng_input", "min"), + Input("lat_input", "max"), + Input("lng_input", "max"), + ) + def invalid_button(lat, lng, lat_min, lng_min, lat_max, lng_max): + return ( + lat < float(lat_min) or lat > float(lat_max) + or lng < float(lng_min) or lng > float(lng_max) + ) + + + @APP.callback( + Output("download-dataframe-csv", "data"), + Input("btn_csv", "n_clicks"), + State("loc_marker", "position"), + State("region", "value"), + State("variable", "value"), + State("start_day", "value"), + State("start_month", "value"), + State("end_day", "value"), + State("end_month", "value"), + State("frost", "value"), + State("wet", "value"), + prevent_initial_call=True, + ) + def send_data_as_csv( + n_clicks, marker_pos, region, variable, + start_day, start_month, end_day, end_month, + frost_threshold, wet_threshold, + ): + lat = marker_pos[0] + lng = marker_pos[1] + start_day = int(start_day) + end_day = int(end_day) + start_month = ac.strftimeb2int(start_month) + end_month = ac.strftimeb2int(end_month) + frost_threshold = float(frost_threshold) + wet_threshold = float(wet_threshold) + model = "Multi-Model-Average" + data_ds, error_msg = local_data( + lat, lng, region, model, variable, + start_day, start_month, end_day, end_month, + frost_threshold, wet_threshold, + ) + if error_msg == None : + lng_units = "E" if (lng >= 0) else "W" + lat_units = "N" if (lat >= 0) else "S" + file_name = ( + f'{data_ds["histo"]["T"].dt.strftime("%m%d")[0].values}-' + f'{data_ds["histo"]["seasons_ends"].dt.strftime("%m%d")[0].values}' + f'_{variable}_{abs(lat)}{lat_units}_{abs(lng)}{lng_units}' + f'.csv' + ) + df = data_ds.to_dataframe() + return dash.dcc.send_data_frame(df.to_csv, file_name) + else : + return None + + + @APP.callback( + Output("local_graph", "figure"), + Input("loc_marker", "position"), + Input("region", "value"), + Input("submit_controls","n_clicks"), + State("model", "value"), + State("variable", "value"), + State("start_day", "value"), + State("end_day", "value"), + State("start_month", "value"), + State("end_month", "value"), + State("start_year", "value"), + State("end_year", "value"), + State("start_year_ref", "value"), + State("end_year_ref", "value"), + State("frost", "value"), + State("wet", "value"), + ) + def local_plots( + marker_pos, region, n_clicks, model, variable, + start_day, end_day, start_month, end_month, + start_year, end_year, start_year_ref, end_year_ref, + frost_threshold, wet_threshold, + ): + lat = marker_pos[0] + lng = marker_pos[1] + start_day = int(start_day) + end_day = int(end_day) + start_month = ac.strftimeb2int(start_month) + end_month = ac.strftimeb2int(end_month) + frost_threshold = float(frost_threshold) + wet_threshold = float(wet_threshold) + data_ds, error_msg = local_data( + lat, lng, region, model, variable, + start_day, start_month, end_day, end_month, + frost_threshold, wet_threshold, + ) + if error_msg != None : + local_graph = pingrid.error_fig(error_msg) + else : + if (end_month < start_month) : + start_format = "%d %b %Y - " + else: + start_format = "%d %b-" + local_graph = pgo.Figure() + data_color = { + "histo": "blue", "picontrol": "green", + "ssp126": "yellow", "ssp370": "orange", "ssp585": "red", + } + lng_units = "˚E" if (lng >= 0) else "˚W" + lat_units = "˚N" if (lat >= 0) else "˚S" + for var in data_ds.data_vars: + local_graph.add_trace(plot_ts( + data_ds[var].mean("M", keep_attrs=True), var, data_color[var], + start_format, data_ds[var].attrs["units"] + )) + add_period_shape( + local_graph, + data_ds, + start_year_ref, + end_year_ref, + "blue", + "RoyalBlue", + "reference period", + ) + add_period_shape( + local_graph, + data_ds, + start_year, + end_year, + "LightPink", + "Crimson", + "projected period", + ) + local_graph.update_layout( + xaxis_title="Time", + yaxis_title=( + # Compared to the monthly case, I presume that groupby dropped + # the long_name, which is fine since it's indeed transformed. Can + # consider resetting it at some point in the calculation + #f'{data_ds["histo"].attrs["long_name"]} ' + f'({data_ds["histo"].attrs["units"]})' + ), + title={ + "text": ( + f'{data_ds["histo"]["T"].dt.strftime("%d %b")[0].values}-' + f'{data_ds["histo"]["seasons_ends"].dt.strftime("%d %b")[0].values}' + f' {variable} seasonal average from model {model} ' + f'at ({abs(lat)}{lat_units}, {abs(lng)}{lng_units})' + ), + "font": dict(size=14), + }, + margin=dict(l=30, r=30, t=30, b=30), + ) + return local_graph + + + @APP.callback( + Output("map_description", "children"), + Input("submit_controls", "n_clicks"), + State("scenario", "value"), + State("model", "value"), + State("variable", "value"), + State("start_day", "value"), + State("end_day", "value"), + State("start_month", "value"), + State("end_month", "value"), + State("start_year", "value"), + State("end_year", "value"), + State("start_year_ref", "value"), + State("end_year_ref", "value"), + ) + def write_map_description( + n_clicks, scenario, model, variable, + start_day, end_day, start_month, end_month, + start_year, end_year, start_year_ref, end_year_ref, + ): + return ( + f'The Map displays the change in ' + f'({start_day} {start_month} - {end_day} {end_month}) season of ' + f'{variable} from {model} model under {scenario} scenario projected for ' + f'{start_year}-{end_year} with respect to historical {start_year_ref}-' + f'{end_year_ref}' + ) + + + @APP.callback( + Output("map_title", "children"), + Input("submit_controls","n_clicks"), + State("scenario", "value"), + State("model", "value"), + State("variable", "value"), + State("start_day", "value"), + State("end_day", "value"), + State("start_month", "value"), + State("end_month", "value"), + State("start_year", "value"), + State("end_year", "value"), + State("start_year_ref", "value"), + State("end_year_ref", "value"), + ) + def write_map_title( + n_clicks, scenario, model, variable, + start_day, end_day, start_month, end_month, + start_year, end_year, start_year_ref, end_year_ref, + ): + return ( + f'({start_day} {start_month} - {end_day} {end_month}) ' + f'{start_year}-{end_year} ' + f'{scenario} {model} {variable} change from ' + f'{start_year_ref}-{end_year_ref}' + ) + + + APP.clientside_callback( + """function(start_year, end_year, start_year_ref, end_year_ref) { + if (start_year && end_year) { + invalid_start_year = (start_year > end_year) + invalid_end_year = invalid_start_year + } else { + invalid_start_year = !start_year + invalid_end_year = !end_year + } + if (start_year_ref && end_year_ref) { + invalid_start_year_ref = (start_year_ref > end_year_ref) + invalid_end_year_ref = invalid_start_year_ref + } else { + invalid_start_year_ref = !start_year_ref + invalid_end_year_ref = !end_year_ref + } + return [ + invalid_start_year, invalid_end_year, + invalid_start_year_ref, invalid_end_year_ref, + ( + invalid_start_year || invalid_end_year + || invalid_start_year_ref || invalid_end_year_ref + ), + ] + } + """, + Output("start_year", "invalid"), + Output("end_year", "invalid"), + Output("start_year_ref", "invalid"), + Output("end_year_ref", "invalid"), + Output("submit_controls", "disabled"), + Input("start_year", "value"), + Input("end_year", "value"), + Input("start_year_ref", "value"), + Input("end_year_ref", "value"), + ) + + + def seasonal_change( + scenario, model, variable, region, + start_day, end_day, start_month, end_month, + start_year, end_year, start_year_ref, end_year_ref, + frost_threshold, wet_threshold, + ): + model = [model] if model != "Multi-Model-Average" else [ + "GFDL-ESM4", "IPSL-CM6A-LR", "MPI-ESM1-2-HR","MRI-ESM2-0", "UKESM1-0-LL" + ] + data_var = select_var(variable) + ref = xr.concat([ + ac.seasonal_wwc( + ac.daily_tobegroupedby_season( + ac.read_data( + "historical", m, data_var, region, + time_res="daily", unit_convert=True + ).sel( + T=slice(str(start_year_ref), str(end_year_ref)) + ).to_dataset(), + start_day, start_month, end_day, end_month, + ), + variable, frost_threshold, wet_threshold, + ).mean(dim="T", keep_attrs=True) for m in model + ], "M").mean("M", keep_attrs=True) + data = xr.concat([ + ac.seasonal_wwc( + ac.daily_tobegroupedby_season( + ac.read_data( + scenario, m, data_var, region, + time_res="daily", unit_convert=True + ).sel( + T=slice(str(start_year), str(end_year)) + ).to_dataset(), + start_day, start_month, end_day, end_month, + ), + variable, frost_threshold, wet_threshold, + ).mean(dim="T", keep_attrs=True) for m in model + ], "M").mean("M", keep_attrs=True) + #Tedious way to make a subtraction only to keep attributes (units) + return xr.apply_ufunc( + np.subtract, data, ref, dask="allowed", keep_attrs="drop_conflicts", + ).rename({"X": "lon", "Y": "lat"}) + + + def map_attributes(data, variable): + data = data[select_var(variable)] + if data.name == "pr": + colorscale = CMAPS["prcp_anomaly"] + if data.name in ["tasmin", "tasmax"]: + colorscale = CMAPS["temp_anomaly"] + if variable in ["frost_days", "dry_days"]: + colorscale = colorscale.reversed() + map_amp = np.max(np.abs(data)).values + colorscale = colorscale.rescaled(-1*map_amp, map_amp) + return colorscale, colorscale.scale[0], colorscale.scale[-1] + + + @APP.callback( + Output("colorbar", "colorscale"), + Output("colorbar", "min"), + Output("colorbar", "max"), + Output("colorbar", "unit"), + Input("region", "value"), + Input("submit_controls","n_clicks"), + State("scenario", "value"), + State("model", "value"), + State("variable", "value"), + State("start_day", "value"), + State("end_day", "value"), + State("start_month", "value"), + State("end_month", "value"), + State("start_year", "value"), + State("end_year", "value"), + State("start_year_ref", "value"), + State("end_year_ref", "value"), + State("frost", "value"), + State("wet", "value"), + ) + def draw_colorbar( + region, n_clicks, scenario, model, variable, + start_day, end_day, start_month, end_month, + start_year, end_year, start_year_ref, end_year_ref, + frost_threshold, wet_threshold, + ): + start_day = int(start_day) + end_day = int(end_day) + start_month = ac.strftimeb2int(start_month) + end_month = ac.strftimeb2int(end_month) + frost_threshold = float(frost_threshold) + wet_threshold = float(wet_threshold) + data = seasonal_change( + scenario, model, variable, region, + start_day, end_day, start_month, end_month, + start_year, end_year, start_year_ref, end_year_ref, + frost_threshold, wet_threshold, + ) + colorbar, min, max = map_attributes(data, variable) + return ( + colorbar.to_dash_leaflet(), min, max, + data[select_var(variable)].attrs["units"] + ) + + + @APP.callback( + Output("layers_control", "children"), + Output("map_warning", "is_open"), + Input("region", "value"), + Input("submit_controls", "n_clicks"), + State("scenario", "value"), + State("model", "value"), + State("variable", "value"), + State("start_day", "value"), + State("end_day", "value"), + State("start_month", "value"), + State("end_month", "value"), + State("start_year", "value"), + State("end_year", "value"), + State("start_year_ref", "value"), + State("end_year_ref", "value"), + State("frost", "value"), + State("wet", "value"), + ) + def make_map( + region, n_clicks, scenario, model, variable, + start_day, end_day, start_month, end_month, + start_year, end_year, start_year_ref, end_year_ref, + frost_threshold, wet_threshold, + ): + try: + send_alarm = False + url_str = ( + f"{TILE_PFX}/{{z}}/{{x}}/{{y}}/{region}/{scenario}/{model}/" + f"{variable}/{start_day}/{end_day}/{start_month}/{end_month}/" + f"{start_year}/{end_year}/{start_year_ref}/{end_year_ref}/" + f"{frost_threshold}/{wet_threshold}" + ) + except: + url_str= "" + send_alarm = True + return mru.layers_controls( + url_str, f"wwc_change_{region}", "WWC Change", + GLOBAL_CONFIG["datasets"][f"shapes_adm_{region}"], GLOBAL_CONFIG, + adm_id_suffix=region, + ), send_alarm + + + @FLASK.route( + ( + f"{TILE_PFX}///////" + f"/////" + f"////" + f"//" + ), + endpoint=f"{config['core_path']}" + ) + def fcst_tiles(tz, tx, ty, + region, scenario, model, variable, + start_day, end_day, start_month, end_month, + start_year, end_year, start_year_ref, end_year_ref, + frost_threshold, wet_threshold, + ): + start_day = int(start_day) + end_day = int(end_day) + start_month = ac.strftimeb2int(start_month) + end_month = ac.strftimeb2int(end_month) + frost_threshold = float(frost_threshold) + wet_threshold = float(wet_threshold) + data = seasonal_change( + scenario, model, variable, region, + start_day, end_day, start_month, end_month, + start_year, end_year, start_year_ref, end_year_ref, + frost_threshold, wet_threshold, + ) + ( + data[select_var(variable)].attrs["colormap"], + data[select_var(variable)].attrs["scale_min"], + data[select_var(variable)].attrs["scale_max"], + ) = map_attributes(data, variable) + resp = pingrid.tile(data[select_var(variable)], tx, ty, tz) + return resp diff --git a/pingrid/impl.py b/pingrid/impl.py index 20af3351..5fbe5dbb 100644 --- a/pingrid/impl.py +++ b/pingrid/impl.py @@ -883,8 +883,10 @@ def produce_shape_tile( _TEMP_ANOMALY_CS = ColorScale( "temp_anomaly", - [PURPLE, CYAN, WHITE, WHITE, YELLOW, RED, FIREBRICK], - [-10, -1, -1, 1, 1, 10, 10], + # The duplicate first anchor makes the scale symmetric and thus allows to use the + # reversed method + [PURPLE, PURPLE, CYAN, WHITE, WHITE, YELLOW, RED, FIREBRICK], + [-10, -10, -1, -1, 1, 1, 10, 10], ) _RAIN_PNE_CS = _RAIN_POE_CS.reversed(name="rain_pne")