Fitting IDF curves with Scipy and Pandas

# -*- coding: utf-8 -*-
"""
Created on Tue Jul 10 11:19:37 2018

@author: rodrigo.goncalves
"""
import pandas as pd
from scipy.optimize import minimize
 
# This is the IDF function, returning the sum of squared errors (SSE)
def func2(par, res):
    p1 =  (par[0] * res.index.values  **par[1])
    p2 = ((res.columns.values.astype(float)+par[2])**par[3])
    f = pd.DataFrame([p/p2 for p in p1],index=res.index,columns=res.columns)
    erroTotQ=((f-res)**2).sum(axis=1).sum()    
    return erroTotQ
 
# copy your rainfall intensities table from excel
# with column headers with rainfall durations
# and row names with Return period value in years
dfInt=pd.read_clipboard()

#initial guess
param1 = [5000, 0.1, 10, 0.9]

res2 = minimize(func2, param1, args=(dfInt,), method='Nelder-Mead')
print(res2)
cs=['K=','a=','b=','c=']
dfResult=pd.DataFrame(res2.x,index=cs).transpose()
print(dfResult)
dfResult.to_clipboard(index=None)

SciPy minimize example - Fitting IDF Curves

SciPy (pronounced “Sigh Pie”) is an open source Python library used by scientists, analysts, and engineers doing scientific computing and technical computing.

SciPy contains modules for optimization, linear algebra, integration, interpolation, special functions, FFT, signal and image processing, ODE solvers and other tasks common in science and engineering.

In this post I will show how to use a powerful function of SciPy - minimize.

Minimize has some methods of minimizing functions. Its official documentation is shown here.

The example used is to find coefficients of a standard rainfall IDF (intensity-duration-frequency) equation. The format of the equation is as following:

We have to input the intensity data to be fitted, and the equation as a function.

import numpy as np
from scipy.optimize import minimize

# This is the IDF function, for return periods of 10 and 25 years
def func2(par, res):
    f10 =  (par[0] * 10  **par[1])/((res[0,:]+par[2])**par[3])
    f25 =  (par[0] * 25  **par[1])/((res[0,:]+par[2])**par[3])
    erroTotQ = np.sum((f10-res[1,:])**2+(f25-res[2,:])**2)
    return erroTotQ

#durations array - in minutes
d=[5,10,15,20,25,30,35,40,60,120,180,360,540,720,900,1260,1440]

# Rainfall intensities, mm/h, same lenght as minutes
r10=[236.1,174.0,145.6,128.3,116.3,107.3,100.3,94.5,79.1,48.4,34.2,18.9,13.4,10.5,8.7,6.5,5.8]
r25=[294.5,217.1,181.6,160.0,145.1,133.9,125.1,118.0,98.7,60.4,42.7,23.6,16.7,13.1,10.8,8.1,7.2]

# the following line only gather the duration and the intensities as a numpy array,
# in order to pass all of the idf constants as one single parameter, named "valid"
valid = np.vstack((d, r10, r25))

#initial guess
param1 = [5000, 0.1, 10, 0.9]

res2 = minimize(func2, param1, args=(valid,), method='Nelder-Mead')
np.set_printoptions(formatter={'float': '{: 0.3f}'.format})
print res2