From e5562fcbf61678ce16dbadc021be7d25e15f8759 Mon Sep 17 00:00:00 2001 From: remic Date: Fri, 10 Mar 2023 12:24:30 -0500 Subject: [PATCH 1/6] total rain between onset and cessation, multiplying size of problem by munber of years --- enacts/config-defaults.yaml | 26 +++++++++++------------ enacts/config-zmd.yaml | 5 ++++- enacts/onset/maproom.py | 41 +++++++++++++++++++++++++++++-------- 3 files changed, 50 insertions(+), 22 deletions(-) diff --git a/enacts/config-defaults.yaml b/enacts/config-defaults.yaml index 2f8877fbf..a6c030746 100644 --- a/enacts/config-defaults.yaml +++ b/enacts/config-defaults.yaml @@ -66,19 +66,19 @@ onset: menu_label: How likely is the rainy season lasting less than... description: The map shows the probability of the rainy season to last less than a certain number of days through all the years of available data. - #total_mean: - # menu_label: How much rain falls in the rainy season? - # description: The map shows the average precipitation (in mm) during the rainy season, - # defined as the period between onset and cessation dates, - # over all years of available data. - #total_stddev: - # menu_label: How uncertain is the amount of rain in the rainy season? - # description: The map shows the standard deviation from the average precipitation in the season, - # over all years of available data. - #total_pe: - # menu_label: How likely is it raining less than... - # description: The map shows the probability of precipitation in the rainy season - # being less than a certain amount through all the years of available data. + total_mean: + menu_label: How much rain falls in the rainy season? + description: The map shows the average precipitation (in mm) during the rainy season, + defined as the period between onset and cessation dates, + over all years of available data. + total_stddev: + menu_label: How uncertain is the amount of rain in the rainy season? + description: The map shows the standard deviation from the average precipitation in the season, + over all years of available data. + total_pe: + menu_label: How likely is it raining less than... + description: The map shows the probability of precipitation in the rainy season + being less than a certain amount through all the years of available data. # Water Balance Monitoring diff --git a/enacts/config-zmd.yaml b/enacts/config-zmd.yaml index d8b239712..9d817ebd0 100644 --- a/enacts/config-zmd.yaml +++ b/enacts/config-zmd.yaml @@ -124,4 +124,7 @@ onset: map_max: 180 length_stddev: map_max: 40 - + total_mean: + map_max: 1100 + total_stddev: + map_max: 300 diff --git a/enacts/onset/maproom.py b/enacts/onset/maproom.py index d4f4dfd8f..48162f62f 100644 --- a/enacts/onset/maproom.py +++ b/enacts/onset/maproom.py @@ -860,15 +860,19 @@ def onset_tile(tz, tx, ty): cess_soil_moisture, cess_dry_spell, ) - if "length" in map_choice: - if cess_dates["T"][0] < onset_dates["T"][0]: - cess_dates = cess_dates.isel({"T": slice(1, None)}) + if cess_dates["T"][0] < onset_dates["T"][0]: + cess_dates = cess_dates.isel({"T": slice(1, None)}) if cess_dates["T"].size != onset_dates["T"].size: - onset_dates = onset_dates.isel({"T": slice(None, -1)}) - seasonal_length = ( - (cess_dates["T"] + cess_dates["cess_delta"]).drop_indexes("T") - - (onset_dates["T"] + onset_dates["onset_delta"]).drop_indexes("T") - ) #.astype("timedelta64[D]") + onset_dates = onset_dates.isel({"T": slice(None, -2)}) + cess_dates = (cess_dates["T"] + cess_dates["cess_delta"]).drop_indexes("T") + onset_dates = (onset_dates["T"] + onset_dates["onset_delta"]).drop_indexes("T") + if "length" in map_choice: + seasonal_length = cess_dates - onset_dates + else: + seasonal_total = precip_tile.where( + (precip_tile["T"] >= onset_dates.rename({"T": "Ty"})) + & (precip_tile["T"] <= cess_dates.rename({"T": "Ty"})) + ).sum("T") if map_choice == "mean": map_data = onset_dates.onset_delta.mean("T") map_max = np.timedelta64(search_days, 'D') @@ -893,6 +897,17 @@ def onset_tile(tz, tx, ty): map_data = (seasonal_length < np.timedelta64(prob_exc_thresh_length, 'D')).mean("T") * 100 map_max = 100 colormap = pingrid.CORRELATION_COLORMAP + if map_choice == "total_mean": + map_data = seasonal_total.mean("Ty") + map_max = CONFIG["map_text"][map_choice]["map_max"] + colormap = pingrid.RAINFALL_COLORMAP + if map_choice == "total_stddev": + map_data = seasonal_total.std(dim="Ty", skipna=True) + map_max = CONFIG["map_text"][map_choice]["map_max"] + if map_choice == "total_pe": + map_data = (seasonal_total < prob_exc_thresh_tot).mean("Ty") * 100 + map_max = 100 + colormap = pingrid.CORRELATION_COLORMAP map_data.attrs["colormap"] = colormap map_data = map_data.rename(X="lon", Y="lat") map_data.attrs["scale_min"] = map_min @@ -935,6 +950,16 @@ def set_colorbar(search_start_day, search_start_month, search_days, map_choice): tick_freq = 5 map_max = CONFIG["map_text"][map_choice]["map_max"] unit = "days" + if map_choice == "total_mean": + colorbar = pingrid.to_dash_colorscale(pingrid.RAINFALL_COLORMAP) + tick_freq = 100 + map_max = CONFIG["map_text"][map_choice]["map_max"] + unit = "mm" + if map_choice == "total_stddev": + tick_freq = 30 + ) + map_max = CONFIG["map_text"][map_choice]["map_max"] + unit = "mm" if map_choice == "monit": precip = rr_mrg.precip.isel({"T": slice(-366, None)}) search_start_dm = calc.sel_day_and_month( From 5ae8a45b6279ae72af4610135636ae59efad6f94 Mon Sep 17 00:00:00 2001 From: remic Date: Thu, 23 Mar 2023 13:42:09 -0400 Subject: [PATCH 2/6] More efficient seasonal totals with masking and grouping --- enacts/onset/maproom.py | 62 +++++++++++++++++++++++++++++++++-------- 1 file changed, 51 insertions(+), 11 deletions(-) diff --git a/enacts/onset/maproom.py b/enacts/onset/maproom.py index 48162f62f..d9cbec523 100644 --- a/enacts/onset/maproom.py +++ b/enacts/onset/maproom.py @@ -21,6 +21,7 @@ from shapely import wkb from shapely.geometry.multipolygon import MultiPolygon import datetime +import xarray as xr from globals_ import FLASK, GLOBAL_CONFIG @@ -864,15 +865,55 @@ def onset_tile(tz, tx, ty): cess_dates = cess_dates.isel({"T": slice(1, None)}) if cess_dates["T"].size != onset_dates["T"].size: onset_dates = onset_dates.isel({"T": slice(None, -2)}) - cess_dates = (cess_dates["T"] + cess_dates["cess_delta"]).drop_indexes("T") - onset_dates = (onset_dates["T"] + onset_dates["onset_delta"]).drop_indexes("T") + cess_dates = cess_dates["T"] + cess_dates["cess_delta"] + onset_dates = onset_dates["T"] + onset_dates["onset_delta"] if "length" in map_choice: - seasonal_length = cess_dates - onset_dates + seasonal_length = cess_dates.drop_indexes("T") - onset_dates.drop_indexes("T") else: - seasonal_total = precip_tile.where( - (precip_tile["T"] >= onset_dates.rename({"T": "Ty"})) - & (precip_tile["T"] <= cess_dates.rename({"T": "Ty"})) - ).sum("T") + # Seasonal total between onset and cessation dates is + # a sum over yearly groups delimited by onset search dates + # of the daily rain masked for days before onset and after cessation. + # Bins can't vary with X and Y so we mask first and then we group + seasonal_total = ( + precip_tile * xr.where( + # Onset mask + xr.concat( + [ + onset_dates, + (precip_tile.isel(T=-1)*np.nan).astype("datetime64[ns]"), + ], + "T", + # Assigns Onset date to all Ts of a group + # Then masks those before onset date + ).resample(T="1D").ffill() <= precip_tile["T"], + 1, + np.nan, + ) * xr.where( + # Cessation mask + xr.concat( + [ + # We are going to use onset_dates' T to group + cess_dates.assign_coords(T=onset_dates["T"]), + (precip_tile.isel(T=-1)*np.nan).astype("datetime64[ns]"), + ], + "T", + # Assigns Cessation date to all Ts of a group + # Then masks those after cessation date + ).resample(T="1D").ffill() >= precip_tile["T"], + 1, + np.nan, + ) + # Dropping Ts that are all NA (is not necessary: helps performance?) + ).dropna(dim="T", how="all").groupby_bins( + "T", + [*(onset_dates["T"].data), np.datetime64(datetime.datetime( + onset_dates["T"][-1].dt.year.values + 1, + onset_dates["T"][-1].dt.month.values, + onset_dates["T"][-1].dt.day.values, + ))], + right=False, + include_lowest=True + ).sum().rename({"T_bins": "T"}) if map_choice == "mean": map_data = onset_dates.onset_delta.mean("T") map_max = np.timedelta64(search_days, 'D') @@ -898,14 +939,14 @@ def onset_tile(tz, tx, ty): map_max = 100 colormap = pingrid.CORRELATION_COLORMAP if map_choice == "total_mean": - map_data = seasonal_total.mean("Ty") + map_data = seasonal_total.mean("T") map_max = CONFIG["map_text"][map_choice]["map_max"] colormap = pingrid.RAINFALL_COLORMAP if map_choice == "total_stddev": - map_data = seasonal_total.std(dim="Ty", skipna=True) + map_data = seasonal_total.std(dim="T", skipna=True) map_max = CONFIG["map_text"][map_choice]["map_max"] if map_choice == "total_pe": - map_data = (seasonal_total < prob_exc_thresh_tot).mean("Ty") * 100 + map_data = (seasonal_total < prob_exc_thresh_tot).mean("T") * 100 map_max = 100 colormap = pingrid.CORRELATION_COLORMAP map_data.attrs["colormap"] = colormap @@ -957,7 +998,6 @@ def set_colorbar(search_start_day, search_start_month, search_days, map_choice): unit = "mm" if map_choice == "total_stddev": tick_freq = 30 - ) map_max = CONFIG["map_text"][map_choice]["map_max"] unit = "mm" if map_choice == "monit": From 31786327b97bea82dc20f87cfb6cc5cdafb9f4d0 Mon Sep 17 00:00:00 2001 From: remic Date: Fri, 24 Mar 2023 11:44:31 -0400 Subject: [PATCH 3/6] tot time series. Need to cover case where onset or cess is NaN --- enacts/onset/layout.py | 8 ++ enacts/onset/maproom.py | 174 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 182 insertions(+) diff --git a/enacts/onset/layout.py b/enacts/onset/layout.py index c96ac96e4..9a0763e47 100644 --- a/enacts/onset/layout.py +++ b/enacts/onset/layout.py @@ -437,6 +437,14 @@ def results_layout(): id="length_tab", label="Season Length", ), + dbc.Tab( + [ + dbc.Spinner(dcc.Graph(id="tot_plot")), + dbc.Spinner(dcc.Graph(id="tot_prob_exc")), + ], + id="tot_tab", + label="Season Total", + ), ], className="mt-4", ) diff --git a/enacts/onset/maproom.py b/enacts/onset/maproom.py index d9cbec523..7b1575dba 100644 --- a/enacts/onset/maproom.py +++ b/enacts/onset/maproom.py @@ -769,6 +769,180 @@ def length_plots( return length_graph, cdf_graph, tab_style +@APP.callback( + Output("tot_plot", "figure"), + Output("tot_prob_exc", "figure"), + Output("tot_tab","tab_style"), + Input("loc_marker", "position"), + Input("search_start_day", "value"), + Input("search_start_month", "value"), + Input("search_days", "value"), + Input("wet_threshold", "value"), + Input("running_days", "value"), + Input("running_total", "value"), + Input("min_rainy_days", "value"), + Input("dry_days", "value"), + Input("dry_spell", "value"), + Input("cess_start_day", "value"), + Input("cess_start_month", "value"), + Input("cess_search_days", "value"), + Input("cess_soil_moisture", "value"), + Input("cess_dry_spell", "value"), +) +def tot_plots( + marker_pos, + search_start_day, + search_start_month, + search_days, + wet_threshold, + running_days, + running_total, + min_rainy_days, + dry_days, + dry_spell, + cess_start_day, + cess_start_month, + cess_search_days, + cess_soil_moisture, + cess_dry_spell, +): + if not CONFIG["ison_cess_date_hist"]: + tab_style = {"display": "none"} + return {}, {}, tab_style + else: + tab_style = {} + lat = marker_pos[0] + lng = marker_pos[1] + try: + precip = pingrid.sel_snap(rr_mrg.precip, lat, lng) + isnan = np.isnan(precip).sum() + if isnan > 0: + error_fig = pingrid.error_fig(error_msg="Data missing at this location") + germ_sentence = "" + return error_fig, error_fig, tab_style + except KeyError: + error_fig = pingrid.error_fig(error_msg="Grid box out of data domain") + germ_sentence = "" + return error_fig, error_fig, tab_style + precip.load() + try: + onset_delta = calc.seasonal_onset_date( + precip, + int(search_start_day), + calc.strftimeb2int(search_start_month), + int(search_days), + int(wet_threshold), + int(running_days), + int(running_total), + int(min_rainy_days), + int(dry_days), + int(dry_spell), + time_coord="T", + ) + isnan = np.isnan(onset_delta["onset_delta"]).mean() + if isnan == 1: + error_fig = pingrid.error_fig(error_msg="No onset dates were found") + return error_fig, error_fig, tab_style + except TypeError: + error_fig = pingrid.error_fig( + error_msg="Please ensure all onset input boxes are filled." + ) + return error_fig, error_fig, tab_style + try: + soil_moisture = calc.water_balance(precip, 5, 60, 0)["soil_moisture"] + cess_delta = calc.seasonal_cess_date( + soil_moisture, + int(cess_start_day), + calc.strftimeb2int(cess_start_month), + int(cess_search_days), + int(cess_soil_moisture), + int(cess_dry_spell), + time_coord="T", + ) + isnan = np.isnan(cess_delta["cess_delta"]).mean() + if isnan == 1: + error_fig = pingrid.error_fig(error_msg="No cessation dates were found") + return error_fig, error_fig, tab_style + except TypeError: + error_fig = pingrid.error_fig( + error_msg="Please ensure all cessation input boxes are filled" + ) + return error_fig, error_fig, tab_style + if cess_delta["T"][0] < onset_delta["T"][0]: + cess_delta = cess_delta.isel({"T": slice(1, None)}) + if cess_delta["T"].size != onset_delta["T"].size: + onset_delta = onset_delta.isel({"T": slice(None, -2)}) + toto = np.stack([ + (onset_delta["T"] + onset_delta["onset_delta"]).data, + (cess_delta["T"] + cess_delta["cess_delta"]).data, + ]) + toto = toto[~np.isnan(toto).any(axis=1)] + seasonal_edges = np.stack([ + (onset_delta["T"] + onset_delta["onset_delta"]).data, + (cess_delta["T"] + cess_delta["cess_delta"]).data, + ]).flatten("F") + try: + # Total seasonal rainfall between onset and cessation dates is + # every other sum over groups defined by onset and cessation dates + # as bins edges + seasonal_total = precip.groupby_bins( + "T", + seasonal_edges, + include_lowest=True + ).sum().rename({"T_bins" :"T"})[0:-1:2] + isnan = np.isnan(seasonal_total).mean() + if isnan == 1: + error_fig = pingrid.error_fig( + error_msg="Onset or cessation not found for any season" + ) + return error_fig, error_fig, tab_style + except TypeError: + error_fig = pingrid.error_fig( + error_msg="Please ensure all onset/cessation input boxes are filled" + ) + return error_fig, error_fig, tab_style + tot_graph = pgo.Figure() + tot_graph.add_trace( + pgo.Scatter( + x=onset_delta["T"].dt.year.values, + y=seasonal_total.squeeze().values, + customdata=np.stack(( + [seasonal_total["T"].data[i].left + for i in np.arange(seasonal_total["T"].size)], + [seasonal_total["T"].data[i].right + for i in np.arange(seasonal_total["T"].size)], + ), axis=-1), + hovertemplate="%{y:d} mm from %{customdata[0]|%-d %b %Y} to %{customdata[1]|%-d %b %Y}", + name="", + line=pgo.scatter.Line(color="blue"), + ) + ) + tot_graph.update_traces(mode="lines", connectgaps=False) + tot_graph.update_layout( + xaxis_title="Year", + yaxis_title=f"Total Seasonal Rainfall in mm ", + title=f"Total Seasonal Rainfall at ({round_latLng(lat)}N,{round_latLng(lng)}E)", + ) + quantiles = np.arange(1, seasonal_total["T"].size + 1) / (seasonal_total["T"].size + 1) + tot_quant = seasonal_total.quantile(quantiles, dim="T").squeeze() + cdf_graph = pgo.Figure() + cdf_graph.add_trace( + pgo.Scatter( + x=tot_quant.values, + y=(1 - quantiles), + hovertemplate="%{y:.0%} chance of exceeding %{x:d} mm", + name="", + line=pgo.scatter.Line(color="blue"), + ) + ) + cdf_graph.update_traces(mode="lines", connectgaps=False) + cdf_graph.update_layout( + xaxis_title=f"Total Seasonal Rainfall in mm", + yaxis_title="Probability of exceeding", + ) + return tot_graph, cdf_graph, tab_style + + @FLASK.route(f"{TILE_PFX}///") def onset_tile(tz, tx, ty): parse_arg = pingrid.parse_arg From 0eff3a7faa352005056d5cec37178233302adddf Mon Sep 17 00:00:00 2001 From: remic Date: Thu, 30 Mar 2023 14:40:16 -0400 Subject: [PATCH 4/6] correct pairing of onset and cessation everywhere --- enacts/onset/maproom.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/enacts/onset/maproom.py b/enacts/onset/maproom.py index 7b1575dba..55caf4715 100644 --- a/enacts/onset/maproom.py +++ b/enacts/onset/maproom.py @@ -870,8 +870,8 @@ def tot_plots( return error_fig, error_fig, tab_style if cess_delta["T"][0] < onset_delta["T"][0]: cess_delta = cess_delta.isel({"T": slice(1, None)}) - if cess_delta["T"].size != onset_delta["T"].size: - onset_delta = onset_delta.isel({"T": slice(None, -2)}) + if cess_delta["T"].size != onset_delta["T"].size: + onset_delta = onset_delta.isel({"T": slice(None, -1)}) toto = np.stack([ (onset_delta["T"] + onset_delta["onset_delta"]).data, (cess_delta["T"] + cess_delta["cess_delta"]).data, @@ -1037,8 +1037,8 @@ def onset_tile(tz, tx, ty): ) if cess_dates["T"][0] < onset_dates["T"][0]: cess_dates = cess_dates.isel({"T": slice(1, None)}) - if cess_dates["T"].size != onset_dates["T"].size: - onset_dates = onset_dates.isel({"T": slice(None, -2)}) + if cess_dates["T"].size != onset_dates["T"].size: + onset_dates = onset_dates.isel({"T": slice(None, -1)}) cess_dates = cess_dates["T"] + cess_dates["cess_delta"] onset_dates = onset_dates["T"] + onset_dates["onset_delta"] if "length" in map_choice: From 28c406d9d7cd896e054e58d9478b1fe02b5a76b6 Mon Sep 17 00:00:00 2001 From: remic Date: Thu, 30 Mar 2023 14:41:10 -0400 Subject: [PATCH 5/6] bar plot for total --- enacts/onset/maproom.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/enacts/onset/maproom.py b/enacts/onset/maproom.py index 55caf4715..4763a7e05 100644 --- a/enacts/onset/maproom.py +++ b/enacts/onset/maproom.py @@ -903,7 +903,7 @@ def tot_plots( return error_fig, error_fig, tab_style tot_graph = pgo.Figure() tot_graph.add_trace( - pgo.Scatter( + pgo.Bar( x=onset_delta["T"].dt.year.values, y=seasonal_total.squeeze().values, customdata=np.stack(( @@ -914,10 +914,8 @@ def tot_plots( ), axis=-1), hovertemplate="%{y:d} mm from %{customdata[0]|%-d %b %Y} to %{customdata[1]|%-d %b %Y}", name="", - line=pgo.scatter.Line(color="blue"), ) ) - tot_graph.update_traces(mode="lines", connectgaps=False) tot_graph.update_layout( xaxis_title="Year", yaxis_title=f"Total Seasonal Rainfall in mm ", From 2e599496c22ebb6bf9d6ad4cc4e66b55d1469522 Mon Sep 17 00:00:00 2001 From: remic Date: Thu, 30 Mar 2023 16:04:46 -0400 Subject: [PATCH 6/6] dealing with missing onset/cess, correct indexing, better use of for --- enacts/onset/maproom.py | 26 +++++++++++--------------- 1 file changed, 11 insertions(+), 15 deletions(-) diff --git a/enacts/onset/maproom.py b/enacts/onset/maproom.py index 4763a7e05..7e8a952b6 100644 --- a/enacts/onset/maproom.py +++ b/enacts/onset/maproom.py @@ -872,15 +872,13 @@ def tot_plots( cess_delta = cess_delta.isel({"T": slice(1, None)}) if cess_delta["T"].size != onset_delta["T"].size: onset_delta = onset_delta.isel({"T": slice(None, -1)}) - toto = np.stack([ - (onset_delta["T"] + onset_delta["onset_delta"]).data, - (cess_delta["T"] + cess_delta["cess_delta"]).data, - ]) - toto = toto[~np.isnan(toto).any(axis=1)] - seasonal_edges = np.stack([ - (onset_delta["T"] + onset_delta["onset_delta"]).data, - (cess_delta["T"] + cess_delta["cess_delta"]).data, - ]).flatten("F") + seasonal_edges = xr.concat( + [ + (onset_delta["T"] + onset_delta["onset_delta"]).drop_vars("T"), + (cess_delta["T"] + cess_delta["cess_delta"]).drop_vars("T"), + ], + dim="edges" + ).where(lambda x: ~np.isnan(x).any(axis=0), drop=True).data.flatten("F") try: # Total seasonal rainfall between onset and cessation dates is # every other sum over groups defined by onset and cessation dates @@ -889,7 +887,7 @@ def tot_plots( "T", seasonal_edges, include_lowest=True - ).sum().rename({"T_bins" :"T"})[0:-1:2] + ).sum().rename({"T_bins" :"T"})[0::2] isnan = np.isnan(seasonal_total).mean() if isnan == 1: error_fig = pingrid.error_fig( @@ -904,13 +902,11 @@ def tot_plots( tot_graph = pgo.Figure() tot_graph.add_trace( pgo.Bar( - x=onset_delta["T"].dt.year.values, + x=[interval.left.year for interval in seasonal_total["T"].data], y=seasonal_total.squeeze().values, customdata=np.stack(( - [seasonal_total["T"].data[i].left - for i in np.arange(seasonal_total["T"].size)], - [seasonal_total["T"].data[i].right - for i in np.arange(seasonal_total["T"].size)], + [interval.left for interval in seasonal_total["T"].data], + [interval.right for interval in seasonal_total["T"].data], ), axis=-1), hovertemplate="%{y:d} mm from %{customdata[0]|%-d %b %Y} to %{customdata[1]|%-d %b %Y}", name="",