From 1e8a2b57144b743e4115684f55f44c91785ea4ec Mon Sep 17 00:00:00 2001 From: Dan Howe Date: Sun, 15 May 2022 17:13:00 +0200 Subject: [PATCH] Allow SLR to be generated with an external function --- .../probabilistic_assessment.py | 88 +++++++++++-------- 1 file changed, 52 insertions(+), 36 deletions(-) diff --git a/probabilistic-analysis/probabilistic_assessment.py b/probabilistic-analysis/probabilistic_assessment.py index c52fb91..98a97d1 100644 --- a/probabilistic-analysis/probabilistic_assessment.py +++ b/probabilistic-analysis/probabilistic_assessment.py @@ -169,6 +169,7 @@ def get_ongoing_recession(n_runs, start_year, end_year, sea_level_rise, 'min' (array_like): minimum value 'mode' (array_like): most likely value 'max' (array_like): maximum value + 'function' (str): optional external function ('package.function') bruun_factor (dict): 'min' (float): minimum value 'mode' (float): most likely value @@ -180,6 +181,12 @@ def get_ongoing_recession(n_runs, start_year, end_year, sea_level_rise, Returns: the simulated recession distance (m) + + Notes: + 'sea_level_rise' is calculated with a triangular probability distribution + by default. Alternatively 'sea_level_rise' can be calculated using an + external function, to which the arguments 'n_runs', 'start_year', and + 'end_year' are passed. """ # Get time interval from input file @@ -187,45 +194,54 @@ def get_ongoing_recession(n_runs, start_year, end_year, sea_level_rise, n_years = len(years) # Interpolate sea level rise projections (m) - slr_mode = np.interp(years, - xp=sea_level_rise['year'], - fp=sea_level_rise['mode'])[:, np.newaxis] - - try: - slr_min = np.interp(years, - xp=sea_level_rise['year'], - fp=sea_level_rise['min'])[:, np.newaxis] - except ValueError: - # Use mode for deterministic beaches - slr_min = slr_mode - - try: - slr_max = np.interp(years, - xp=sea_level_rise['year'], - fp=sea_level_rise['max'])[:, np.newaxis] - except ValueError: - # Use mode for deterministic beaches - slr_max = slr_mode - - # Initialise sea level rise array - slr = np.zeros([n_years, n_runs]) - - for i in range(n_years): - # Use triangular distribution for SLR in each year (m) + if sea_level_rise['function']: # Get slr from separate function + # Get names of package/script and function + pkg, func_name = sea_level_rise['function'].split('.') + + # Import function from package + func = getattr(__import__(pkg), func_name) + slr = func(n_runs=n_runs, start_year=start_year, end_year=end_year) + + else: # Use triangular distribution + slr_mode = np.interp(years, + xp=sea_level_rise['year'], + fp=sea_level_rise['mode'])[:, np.newaxis] + try: - slr[i, :] = np.random.triangular(left=slr_min[i], - mode=slr_mode[i], - right=slr_max[i], - size=n_runs) + slr_min = np.interp(years, + xp=sea_level_rise['year'], + fp=sea_level_rise['min'])[:, np.newaxis] except ValueError: - # Use constant value if slr_min == slr_max - slr[i, :] = np.ones([1, n_runs]) * slr_mode[i] + # Use mode for deterministic beaches + slr_min = slr_mode - # Sort each row, so SLR follows a smooth trajectory for each model run - slr = np.sort(slr, axis=1) - - # Shuffle columns, so the order of model runs is randomised - slr = np.random.permutation(slr.T).T + try: + slr_max = np.interp(years, + xp=sea_level_rise['year'], + fp=sea_level_rise['max'])[:, np.newaxis] + except ValueError: + # Use mode for deterministic beaches + slr_max = slr_mode + + # Initialise sea level rise array + slr = np.zeros([n_years, n_runs]) + + for i in range(n_years): + # Use triangular distribution for SLR in each year (m) + try: + slr[i, :] = np.random.triangular(left=slr_min[i], + mode=slr_mode[i], + right=slr_max[i], + size=n_runs) + except ValueError: + # Use constant value if slr_min == slr_max + slr[i, :] = np.ones([1, n_runs]) * slr_mode[i] + + # Sort each row, so SLR follows a smooth trajectory for each model run + slr = np.sort(slr, axis=1) + + # Shuffle columns, so the order of model runs is randomised + slr = np.random.permutation(slr.T).T # Shift sea level so it is zero in the start year slr -= slr[0, :].mean(axis=0)