Estimar parámetros de la curva de calibración de caudales o descarga con R

#################################################################################
#
#
# Estimación de los parámetros de l curva de descarga
# Ajusta y plotea la grafica de curva descarga siguiendo una funcion exponencial
#
#  Walter Arnoldo Bardales Espinoza
#
#################################################################################

# remover todos los objetos de almacenamiento
rm(list = ls())

# cargar libreria
require(Hmisc)
library(minpack.lm)

# leer datos de aforos y niveles
dato.aforos = read.table("D:/Documents/Curva_calibracion_de_caudales/est_Matucuy.txt", header = TRUE, sep = "\t",na.strings = "-99.9")


plot(Caudal~Nivel)

# Para usar los nombres de los campos de las columnas
attach(dato.aforos)


# Estima los valores de los parametros iniciales de la curva
# Estableciendo h0 a un valor por debajo del minimo observado
rango=range(Nivel)
h.min=rango[1]
h.max=rango[2]
ho.inicial =  h.min - 0.1*(h.max - h.min)

# modelo lineal para estimar los valores inicial de k y n
modelo.log = lm( log10(Caudal) ~ log10(Nivel - ho.inicial) )
a.inicial = 10^modelo.log$coefficients[[1]]
b.inicial = modelo.log$coefficients[[2]]

# Utilizar los coeficientes calculados inicialmente para estimar los parametros de la regresion no lineal
curva.descarga = nlsLM( Caudal ~ k*(Nivel-ho)^n, start = list(k=a.inicial,n=b.inicial, ho=ho.inicial), control = list(maxiter=500))

# Resumen de los estadisticos de la curva de descarga
summary(curva.descarga)

# Grafico de la curva de descarga o calibracion
xmin = min(Nivel)
xmax = max(Nivel)
dx = 0.01*(xmax - xmin)
x = seq(xmin,xmax,dx)
coef.model = coef(curva.descarga)
qest = coef.model[1]*(x - coef.model[3])^coef.model[2]


plot(Nivel, Caudal,
     ylim = c(0, max(Caudal)+0.2*(max(Caudal)-min(Caudal))),
     xlim = c(min(Nivel)-0.1*(rango[2]-rango[1]), max(Nivel)+0.2*(rango[2]-rango[1])),
     xlab = "Nivel (m)",
     ylab = expression("Q ("*m^3*s^{-1}*")"),
     pch = 21, bg = "black",
     main = "Curva de calibracion de caudales de la estacion Matucuy, Periodo 2002 - 2015"
)
lines(x, qest, col = "red")
text(min(Nivel)+0.1*(max(Nivel)-min(Nivel)),max(Caudal), paste0("Q = ",round(coef.model[1],3),"*(H-(",round(coef.model[3],3),"))^",round(coef.model[2],3)))
text(min(Nivel)+0.1*(max(Nivel)-min(Nivel)),0.95*max(Caudal),paste0("RMS =", round(mean((Caudal-(coef.model[1]*(Nivel - coef.model[3])^coef.model[2]))^2))))

#################################################################### Fin código

Descargar datos y scripts R