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)

Solve Manning's Equation with Python Scipy library

from scipy.optimize import root

def manningC(d, args):
    Q, w,h,sSlopeL,sSlopeR,nMann,lSlope = args
    #left side slope can be different from right side slope
    area = ((((d*sSlopeL)+(d*sSlopeR)+w)+w)/2)*d
    # wet perimeter
    wPer = w+(d*(sSlopeL*sSlopeL+1)**0.5)+(d*(sSlopeR*sSlopeR+1)**0.5)
    #Hydraulic Radius
    hR = area/ wPer
    # following formula must be zero
    # manipulation of Manning's formula
    mannR = (Q*nMann/lSlope**0.5)-(area*hR**(2.0/3.0))
    return mannR
    
###### MAIN CODE
# the following are input data to our open channel manning calculation
# flow, width, height, left side slope, right side slope, 
# Manning coefficient, longitudinal slope
args0 = [2.5,2,.5,1.0,1.0,.015,.005]
initD = .00001  # initial water depth value

# then we call the root scipy function to the manningC
sol =root(manningC,initD, args=(args0,))    
# print the root found
print(sol.x)

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