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

Relleno de datos de lluvia mensual de una estación con CHIRPS

Anteriormente, generamos una función para descargar los raster de lluvia estimada a partir de satélite (CHIRPS V 2.0), ahora lo que vamos hacer es extraer los datos para un punto en especifico o la ubicación de una estación climática, el periodo de tiempo de la serie de datos chirps debe ser la misma que la serie de datos de la estación (no importa que existan datos faltantes, solo especificar el datos faltante como -99.9) para que sean comparables.

Luego se creara un modelo de ajuste entre los datos estimados y los datos reales, esto con el fin de que los datos a rellenar sean acordes a la tendencia de los datos observados o pluviometricos.

La importancia de completar series de datos radica darle continuidad a la serie y que esto me permita realizar mejores análisis, crear modelos hidrológicos estables, balances hídricos, etc.

####################################################################

library(raster)

# Ruta del directorio de trabajo
setwd("D:/Documents/chirps/")

# Directorio de los datos de chirp V-2.0
ruta_chirps = "D:/Documents/chirps/raster/"

# Periodo de datos a analizar
fecha_inicial = c("1981/1/1")
fecha_final = c("2018/5/1")

# Coordenada geografica de las estacion (longitud, Latitud)
coor.x = c(-89.8106)
coor.y = c(15.6083)

# Convertir las coordenadas de las estaciones a puntos
ubicacion = data.frame(lon=coor.x, lat=coor.y)
coordinates(ubicacion) = c("lon", "lat")
proj4string(ubicacion) = CRS("+init=epsg:4326") # WGS 84

# Cargar los datos raster
lista_chirps = paste0(ruta_chirps,"chirps-v2.0.", format(seq(as.Date(fecha_inicial), as.Date(fecha_final), "month"), "%Y.%m"), ".tif")
ras_chirps = stack(lista_chirps)

# Extraer datos de lluvia de la estacion
chirps = as.vector((extract(ras_chirps, ubicacion)))

# Cargar datos pluviometricos (reales u observados) de la estacion
estacion = scan("estaciones/Cahabon_1_1981_5_2018.txt", sep = "\t", na.strings = "-99.9")
serie = estacion


# Modelo de ajuste de datos chirps a pluviometicos
model = lm(estacion~chirps)

# estadisticos del modelo de ajuste lineal
summary(model)

# coeficientes del modelo de ajuste lineal
alfa = coef(model)[[2]]
beta = coef(model)[[1]]

# Grafica
plot(estacion~chirps)
abline(b=1,a=1, col="red")

# posiciones de datos faltantes pluviometricos
faltantes = which(is.na(serie))

# Datos chirps ajustados para relleno
relleno = round(alfa*chirps[faltantes]+beta,1)

# Completacion de la serie real
serie[faltantes] = relleno

# Grafico de serie completa (negro) y serie incompleta (rojo)
plot(serie,type="l", xlab="Lluvia [mm]", main="Comparacion de serie completa e incompleta de lluvia")
lines(estacion, col="red")

# Crear tabla de fechas y serie completa
fecha = format(seq(as.Date(fecha_inicial), as.Date(fecha_final), "month"), "%Y.%m")
serie.completa = data.frame(fecha,serie)

write.csv(serie.completa, "Cahabon_rellenada.csv", row.names = FALSE)

#############################################################

Descargar
Archivo_R

Delimitacion de cuencas con qgis y grass

https://youtu.be/ytsEYNy4Xn4

Extraer datos de lluvia en una cuenca usando R

La falta de datos de lluvia en una cuenca, no es impedimento para poder realizar un balance hídrico hoy en día, ya que existen datos de hidroestimación o datos de lluvia estimada a partir de satelites o radares.  Estos datos tienen la ventaja que proporcionan la distribución espacial de la lluvia, pero su desventaja es que tienen a subestimar o sobrestimar en algunos eventos de lluvia y esto depende del tipo de lluvia; pero mediante un ajuste por regresión se puede reducir el error.
 También, pueden usarse para rellenar datos mensuales, trimestrales, anuales de lluvia en las series observadas, mediante un ajuste.

## Script para descargar datos chirps 2.0 y extraer la lluvia media de la cuenca
## Walter Bardales
## Fecha: 01/07/2018


#################################################################################

library(R.utils)

# Funcion de descarga Chirps V2.0
download.chirps2.0 = function(Ruta, fecha_inicial, fecha_final){
fecha=format(seq(as.Date(fecha_inicial), as.Date(fecha_final), "month"), format="%Y.%m")
for (i in 1:length(fecha)){
wurl=paste0("ftp://ftp.chg.ucsb.edu/pub/org/chg/products/CHIRPS-2.0/camer-carib_monthly/tifs/","chirps-v2.0.",fecha[i],".tif.gz")
download.file(wurl,paste0(Ruta,"chirps-v2.0.",fecha[i],".tif.gz"))
gunzip(paste0(Ruta,"chirps-v2.0.",fecha[i],".tif.gz"))
}
}

#################################################################################

# Cargar libreria de trabajo
library(raster)
library(rgdal)

# Directorio de trabajo
setwd("D:/Documents/chirps/")

# Cargar funcion de descarga chirps V2.0
source("Funcion_descarga_chirps.R")

# Directorio de la descarga
Ruta = "D:/Documents/chirps/raster/"

# Periodo de descarga
fecha_inicial = c("1983/5/1")
fecha_final = c("1983/8/1")

# Descargar archivos chirps y descomprimir
download.chirps2.0(Ruta,fecha_inicial, fecha_final)

# Desactivar el paquete R.utils por tener conflicto con el paquete Raster
detach("package:R.utils", unload=TRUE)

# Cargar los datos de chirps 2.0
lista_archivos = paste0("raster/","chirps-v2.0.",format(seq(as.Date(fecha_inicial), as.Date(fecha_final), "month"),"%Y.%m"),".tif")
chirps = stack(lista_archivos)

# Cargar el poligono de la cuenca en coordenadas geograficas
cuenca = readOGR("shapefile/La_Pasion.shp") # Indicar el directorio de la cuenca

# Extraer los datos
ppcuenca = t(extract(chirps, cuenca, fun=mean))
row.names(ppcuenca) = NULL # Crear serie de datos y guardar
lluvia_media = data.frame(fecha=format(seq(as.Date(fecha_inicial), as.Date(fecha_final), "month"),"%Y.%m"), Lluvia=ppcuenca)

write.csv(lluvia_media, "Lluvia_La_Pasion.csv", row.names = FALSE) # Cambiar el nombre del archivo de salida

#######################################################################

Descargar scripts

Archivo

Graficar imagenes goes con R

library(ncdf4)
library(raster)
library(rgdal)
library(fields)
library(RColorBrewer)

#Directorio de trabajo
setwd("D:/Documents/Goes_R/netcdf")

# Shapefile
shp1 = readOGR("D:/Documents/GIS DataBase/Shape_WGS84/departamentos.shp")
shp2 = readOGR("D:/Documents/GIS DataBase/Shape_WGS84/Centro_America.shp")

# Proyeccion de goes
geost = CRS("+proj=geos +a=6378137. +b=6356752.31414 +lon_0=-75 +f=0.00335281068119356027 +h=35786023. +sweep=x +no_defs")

# Coordenadas geograficas
geo = CRS("+proj=longlat +datum=WGS84")

# Leer archivo netcdf y voltearlo verticalmente
goes_16 = flip(raster("OR_ABI-L2-CMIPF-M3C14_G16_s20181742345453_e20181742356220_c20181742356317.nc"),direction = "y")
crs(goes_16) = geost # Asignar la proyeccion del satelite geoestacionario
extent(goes_16) = c(-5434894.885056, 5434894.885056, -5434894.885056, 5434894.885056) # Asignarle la extension

# Recortar el area de interes (Centro america)
cgoes_16 = crop(goes_16,extent(c(-3424153, -231218.2, 791968.4, 3216755)))

# Reproyectar de GRS80 a coordenadas WGS84
cgoes_16 = projectRaster(cgoes_16, crs=geo)

# Plotear y agregar el mapa de centro america y Departamentos Guatemala
plot(cgoes_16 - 273.15, xlab="Longitud", ylab="Latitud", main="Banda 14")

plot(cgoes_16-273.15, xlab="Longitud", ylab="Latitud", main="Banda 14", xlim=c(-105, -67), ylim=c(8,30), col=rainbow(20, alpha = .6))
lines(shp2)
lines(shp1)

Prueba de infiltración

Prueba de infiltración y ajuste con el modelo de Horton o Kostiakov

Infiltración (I):

La teoría de infiltración está basada en un análisis del movimiento del agua en el suelo bajo condiciones de no saturación. Durante la infiltración, la fase líquida y la fase gaseosa coexisten en la masa de suelo, con excepción de la zona de contacto entre el suelo y el agua en la superficie del terreno.

La infiltración, velocidad de infiltración o intensidad de entrada en el suelo, se puede definir como la velocidad de penetración del agua en el perfil del suelo, cuando la superficie del terreno está cubierta por una capa de agua poco profunda.

La infiltración tiene dimensión de velocidad (L T^-1), como la lámina de agua (L) admitida por el suelo en una unidad de tiempo (T) o como la cantidad de agua absorbida por la unidad de superficie del terreno en la unidad de tiempo (L³ T^-1 L^-2), respectivamente. Si las mismas unidades se usan en ambos casos, las expresiones son dimensionalmente equivalentes (L T^-1). En la primera forma la expresión común de velocidad de entrada es cm h^-1 o cm min^-1.  En la segunda forma, generalmente se expresa como m³ min^-1 m^-2. Es común, por varios autores, señalar que la infiltración obedece a un fenómeno de movimiento vertical del agua en el perfil de suelo.  A la disminución de la velocidad de infiltración después de cierto período de tiempo alcanzando una velocidad que tiende a ser constante y se denomina infiltración básica (I_b), la cual se considera cercana a la conductividad hidráulica o permeabilidad (K).

Factores que afectan la infiltración

·       Textura del suelo: Proporción de arenas, limos y arcilla

·       Estructura del suelo: Es la forma en la que se encuentra posicionadas las partículas del suelo o pedón (granular, migajosa, bloques, prismática o columnar o laminar)

·       Contenido de humedad del suelo (Capacidad de campo, punto de marchitez permanente y saturación)

Métodos para medir la velocidad de infiltración:

Las técnicas más utilizadas para determinar tasas de infiltración y conductividad hidráulica saturada, son: Lagunas de infiltración, permeámetros, infiltrómetros cilíndricos y el método de Porchet.

Porchet

Para la estimación de la tasa de infiltración en terreno se puede utilizar el método de Porchet, el cual consiste en excavar un cilindro de radio (R) y se llenarlo con agua hasta una altura (h).

La superficie por la cual pasa el agua al infiltrarse es:

S=2π*R*h+ π*R^2

S=πR(2h+R)

La velocidad de infiltración (f) que pasa por la superficie (S) del agujero cilíndrico es proporcional a la variación del volumen de agua del agujero cilíndrico para un intervalo de tiempo infinitamente pequeño (dt):

S*f =  dV/dt

Como el radio permanece constante, la expresión se puede simplificar de la siguiente forma:

f=  (πR^2)/S dh/dt

Sustituir la superficie del agujero cilíndrico

f=  (πR^2)/(πR(2h+R))  dh/dt

Al simplificar

f=  R/(2h+R)  dh/dt

Se parar los diferenciales

f dt=  R/((2h+R))  dh

Al integrar y valuar en los intervalos de tiempo (t1, t2) y el nivel de agua (h1 y h2), siendo t12

y h21

f (t2-t1)=  R/2*(Ln(2*h1+R)- Ln(2*h2+R))

Al simplificar

f=  R/(2(t2-t1))*Ln((2*h1+R)/(2*h2+R))

Métodos de ajuste de la ecuación de infiltración

Horton:

          Horton estudio la capacidad de infiltración, en su estudio determinó que la capacidad de infiltración tiene la siguiente forma.

fp= fc + (fo-fc )*e^(-kt)  

Donde:

fp= Capacidad de infiltración (mm/hr)

fc= Capacidad final o equilibrio (mm/hr)
fo= Capacidad inicial (mm/hr)

K= Coeficiente de decrecimiento de la capacidad de infiltración en el tiempo

Kostiakov

          Este investigador propuso una ecuación empírica, basada en un modelo potencial

f=a*t^b

Donde:

f= velocidad de infiltración (mm/min)

a y b = parámetros de ajuste

t= Tiempo (min)

Introduccion a la estadistica climatica 1