Calculating the R-factor according to the ABAG

rfactor.py

# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import pylab as plt
import datetime as dt
from glob import glob
import sys

def calculate_r_factor(t0, ri_mm_h):
    R_i = []
    R_cum = []
    R_per_year = {}
    tseries = []
    Ecum = I30max = Rcum = 0.0
    year = 2012
    zero_counter = 0
    for i, I30 in enumerate(ri_mm_h):
        if I30 > 0:
            zero_counter = 0
            I30max = max(I30max, I30)
            Ni = I30 * 0.5 # mm/30min
            if I30 < 0.05:
                Ei = 0
            elif I30 < 76.2:
                Ei = (11.89 + 8.73 * np.log10(I30)) * Ni * 1e-3
            else:
                Ei = 28.33 * Ni * 1e-3
            Ecum += Ei
        else: # no rain
            zero_counter += 1
        if zero_counter == 12:
            # Rain is over
            t = t0 + dt.timedelta(hours=(i - 12) / 2)
            Ri = Ecum * I30max
            if year != t.year:
                Rcum = 0
                year = t.year
            Rcum += Ri
            R_per_year[t.year] = R_per_year.get(t.year, 0.0) + Ri
            R_i.append(Ri)
            R_cum.append(Rcum)
            tseries.append(plt.date2num(t))
            Ecum = I30max = 0.0
    y, Ry = zip(*R_per_year.items())
    R_yearly = pd.Series(data=Ry, index=[pd.Period(_y) for _y in y])
        
    R_i = pd.Series(data=R_i, index=plt.num2date(tseries))
    R_cum = pd.Series(data=R_cum, index=plt.num2date(tseries))

    return R_i, R_cum, R_yearly

def load_dwd(*filenames):
    def conv(s):
        return dt.datetime.strptime(s, '%Y%m%d%H%M')
    
    dataframes = []
    for fn in sorted(filenames):
        print(fn)
        df = pd.read_csv(fn, converters={1: conv}, sep=';',skipinitialspace=True)
        df.index = df.MESS_DATUM
        df = df['RWS_10']
        df[df == -999] = 0.0
        df *= 6 # Convert from mm/10min to mm/h
        if dataframes:
            df = df[df.index > dataframes[-1].index.max()]
        dataframes.append(df)
    return pd.concat(dataframes).resample('30min').mean()
        
def load_schwingbach(filename):
    data = np.genfromtxt(filename, delimiter=",", 
                         filling_values=0.0, 
                         skip_header=1, usecols=(1,)
                         )
    data[np.isnan(data)] = 0.0

    ri_mm_h = data / 24
    return ri_mm_h

        
def plot_ri(series):
    
    plt.figure(figsize=(6, 4.5))
    series.plot()
    plt.ylabel('Rain intensity mm/h')
    plt.xlabel('')
    plt.grid()
    plt.savefig('rainintensity.png', dpi=300)
    
def plot_rf(R_i, R_cum, tseries):
    
    plt.figure(figsize=(6, 4.5))
    R_i.plot('line', style='g-')
    R_cum.plot('line', style='r-')
    plt.ylabel('R - factor N/h')
    plt.xlabel('')
    plt.grid()
    plt.savefig('rfactor.png', dpi=300)
 
if __name__ == '__main__':
    filenames = sys.argv[1:]
    if 'series' not in globals():
       series = load_dwd(*filenames)
    R_i, R_cum, R_per_year = calculate_r_factor(series.index[0].to_pydatetime(), series)