From a9f251e70ca2790c52aba481cbf092e57bdb65f2 Mon Sep 17 00:00:00 2001 From: Taher Chegini Date: Sun, 26 May 2024 09:00:21 -0400 Subject: [PATCH] Revert "baseflow recession added" --- hydrosignatures/hydrosignatures.py | 222 ----------------------------- 1 file changed, 222 deletions(-) diff --git a/hydrosignatures/hydrosignatures.py b/hydrosignatures/hydrosignatures.py index 019d1e9..e960a15 100644 --- a/hydrosignatures/hydrosignatures.py +++ b/hydrosignatures/hydrosignatures.py @@ -493,228 +493,6 @@ def extract_extrema(ts: pd.Series, var_name: str, n_pts: int) -> pd.DataFrame: ts_df.loc[ts_df.iloc[ilocs_min].index, "trough"] = True return ts_df -class BaseflowRecession: - def __init__(self, streamflow, ndmin=7, ndmax=30, nregmx=300): - self.streamflow = streamflow - self.ndmin = ndmin - self.ndmax = ndmax - self.ndays = len(streamflow) - self.nregmx = nregmx - - self.alpha = np.zeros(self.ndays) - self.ndreg = np.zeros(self.ndays, dtype=int) - self.q0 = np.zeros(self.ndays) - self.q10 = np.zeros(self.ndays) - self.bfdd = np.zeros(self.nregmx) - self.qaveln = np.zeros(self.nregmx) - self.florec = np.zeros((self.nregmx, 200)) - self.aveflo = np.zeros(self.nregmx) - self.npreg = np.zeros(self.nregmx, dtype=int) - self.idone = np.zeros(self.nregmx, dtype=int) - self.icount = np.zeros(200, dtype=int) - - def process_streamflow_data(self): - self.calculate_alpha() - self.process_alpha_data() - self.estimate_master_recession_curve(self.npr) - - def calculate_alpha(self): - """Calculate the recession constant alpha and baseflow days.""" - nd = 0 - self.npr = 0 - - for i in range(1, self.ndays): - if np.isnan(self.streamflow[i]) or np.isnan(self.streamflow[i - 1]): - continue - if self.streamflow[i] <= 0: - nd = 0 - else: - if self.streamflow[i] / self.streamflow[i - 1] < 0.99: - if nd >= self.ndmin: - if not np.isnan(self.streamflow[i - nd]) and not np.isnan(self.streamflow[i - 1]): - self.alpha[i] = np.log(self.streamflow[i - nd] / self.streamflow[i - 1]) / nd - self.ndreg[i] = nd - nd = 0 - else: - nd += 1 - if nd >= self.ndmax: - if not np.isnan(self.streamflow[i - nd + 1]) and not np.isnan(self.streamflow[i]): - self.alpha[i] = np.log(self.streamflow[i - nd + 1] / self.streamflow[i]) / nd - self.ndreg[i] = nd - nd = 0 - - def process_alpha_data(self): - """Process the alpha values to compute the recession constant.""" - for i in range(1, self.ndays): - if self.alpha[i] > 0: - if self.npr < self.nregmx - 1: - self.npr += 1 - if not np.isnan(self.streamflow[i - 1]) and not np.isnan(self.streamflow[i - self.ndreg[i]]): - self.q10[self.npr] = self.streamflow[i - 1] - self.q0[self.npr] = self.streamflow[i - self.ndreg[i]] - if self.q0[self.npr] - self.q10[self.npr] > 0.001: - self.bfdd[self.npr] = self.ndreg[i] / (np.log(self.q0[self.npr]) - np.log(self.q10[self.npr])) - self.qaveln[self.npr] = np.log((self.q0[self.npr] + self.q10[self.npr]) / 2) - self.fill_florec(self.npr, i) - else: - self.npr -= 1 - - def fill_florec(self, npr, i): - """Fill the florec array with log-transformed flow data.""" - kk = 0 - for k in range(self.ndreg[i]): - if not np.isnan(self.streamflow[i - k]): - x = np.log(self.streamflow[i - k]) - if x > 0: - kk += 1 - if kk < self.florec.shape[0] and npr < self.florec.shape[1]: - self.florec[kk, npr] = x - if kk == 0 and npr > 0: - npr -= 1 - - def estimate_master_recession_curve(self, npr): - """Estimate the master recession curve.""" - if npr > 1: - n_points, sumx, sumy, sumxy, sumx2 = 0, 0, 0, 0, 0 - for i in range(1, npr + 1): - if not np.isnan(self.qaveln[i]) and not np.isnan(self.bfdd[i]): - n_points += 1 - x = self.qaveln[i] - sumx += x - sumy += self.bfdd[i] - sumxy += x * self.bfdd[i] - sumx2 += x * x - - if n_points > 1: - ssxy = sumxy - (sumx * sumy) / n_points - ssxx = sumx2 - (sumx * sumx) / n_points - slope_ = ssxy / ssxx - yint = sumy / n_points - slope_ * sumx / n_points - - self.fill_aveflo(slope_, yint, npr) - self.calculate_final_alf(n_points, sumx, sumy, sumxy, sumx2, npr) - else: - self.alf = np.nan - self.bfd = np.nan - else: - self.alf = np.nan - self.bfd = np.nan - - def linear_regression_fill(self, slope_, yint, now): - """Perform linear regression to fill aveflo array.""" - n_points, sumx, sumy, sumxy, sumx2 = 0, 0, 0, 0, 0 - for i in range(1, self.nregmx + 1): - if i < self.aveflo.size and not np.isnan(self.aveflo[i]) and self.aveflo[i] > 0: - n_points += 1 - x = float(i) - sumx += x - sumy += self.aveflo[i] - sumxy += x * self.aveflo[i] - sumx2 += x * x - - if n_points > 1: - ssxy = sumxy - (sumx * sumy) / n_points - ssxx = sumx2 - (sumx * sumx) / n_points - if ssxx != 0: - slope_ = ssxy / ssxx - else: - slope_ = np.nan - - if np.isfinite(slope_): - yint = sumy / n_points - slope_ * sumx / n_points - if now < self.florec.shape[1]: - if slope_ != 0: - self.icount[now] = int((self.florec[1, now] - yint) / slope) - else: - self.icount[now] = 0 - else: - self.icount[now] = 0 - else: - self.icount[now] = 0 - - def fill_aveflo(self, slope_, yint, npr): - """Fill the aveflo array and adjust counts.""" - for j in range(1, npr + 1): - amn = np.inf - for i in range(1, npr + 1): - if self.idone[i] == 0 and i < self.florec.shape[1]: - if self.florec[1, i] < amn: - amn = self.florec[1, i] - now = i - if now >= self.florec.shape[1]: - continue - self.idone[now] = 1 - - igap = 0 - if j == 1: - self.icount[now] = 1 - igap = 1 - else: - for i in range(1, self.nregmx + 1): - if i < self.aveflo.size and self.florec[1, now] <= self.aveflo[i]: - self.icount[now] = i - igap = 1 - break - - if igap == 0: - self.linear_regression_fill(slope_, yint, now) - - for i in range(1, self.ndmax + 1): - k = self.icount[now] + i - 1 - if k < 0 or k >= self.aveflo.size: - continue - if self.florec[i, now] > 0.0001: - self.aveflo[k] = (self.aveflo[k] * self.npreg[k] + self.florec[i, now]) / (self.npreg[k] + 1) - if self.aveflo[k] <= 0: - self.aveflo[k] = slope_ * i + yint - self.npreg[k] += 1 - - def calculate_final_alf(self, n_points, sumx, sumy, sumxy, sumx2, npr): - """Calculate the final alpha factor (alf) and baseflow days (bfd).""" - n_points, sumx, sumy, sumxy, sumx2 = 0, 0, 0, 0, 0 - for j in range(1, min(npr + 1, self.florec.shape[1])): - for i in range(1, min(self.ndmax + 1, self.florec.shape[0])): - if not np.isnan(self.florec[i, j]) and self.florec[i, j] > 0: - n_points += 1 - x = float(self.icount[j] + i) - sumx += x - sumy += self.florec[i, j] - sumxy += x * self.florec[i, j] - sumx2 += x * x - - if n_points > 1: - ssxy = sumxy - (sumx * sumy) / n_points - ssxx = sumx2 - (sumx * sumx) / n_points - if ssxx != 0: - self.alf = ssxy / ssxx - self.bfd = 2.3 / self.alf - else: - self.alf = np.nan - self.bfd = np.nan - else: - self.alf = np.nan - self.bfd = np.nan - - def get_results(self): - """Return the results as a dictionary.""" - # Extract the data points for the master recession curve - days = [] - flow_rates = [] - for i in range(self.florec.shape[1]): - if np.any(self.florec[:, i] > 0): - valid_indices = np.where(self.florec[:, i] > 0)[0] - days.extend(valid_indices + 1) - flow_rates.extend(self.florec[valid_indices, i]) - - return { - 'alpha': self.alpha, - 'baseflow_days': self.bfd, - 'alf': self.alf, - 'npr': self.npr, - 'master_rc_days': days, - 'master_rc_flow_rates': flow_rates - } - @dataclass(frozen=True) class SignaturesBool: