Skip to content

Commit

Permalink
Revert "baseflow recession added"
Browse files Browse the repository at this point in the history
  • Loading branch information
@cheginit
cheginit authored May 26, 2024
1 parent 0773f35 commit a9f251e
Showing 1 changed file with 0 additions and 222 deletions.
222 changes: 0 additions & 222 deletions hydrosignatures/hydrosignatures.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down

0 comments on commit a9f251e

Please sign in to comment.