diff --git a/ccsr/config-defaults.yaml b/ccsr/config-defaults.yaml index 1705a01a..a25e2b7a 100644 --- a/ccsr/config-defaults.yaml +++ b/ccsr/config-defaults.yaml @@ -14,11 +14,20 @@ maprooms: flex_fcst: # App set up - - title: flex_cst_precip + - title: flex_fcst_precip core_path: flex_fcst_precip forecast_path: /Data/data21/iri/FD.3 variable: Precipitation y_transform: False method: elr # or pycpt var: precip - issue_format: "%b %Y" \ No newline at end of file + issue_format: "%b %Y" + + terc_fcst: + - title: terc_fcst_precip + core_path: terc_fcst_precip + forecast_path: /Data/data21/iri/FD/NMME_Seasonal_Forecast/Precipitation_ELR + variable: Precipitation + var: precip + issue_format: "%b %Y" + \ No newline at end of file diff --git a/ccsr/config-iri.yaml b/ccsr/config-iri.yaml index 990209e6..73765c3e 100644 --- a/ccsr/config-iri.yaml +++ b/ccsr/config-iri.yaml @@ -35,3 +35,11 @@ maprooms: method: elr # or pycpt var: precip issue_format: "%b %Y" + + terc_fcst: + - title: terc_fcst_precip + core_path: terc_fcst_precip + forecast_path: /Data/data21/iri/FD/NMME_Seasonal_Forecast/Precipitation_ELR + variable: Precipitation + var: precip + issue_format: "%b %Y" diff --git a/ccsr/flex_fcst/layout.py b/ccsr/flex_fcst/layout.py index 88a6b9e7..f7e9a7f3 100644 --- a/ccsr/flex_fcst/layout.py +++ b/ccsr/flex_fcst/layout.py @@ -1,8 +1,5 @@ -from dash import dcc from dash import html -import dash_bootstrap_components as dbc -import dash_leaflet as dlf -from fieldsets import Block, Select, PickPoint, Month, Number +from fieldsets import Block, Select, PickPoint, Number import layout_utilities as lou from globals_ import GLOBAL_CONFIG diff --git a/ccsr/flex_fcst_calc.py b/ccsr/flex_fcst_calc.py index 4d07b0fa..a1d3a9f9 100644 --- a/ccsr/flex_fcst_calc.py +++ b/ccsr/flex_fcst_calc.py @@ -4,6 +4,7 @@ import pingrid import numpy as np import datetime +import predictions def get_targets_dict(fcst_conf, fcst_ds, start_date): @@ -11,7 +12,7 @@ def get_targets_dict(fcst_conf, fcst_ds, start_date): return pycpt.pycptv2_targets_dict(fcst_ds, start_date=start_date) else: S = datetime.datetime( - int(start_date[4:8]), ife.strftimeb2int(start_date[0:3]), 1 + int(start_date[4:8]), predictions.strftimeb2int(start_date[0:3]), 1 ) return ife.targets_dict(fcst_ds, S) diff --git a/ccsr/iri_fd_elr.py b/ccsr/iri_fd_elr.py index 97d7ba9d..3a419696 100644 --- a/ccsr/iri_fd_elr.py +++ b/ccsr/iri_fd_elr.py @@ -1,4 +1,3 @@ -from pathlib import Path import xarray as xr import datetime import predictions @@ -6,42 +5,6 @@ import numpy as np -def strftimeb2int(strftimeb): - """Convert month values to integers (1-12) from strings. - - Parameters - ---------- - strftimeb : str - String value representing months of year. - Returns - ------- - strftimebint : int - Integer value corresponding to month. - See Also - -------- - Notes - ----- - Examples - -------- - """ - strftimeb_all = { - "Jan": 1, - "Feb": 2, - "Mar": 3, - "Apr": 4, - "May": 5, - "Jun": 6, - "Jul": 7, - "Aug": 8, - "Sep": 9, - "Oct": 10, - "Nov": 11, - "Dec": 12, - } - strftimebint = strftimeb_all[strftimeb] - return strftimebint - - def get_elr_S(data_path, var): """Returns list of dates of forecast issues. @@ -58,7 +21,7 @@ def get_elr_S(data_path, var): files_list = data_path.glob(f"forecast_mean_{var}_*.nc") return sorted([datetime.datetime( int(f.name[-7:-3]), - strftimeb2int(f.name[-10:-7]), + predictions.strftimeb2int(f.name[-10:-7]), 1, ) for f in files_list], reverse=True) diff --git a/ccsr/layout_utilities.py b/ccsr/layout_utilities.py index 776e4238..ab89fd92 100644 --- a/ccsr/layout_utilities.py +++ b/ccsr/layout_utilities.py @@ -103,8 +103,10 @@ def description(title, subtitle, *elems): def map( default_zoom, layers_control_position="topleft", scale_control_position="bottomright", - cb_nTicks=9, cb_opacity=1, cb_tooltip=True, - cb_position="topright", cb_width=10, cb_height=300, + colorbars={"colorbar" : { + "nTicks": 9, "opacity": 1, "tooltip": True, "position": "topright", + "width": 10, "height": 300 + }}, ): """ A dlf map topped with and H5 title, @@ -134,21 +136,18 @@ def map( dlf.ScaleControl( imperial=False, position=scale_control_position ), + ] + [ dlf.Colorbar( - id="colorbar", - nTicks=cb_nTicks, - opacity=cb_opacity, - tooltip=cb_tooltip, - position=cb_position, - width=cb_width, - height=cb_height, + id=cb, className="p-1", style={ "background": "white", "border-style": "inset", "-moz-border-radius": "4px", "border-radius": "4px", "border-color": "LightGrey", }, - ), + **cbd, + ) + for cb, cbd in colorbars.items() ], id="map", center=None, diff --git a/ccsr/predictions.py b/ccsr/predictions.py index 862aa9ad..c3c131f3 100644 --- a/ccsr/predictions.py +++ b/ccsr/predictions.py @@ -52,3 +52,38 @@ def target_range_formatting(target_start, target_end, time_units): date_range = f"{target_start_str}-{target_end_str}" return date_range + +def strftimeb2int(strftimeb): + """Convert month values to integers (1-12) from strings. + + Parameters + ---------- + strftimeb : str + String value representing months of year. + Returns + ------- + strftimebint : int + Integer value corresponding to month. + See Also + -------- + Notes + ----- + Examples + -------- + """ + strftimeb_all = { + "Jan": 1, + "Feb": 2, + "Mar": 3, + "Apr": 4, + "May": 5, + "Jun": 6, + "Jul": 7, + "Aug": 8, + "Sep": 9, + "Oct": 10, + "Nov": 11, + "Dec": 12, + } + strftimebint = strftimeb_all[strftimeb] + return strftimebint diff --git a/ccsr/terc_fcst/__init__.py b/ccsr/terc_fcst/__init__.py new file mode 100644 index 00000000..0379666a --- /dev/null +++ b/ccsr/terc_fcst/__init__.py @@ -0,0 +1 @@ +from .maproom import register diff --git a/ccsr/terc_fcst/layout.py b/ccsr/terc_fcst/layout.py new file mode 100644 index 00000000..5b7a6369 --- /dev/null +++ b/ccsr/terc_fcst/layout.py @@ -0,0 +1,56 @@ +from dash import html +from fieldsets import Block, PickPoint +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("Forecast", + Block("issue_block", "Issue", html.Div(id="start_div")), + Block("target_block", "Target Period", html.Div(id="lead_div")), + PickPoint(width="8em"), + ), + + lou.description( + "Precipitation Terciles Seasonal Forecast", + """ + The seasonal forecast for above-, below- and near-normal precipitation + from the IRI. + """, + html.P( + """ + The default map shows globally the seasonal precipitation forecast + tercile probability. The historical climatology used is 1982-2010 up + to the forecast issued in August 2021, and is 1991-2020 from the + forecast issued in September 2021. The forecast shown is the latest + forecast made (e.g. Dec 2017) for the next season to come + (e.g. Jan-Mar 2018). Four different seasons are forecasted and it is + also possible to consult forecasts made previously. The forecasts + are directly computed from the extended logistic regression model as + tercile probabilities. + """ + ), + ), + + lou.map(GLOBAL_CONFIG["zoom"], scale_control_position="bottomleft", colorbars={ + "below_cb" : { + "nTicks": 6, "min": 37.5, "max": 100, "tickDecimals": 1, "unit": "%", + "opacity": 1, "tooltip": True, "position": "bottomright", + "width": 10, "height": 140, + }, + "above_cb" : { + "nTicks": 6, "min": 37.5, "max": 100, "tickDecimals": 1, "unit": "%", + "opacity": 1, "tooltip": True, "position": "topright", + "width": 10, "height": 140, + }, + }), + + lou.local_single_tabbed("Tercile Probabilities"), + ) diff --git a/ccsr/terc_fcst/maproom.py b/ccsr/terc_fcst/maproom.py new file mode 100644 index 00000000..77d74979 --- /dev/null +++ b/ccsr/terc_fcst/maproom.py @@ -0,0 +1,262 @@ +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.express as px +import numpy as np +import xarray as xr +import terc_fcst_calc as tfc +import maproom_utilities as mru +from globals_ import GLOBAL_CONFIG +from fieldsets import Select +import datetime +import predictions + + +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": "Forecast"}, + {"name": "viewport", "content": "width=device-width, initial-scale=1.0"}, + ], + ) + APP.title = "Tercile Forecast" + + APP.layout = layout.app_layout() + + @APP.callback( + Output("start_div", "children"), + Output("lead_div", "children"), + Output("below_cb", "colorscale"), + Output("above_cb", "colorscale"), + 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"), + Input("location", "pathname"), + ) + def initialize(path): + S_dates = tfc.get_issues(config) + issue_select = Select("start_date", options=[d.strftime("%b %Y") for d in S_dates]) + fcst_ds = tfc.get_fcst(config) + target_dict = tfc.get_targets_dict(fcst_ds, S_dates[0]) + lead_select = Select("lead_time", + options=[ld["value"] for ld in target_dict], + labels=[ld["label"] for ld in target_dict], + ) + return ( + issue_select, lead_select, + # while the data is classified and mapped according to classification + # colormap, I draw the colorscale from 2 colormaps (for below and above) + # corresponding to those classifications and with the probability values + CMAPS["prcp_terciles_below"].to_dash_leaflet(), + CMAPS["prcp_terciles_above"].to_dash_leaflet(), + ) + mru.initialize_map(fcst_ds) + + + @APP.callback( + Output("lead_time","options"), + Input("start_date","value"), + ) + def target_range_options(start_date): + S = datetime.datetime( + int(start_date[4:8]), predictions.strftimeb2int(start_date[0:3]), 1 + ) + return tfc.get_targets_dict(tfc.get_fcst(config, start_date), S) + + + @APP.callback( + Output("map_title", "children"), + Input("start_date", "value"), + Input("lead_time", "value"), + Input("lead_time", "options"), + ) + def write_map_title(start_date, lead_time, targets): + for option in targets : + if option["value"] == int(lead_time) : + target = option["label"] + return f'{target} {config["variable"]} Forecast issued {start_date}' + + + @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"), + State("lat_input", "value"), + State("lng_input", "value") + ) + def pick_location(n_clicks, click_lat_lng, latitude, longitude): + return mru.picked_location( + tfc.get_fcst(config), [""], click_lat_lng, latitude, longitude + ) + + + @APP.callback( + Output("local_graph", "figure"), + Input("loc_marker", "position"), + Input("start_date", "value"), + Input("lead_time", "value"), + Input("lead_time", "options"), + ) + def local_plots(marker_pos, start_date, lead_time, targets): + lat = marker_pos[0] + lng = marker_pos[1] + fcst_ds = tfc.get_fcst(config, start_date=start_date, lead_time=lead_time) + # Errors handling + error_msg = None + if fcst_ds["prob"] is None: + error_msg = "Data missing for this issue and target" + else: + try: + fcst_ds = pingrid.sel_snap(fcst_ds, lat, lng) + isnan = np.isnan(fcst_ds["prob"]).all() + if isnan: + error_msg = "Data missing at this location" + except KeyError: + error_msg = "Grid box out of data domain" + if error_msg is not None: + local_graph = pingrid.error_fig(error_msg) + else: + for option in targets : + if option["value"] == int(lead_time) : + target = option["label"] + lng_units = "˚E" if (fcst_ds['X'] >= 0) else "˚W" + lat_units = "˚N" if (fcst_ds['Y'] >= 0) else "˚S" + #Get colors from the colormap to assign each bar + if fcst_ds["prob"].isel(cat=0) >= 37.5: + below_color = CMAPS["prcp_terciles_below"].to_rgba_array()[ + int( + 0.5 + + 255 + * (fcst_ds["prob"].isel(cat=0).values - 37.5) + / (100 - 37.5) + ), + 0:3, + ] + below_color = ( + f"rgb({below_color[0]}, {below_color[1]}, {below_color[2]})" + ) + else: + below_color = "lightyellow" + if fcst_ds["prob"].isel(cat=2) >= 37.5: + above_color = CMAPS["prcp_terciles_above"].to_rgba_array()[ + int( + 0.5 + + 255 + * (fcst_ds["prob"].isel(cat=2).values - 37.5) + / (100 - 37.5) + ), + 0:3, + ] + above_color = ( + f"rgb({above_color[0]}, {above_color[1]}, {above_color[2]})" + ) + else: + above_color = "rgb(220, 255, 220)" + local_graph = px.bar( + ( + fcst_ds["prob"] + .assign_coords(Terciles=("cat", ["Below", "Normal", "Above"])) + .to_dataframe() + ), + x="Terciles", y="prob", range_y=[0, 100], color="Terciles", + color_discrete_sequence=[below_color, "grey", above_color], + title=( + f"{target} Forecast issued {start_date} at" + f" ({abs(fcst_ds['Y'].values)}{lat_units}" + f", {abs(fcst_ds['X'].values)}{lng_units})" + ) + ).add_hline( + y=100/3, line_dash="dash", line_color="black", + annotation_text="climatology", annotation_position="bottom right", + ) + y_ticks = [0, 37.5, 42.5, 47.5, 57.5, 67.5, 100] + local_graph.update_layout( + yaxis={ + "tickmode":"array", "tickvals": y_ticks, + "ticktext": [f"{yt} %" for yt in y_ticks]} + ) + return local_graph + + + @APP.callback( + Output("layers_control", "children"), + Input("start_date","value"), + Input("lead_time","value") + ) + def make_map(start_date, lead_time): + url_str = f"{TILE_PFX}/{{z}}/{{x}}/{{y}}/{start_date}/{lead_time}" + return mru.layers_controls( + url_str, "terc_fcst", "Forecast", + GLOBAL_CONFIG["datasets"]["shapes_adm"], GLOBAL_CONFIG, + ) + + @FLASK.route( + f"{TILE_PFX}/////", + endpoint=f"{config['core_path']}" + ) + def fcst_tiles(tz, tx, ty, start_date, lead_time): + # Reading + fcst_ds = tfc.get_fcst(config, start_date=start_date, lead_time=lead_time) + dominant_cat = fcst_ds["prob"].idxmax(dim="cat") + # Classifying each tercile forecast from 0 to 5 + fcst_class = ( + (fcst_ds["prob"] > 37.5) * 1. + + (fcst_ds["prob"] > 42.5) * 1. + + (fcst_ds["prob"] > 47.5) * 1. + + (fcst_ds["prob"] > 57.5) * 1. + + (fcst_ds["prob"] > 67.5) * 1. + ) + # Can I do this without for? + # I am afraid not because xr.where would not know how to broadcast...? + dominant_fcst_class = sum([ + xr.where( + # Where class is dominant + fcst_ds["category"].isel(cat=c) == dominant_cat, + # Translates class by 5c except for climatology that must remain 0 + # to get: + # Below (1, 2, 3, 4, 5) + # Normal (6, 7, 8, 9, 10) + # Above (11, 12, 13, 14, 15) + xr.where( + fcst_class.isel(cat=c) != 0, fcst_class.isel(cat=c)+ 5 * c, 0 + ), + # Elsewhere set to 0 so can squeeze the list with sum + 0, + ) + for c in range(3) + ]) + dominant_fcst_class = xr.apply_ufunc( + # Return all Normal classes to one class: 6 + # Reclassify Above to (7, 8, 9, 10, 11) + np.piecewise, dominant_fcst_class, + [ + dominant_fcst_class.data < 6, + dominant_fcst_class.data > 10, + ], + [lambda x : x, lambda x : x - 4, 6], + ) + dominant_fcst_class.attrs["colormap"] = CMAPS["prcp_terciles"] + dominant_fcst_class.attrs["scale_min"] = 0 + dominant_fcst_class.attrs["scale_max"] = 11 + return pingrid.tile( + dominant_fcst_class.rename({"X": "lon", "Y": "lat"}), tx, ty, tz + ) diff --git a/ccsr/terc_fcst_calc.py b/ccsr/terc_fcst_calc.py new file mode 100644 index 00000000..6d81393e --- /dev/null +++ b/ccsr/terc_fcst_calc.py @@ -0,0 +1,49 @@ +from pathlib import Path +import datetime +import predictions +from dateutil.relativedelta import relativedelta +import xarray as xr + + +def get_issues(fcst_conf): + files_list = Path(fcst_conf["forecast_path"]).glob( + f"forecast_{fcst_conf['var']}_*.nc" + ) + return sorted([datetime.datetime( + int(f.name[-7:-3]), + predictions.strftimeb2int(f.name[-10:-7]), + 1, + ) for f in files_list], reverse=True) + + +def get_fcst(fcst_conf, start_date=None, lead_time=None): + data_path = Path(fcst_conf["forecast_path"]) + var = fcst_conf["var"] + S = None if start_date is None else (start_date[0:3] + start_date[4:8]) + if S is None: + S = get_issues(fcst_conf)[0].strftime('%b%Y') + fcst_ds = xr.open_dataset(data_path / f"forecast_{var}_{S}.nc") + fcst_ds = ( + fcst_ds + .assign_coords({"cat" : fcst_ds["category"]}) + .assign_coords(longitude=(((fcst_ds.longitude + 180) % 360) - 180)) + .sortby("longitude") + .rename({"longitude": "X", "latitude": "Y"}) + ) + if lead_time is not None: + fcst_ds = fcst_ds.sel(lead=int(lead_time)) + return fcst_ds + + +def get_targets_dict(fcst_ds, S_date): + leads = fcst_ds["lead"] + return [ + { + "label": predictions.target_range_formatting( + S_date + relativedelta(months=(int(lead))), + S_date + relativedelta(months=(int(lead) + 2)), + "months", + ), + "value": lead, + } for lead in leads.values + ] \ No newline at end of file diff --git a/pingrid/impl.py b/pingrid/impl.py index 5fbe5dbb..e1c5a8d4 100644 --- a/pingrid/impl.py +++ b/pingrid/impl.py @@ -33,6 +33,7 @@ 'DEEPSKYBLUE', 'FIREBRICK', 'GREEN', + 'GRAY', 'LIMEGREEN', 'MOCCASIN', 'NAVY', @@ -650,6 +651,7 @@ def produce_shape_tile( DARKRED = Color(128, 0, 0) DEEPSKYBLUE = Color(0, 191, 255) FIREBRICK = Color(178, 34, 34) +GRAY = Color(192, 192, 192) GREEN = Color(0, 255, 0) LIMEGREEN = Color(50, 205, 50) MOCCASIN = Color(255, 228, 181) @@ -834,6 +836,50 @@ def produce_shape_tile( 450, 450, 500], ) +_PRCP_TERCILES_CS = ColorScale( + "prcp_terciles", + [ + WHITE, WHITE, + Color(250, 250, 0), Color(250, 250, 0), + Color(232, 184, 51), Color(232, 184, 51), + Color(208, 128, 51), Color(208, 128, 51), + Color(170, 70, 30), Color(170, 70, 30), + Color(120, 50, 0), Color(120, 50, 0), + GRAY, GRAY, + Color(200, 250, 200), Color(200, 250, 200), + Color(150, 250, 150), Color(150, 250, 150), + Color(90, 190, 100), Color(90, 190, 100), + Color(0, 150, 210), Color(0, 150, 210), + Color(0, 51, 255), Color(0, 51, 255), + ], + [0, 0.5, 0.5, 1.5, 1.5, 2.5, 2.5, 3.5, 3.5, 4.5, 4.5, 5.5, + 5.5, 6.5, 6.5, 7.5, 7.5, 8.5, 8.5, 9.5, 9.5, 10.5, 10.5, 11], +) + +_PRCP_TERCILES_BELOW_CS = ColorScale( + "prcp_terciles_below", + [ + Color(250, 250, 0), Color(250, 250, 0), + Color(232, 184, 51), Color(232, 184, 51), + Color(208, 128, 51), Color(208, 128, 51), + Color(170, 70, 30), Color(170, 70, 30), + Color(120, 50, 0), Color(120, 50, 0), + ], + [37.5, 42.5, 42.5, 47.5, 47.7, 57.5, 57.5, 67.5, 67.5, 100], +) + +_PRCP_TERCILES_ABOVE_CS = ColorScale( + "prcp_terciles_above", + [ + Color(200, 250, 200), Color(200, 250, 200), + Color(150, 250, 150), Color(150, 250, 150), + Color(90, 190, 100), Color(90, 190, 100), + Color(0, 150, 210), Color(0, 150, 210), + Color(0, 51, 255), Color(0, 51, 255), + ], + [37.5, 42.5, 42.5, 47.5, 47.7, 57.5, 57.5, 67.5, 67.5, 100], +) + _RAIN_POE_CS = ColorScale( "rain_poe", [BLACK, BROWN, ORANGE, YELLOW, MOCCASIN, MOCCASIN, LIMEGREEN, TURQUOISE, BLUE, PURPLE], @@ -979,6 +1025,9 @@ def produce_shape_tile( _PRECIP_CS, _PRCP_ANOMALY_CS, _PRCP_ANOMALY_BLUE_CS, + _PRCP_TERCILES_CS, + _PRCP_TERCILES_ABOVE_CS, + _PRCP_TERCILES_BELOW_CS, _RAIN_PNE_CS, _RAIN_POE_CS, _RAINBOW_CS,