rogerlew / wepppy

Other
12 stars 7 forks source link

return period rankings #317

Open rogerlew opened 1 month ago

rogerlew commented 1 month ago

script to test jim's return period rankings

import pandas as pd
import numpy as np

def weibull_series(recurrence, years):
    """
    this came from Jim F.'s code. recurrence is a list of recurrence intervals. years is the number
    of years in the simulation. For each RI it determines the rank event index to estimate the return period.

    Not sure where Jim got it.
    """

    print(f"weibull_series")

    recurrence = sorted(recurrence)

    rec = {}
    i = 0
    rankind = years
    orgind = years + 1
    reccount = 0

    while i < len(recurrence) and rankind >= 2.5:
        retperiod = recurrence[i]
        rankind = float(years + 1) / retperiod
        intind = int(rankind) - 1

        if intind < orgind:
            rec[retperiod] = intind
            orgind = intind
            reccount += 1

        i += 1

    return rec

def weibull_series2(recurrence, years, gringorten_correction=False):
    print(f"weibull_series2 (Gringorten Correction: {gringorten_correction})")
    ndays = int(round(years * 365.25))
    ranks = np.array(list(range(1, ndays + 1)))
    if gringorten_correction:
        Ts = (ndays + 1) / (ranks - 0.44)
    else:
        Ts = (ndays + 1) / ranks

    yearTs = Ts / 365.25

    rec = {}
    for rec_interval in sorted(recurrence):
        for _rank, _yearTs in zip(ranks[::-1], yearTs[::-1]):
            _rank_indx = _rank - 1
            if _yearTs >= rec_interval and _rank_indx not in rec.values():
                rec[rec_interval] = _rank_indx
                break

    return rec

def tabulizer(rec, years):
    df = pd.DataFrame(list(rec.items()), columns=['Return Period', 'Rank Index'])
    df['Rank'] = df['Rank Index'] + 1
    df['T'] = (years + 1) / (df['Rank'])
    df['Tg'] = (years + 1) / (df['Rank']) - 0.44
    df['Years'] = years
    print(df.to_string(index=False))
    print()

if __name__ == "__main__":
    tabulizer(weibull_series([2, 5, 10, 20, 50, 100], 30), 30)
    tabulizer(weibull_series2([2, 5, 10, 20, 50, 100], 30), 30)

    tabulizer(weibull_series([2, 5, 10, 20, 50, 100], 48), 48)
    tabulizer(weibull_series2([2, 5, 10, 20, 50, 100], 48), 48)

    tabulizer(weibull_series([2, 5, 10, 20, 50, 100], 50), 50)
    tabulizer(weibull_series2([2, 5, 10, 20, 50, 100], 50), 50)

    tabulizer(weibull_series([2, 5, 10, 20, 50, 100], 100), 100)
    tabulizer(weibull_series2([2, 5, 10, 20, 50, 100], 100), 100)

    tabulizer(weibull_series([2, 5, 10, 20, 50, 100], 123), 123)
    tabulizer(weibull_series2([2, 5, 10, 20, 50, 100], 123), 123)

    tabulizer(weibull_series([2, 5, 10, 20, 50, 100], 124), 124)
    tabulizer(weibull_series2([2, 5, 10, 20, 50, 100], 124), 124)

    tabulizer(weibull_series([2, 5, 10, 20, 50, 100, 500], 500), 500)
    tabulizer(weibull_series2([2, 5, 10, 20, 50, 100, 500], 500), 500)
rogerlew commented 1 month ago
weibull_series
 Return Period  Rank Index  Rank         T        Tg  Years
             2          14    15  2.066667  1.626667     30
             5           5     6  5.166667  4.726667     30
            10           2     3 10.333333  9.893333     30
            20           0     1 31.000000 30.560000     30

weibull_series2 (Gringorten Correction: False)
 Return Period  Rank Index  Rank         T        Tg  Years
             2          14    15  2.066667  1.626667     30
             5           5     6  5.166667  4.726667     30
            10           2     3 10.333333  9.893333     30
            20           0     1 31.000000 30.560000     30

weibull_series
 Return Period  Rank Index  Rank         T        Tg  Years
             2          23    24  2.041667  1.601667     48
             5           8     9  5.444444  5.004444     48
            10           3     4 12.250000 11.810000     48
            20           1     2 24.500000 24.060000     48

weibull_series2 (Gringorten Correction: False)
 Return Period  Rank Index  Rank         T        Tg  Years
             2          23    24  2.041667  1.601667     48
             5           8     9  5.444444  5.004444     48
            10           3     4 12.250000 11.810000     48
            20           1     2 24.500000 24.060000     48

weibull_series
 Return Period  Rank Index  Rank     T    Tg  Years
             2          24    25  2.04  1.60     50
             5           9    10  5.10  4.66     50
            10           4     5 10.20  9.76     50
            20           1     2 25.50 25.06     50
            50           0     1 51.00 50.56     50

weibull_series2 (Gringorten Correction: False)
 Return Period  Rank Index  Rank     T    Tg  Years
             2          24    25  2.04  1.60     50
             5           9    10  5.10  4.66     50
            10           4     5 10.20  9.76     50
            20           1     2 25.50 25.06     50
            50           0     1 51.00 50.56     50

weibull_series
 Return Period  Rank Index  Rank     T    Tg  Years
             2          49    50  2.02  1.58    100
             5          19    20  5.05  4.61    100
            10           9    10 10.10  9.66    100
            20           4     5 20.20 19.76    100
            50           1     2 50.50 50.06    100

weibull_series2 (Gringorten Correction: False)
 Return Period  Rank Index  Rank      T     Tg  Years
             2          49    50   2.02   1.58    100
             5          19    20   5.05   4.61    100
            10           9    10  10.10   9.66    100
            20           4     5  20.20  19.76    100
            50           1     2  50.50  50.06    100
           100           0     1 101.00 100.56    100

weibull_series
 Return Period  Rank Index  Rank         T        Tg  Years
             2          61    62  2.000000  1.560000    123
             5          23    24  5.166667  4.726667    123
            10          11    12 10.333333  9.893333    123
            20           5     6 20.666667 20.226667    123
            50           1     2 62.000000 61.560000    123

weibull_series2 (Gringorten Correction: False)
 Return Period  Rank Index  Rank          T         Tg  Years
             2          60    61   2.032787   1.592787    123
             5          23    24   5.166667   4.726667    123
            10          11    12  10.333333   9.893333    123
            20           5     6  20.666667  20.226667    123
            50           1     2  62.000000  61.560000    123
           100           0     1 124.000000 123.560000    123

weibull_series
 Return Period  Rank Index  Rank          T         Tg  Years
             2          61    62   2.016129   1.576129    124
             5          24    25   5.000000   4.560000    124
            10          11    12  10.416667   9.976667    124
            20           5     6  20.833333  20.393333    124
            50           1     2  62.500000  62.060000    124
           100           0     1 125.000000 124.560000    124

weibull_series2 (Gringorten Correction: False)
 Return Period  Rank Index  Rank          T         Tg  Years
             2          61    62   2.016129   1.576129    124
             5          23    24   5.208333   4.768333    124
            10          11    12  10.416667   9.976667    124
            20           5     6  20.833333  20.393333    124
            50           1     2  62.500000  62.060000    124
           100           0     1 125.000000 124.560000    124

weibull_series
 Return Period  Rank Index  Rank       T      Tg  Years
             2         249   250   2.004   1.564    500
             5          99   100   5.010   4.570    500
            10          49    50  10.020   9.580    500
            20          24    25  20.040  19.600    500
            50           9    10  50.100  49.660    500
           100           4     5 100.200  99.760    500
           500           0     1 501.000 500.560    500

weibull_series2 (Gringorten Correction: False)
 Return Period  Rank Index  Rank       T      Tg  Years
             2         249   250   2.004   1.564    500
             5          99   100   5.010   4.570    500
            10          49    50  10.020   9.580    500
            20          24    25  20.040  19.600    500
            50           9    10  50.100  49.660    500
           100           4     5 100.200  99.760    500
           500           0     1 501.000 500.560    500

The rankings make sense. Jim's algorithm does have a short coming where it doesn't calculate 100 year return periods until there is 124 years of data

rogerlew commented 1 month ago

method revised for complete time or annual maxima

def weibull_series2(recurrence, years, method='cta', gringorten_correction=False):
    """
    Generates a Weibull distribution for return periods based on a given number of years and recurrence intervals.

    Args:
        recurrence (list): A list of recurrence intervals (in years) for which to find corresponding ranks.
        years (float): The total number of years for the time series.
        method (str): The method used to calculate the return periods. Options are 'cta' (default) complete time series analysis or 'am' or annual maxima.
        gringorten_correction (bool): If True, applies the Gringorten correction to the Weibull formula.

    Returns:
        dict: A dictionary where keys are the recurrence intervals and values are the ranks in the time series 
              that correspond to the given recurrence intervals.

    Explanation:
        The function calculates the recurrence times (return periods) using the Weibull formula, which can be corrected 
        using the Gringorten correction if specified. It then identifies and returns the ranks that correspond to 
        the provided recurrence intervals.
    """
    print(f"weibull_series2 (method: {method} gringorten_correction: {gringorten_correction})")

    # Calculate the total number of days based on the number of years
    if method == 'cta':
        n = int(round(years * 365.25))  # 365.25 accounts for leap years
    elif method == 'am':
        n = int(years)

    # Create an array of ranks from 1 to the total number of days
    ranks = np.array(list(range(1, n + 1)))

    # Apply the Weibull formula with or without Gringorten correction
    if gringorten_correction:
        Ts = (n + 1) / (ranks - 0.44)  # Gringorten correction applied
    else:
        Ts = (n + 1) / ranks  # Standard Weibull formula without correction

    # Convert return periods (Ts) to years by dividing by the number of days in a year
    if method == 'cta':
        yearTs = Ts / 365.25
    elif method == 'am':
        yearTs = Ts

    # Initialize a dictionary to store recurrence intervals and their corresponding rank indices
    rec = {}

    # Loop through the sorted recurrence intervals to find corresponding rank indices
    for rec_interval in sorted(recurrence):
        # Reverse the ranks and yearTs arrays to start from the highest values
        for _rank, _yearTs in zip(ranks[::-1], yearTs[::-1]):
            _rank_indx = _rank - 1  # Adjust rank to 0-based index
            # Check if the current return period is greater than or equal to the recurrence interval
            # and if the rank index is not already assigned
            if _yearTs >= rec_interval and _rank_indx not in rec.values():
                rec[rec_interval] = _rank_indx  # Assign the rank index to the recurrence interval
                break  # Move on to the next recurrence interval

    return rec
rogerlew commented 1 month ago
weibull_series
 Return Period  Rank Index  Rank         T        Tg  Years
             2          14    15  2.066667  1.626667     30
             5           5     6  5.166667  4.726667     30
            10           2     3 10.333333  9.893333     30
            20           0     1 31.000000 30.560000     30

weibull_series2 (method: cta gringorten_correction: False)
 Return Period  Rank Index  Rank         T        Tg  Years
             2          14    15  2.066667  1.626667     30
             5           5     6  5.166667  4.726667     30
            10           2     3 10.333333  9.893333     30
            20           0     1 31.000000 30.560000     30

weibull_series2 (method: am gringorten_correction: False)
 Return Period  Rank Index  Rank         T        Tg  Years
             2          14    15  2.066667  1.626667     30
             5           5     6  5.166667  4.726667     30
            10           2     3 10.333333  9.893333     30
            20           0     1 31.000000 30.560000     30

weibull_series
 Return Period  Rank Index  Rank         T        Tg  Years
             2          23    24  2.041667  1.601667     48
             5           8     9  5.444444  5.004444     48
            10           3     4 12.250000 11.810000     48
            20           1     2 24.500000 24.060000     48

weibull_series2 (method: cta gringorten_correction: False)
 Return Period  Rank Index  Rank         T        Tg  Years
             2          23    24  2.041667  1.601667     48
             5           8     9  5.444444  5.004444     48
            10           3     4 12.250000 11.810000     48
            20           1     2 24.500000 24.060000     48

weibull_series2 (method: am gringorten_correction: False)
 Return Period  Rank Index  Rank         T        Tg  Years
             2          14    15  3.266667  2.826667     48
             5           5     6  8.166667  7.726667     48
            10           2     3 16.333333 15.893333     48
            20           0     1 49.000000 48.560000     48

weibull_series
 Return Period  Rank Index  Rank     T    Tg  Years
             2          24    25  2.04  1.60     50
             5           9    10  5.10  4.66     50
            10           4     5 10.20  9.76     50
            20           1     2 25.50 25.06     50
            50           0     1 51.00 50.56     50

weibull_series2 (method: cta gringorten_correction: False)
 Return Period  Rank Index  Rank     T    Tg  Years
             2          24    25  2.04  1.60     50
             5           9    10  5.10  4.66     50
            10           4     5 10.20  9.76     50
            20           1     2 25.50 25.06     50
            50           0     1 51.00 50.56     50

weibull_series2 (method: am gringorten_correction: False)
 Return Period  Rank Index  Rank    T    Tg  Years
             2          14    15  3.4  2.96     50
             5           5     6  8.5  8.06     50
            10           2     3 17.0 16.56     50
            20           0     1 51.0 50.56     50

weibull_series
 Return Period  Rank Index  Rank     T    Tg  Years
             2          49    50  2.02  1.58    100
             5          19    20  5.05  4.61    100
            10           9    10 10.10  9.66    100
            20           4     5 20.20 19.76    100
            50           1     2 50.50 50.06    100

weibull_series2 (method: cta gringorten_correction: False)
 Return Period  Rank Index  Rank      T     Tg  Years
             2          49    50   2.02   1.58    100
             5          19    20   5.05   4.61    100
            10           9    10  10.10   9.66    100
            20           4     5  20.20  19.76    100
            50           1     2  50.50  50.06    100
           100           0     1 101.00 100.56    100

weibull_series2 (method: am gringorten_correction: False)
 Return Period  Rank Index  Rank      T     Tg  Years
             2          49    50   2.02   1.58    100
             5          19    20   5.05   4.61    100
            10           9    10  10.10   9.66    100
            20           4     5  20.20  19.76    100
            50           1     2  50.50  50.06    100
           100           0     1 101.00 100.56    100

weibull_series
 Return Period  Rank Index  Rank         T        Tg  Years
             2          61    62  2.000000  1.560000    123
             5          23    24  5.166667  4.726667    123
            10          11    12 10.333333  9.893333    123
            20           5     6 20.666667 20.226667    123
            50           1     2 62.000000 61.560000    123

weibull_series2 (method: cta gringorten_correction: False)
 Return Period  Rank Index  Rank          T         Tg  Years
             2          60    61   2.032787   1.592787    123
             5          23    24   5.166667   4.726667    123
            10          11    12  10.333333   9.893333    123
            20           5     6  20.666667  20.226667    123
            50           1     2  62.000000  61.560000    123
           100           0     1 124.000000 123.560000    123

weibull_series2 (method: am gringorten_correction: False)
 Return Period  Rank Index  Rank          T         Tg  Years
             2          61    62   2.000000   1.560000    123
             5          23    24   5.166667   4.726667    123
            10          11    12  10.333333   9.893333    123
            20           5     6  20.666667  20.226667    123
            50           1     2  62.000000  61.560000    123
           100           0     1 124.000000 123.560000    123

weibull_series
 Return Period  Rank Index  Rank          T         Tg  Years
             2          61    62   2.016129   1.576129    124
             5          24    25   5.000000   4.560000    124
            10          11    12  10.416667   9.976667    124
            20           5     6  20.833333  20.393333    124
            50           1     2  62.500000  62.060000    124
           100           0     1 125.000000 124.560000    124

weibull_series2 (method: cta gringorten_correction: False)
 Return Period  Rank Index  Rank          T         Tg  Years
             2          61    62   2.016129   1.576129    124
             5          23    24   5.208333   4.768333    124
            10          11    12  10.416667   9.976667    124
            20           5     6  20.833333  20.393333    124
            50           1     2  62.500000  62.060000    124
           100           0     1 125.000000 124.560000    124

weibull_series2 (method: am gringorten_correction: False)
 Return Period  Rank Index  Rank          T         Tg  Years
             2          61    62   2.016129   1.576129    124
             5          24    25   5.000000   4.560000    124
            10          11    12  10.416667   9.976667    124
            20           5     6  20.833333  20.393333    124
            50           1     2  62.500000  62.060000    124
           100           0     1 125.000000 124.560000    124

weibull_series
 Return Period  Rank Index  Rank       T      Tg  Years
             2         249   250   2.004   1.564    500
             5          99   100   5.010   4.570    500
            10          49    50  10.020   9.580    500
            20          24    25  20.040  19.600    500
            50           9    10  50.100  49.660    500
           100           4     5 100.200  99.760    500
           500           0     1 501.000 500.560    500

weibull_series2 (method: cta gringorten_correction: False)
 Return Period  Rank Index  Rank       T      Tg  Years
             2         249   250   2.004   1.564    500
             5          99   100   5.010   4.570    500
            10          49    50  10.020   9.580    500
            20          24    25  20.040  19.600    500
            50           9    10  50.100  49.660    500
           100           4     5 100.200  99.760    500
           500           0     1 501.000 500.560    500

weibull_series2 (method: am gringorten_correction: False)
 Return Period  Rank Index  Rank       T      Tg  Years
             2         249   250   2.004   1.564    500
             5          99   100   5.010   4.570    500
            10          49    50  10.020   9.580    500
            20          24    25  20.040  19.600    500
            50           9    10  50.100  49.660    500
           100           4     5 100.200  99.760    500
           500           0     1 501.000 500.560    500
rogerlew commented 1 month ago

Complete Time Series Analysis (CTA): Return periods are estimated by ranking all daily events in the time series from largest to smallest. Weibull T-statistics are then calculated for each rank. The first event that exceeds the desired recurrence interval is selected as the event corresponding to that return period.

Annual Maxima (AM) Analysis: Return periods are estimated by first identifying the largest event from each year. These annual maxima are ranked, and Weibull T-statistics are calculated for each ranked event. The event that exceeds the recurrence interval is selected in the same way as in the CTA method.