import * as stats from "./stats.js";
/**
* Main class used for hydrological analyses.
* @class
* @name hydro
*/
export default class hydro {
/**
* Create stage-discharge rating curves using various regression methods
* Fundamental tool for converting water level measurements to flow rates
* Supports power law, polynomial, linear, and datum-corrected regressions
*
* @function stageDischarge
* @memberof hydro
* @param {Object} options - Function options
* @param {Object} [options.params] - Parameters
* @param {string} [options.params.type] - Regression type: 'power' (default), 'power_corrected', 'linear', 'polynomial'
* @param {number} [options.params.order=2] - Polynomial order (2 or 3, only for polynomial type)
* @param {number} [options.params.steps=50] - Number of curve points to generate
* @param {boolean} [options.params.onlyCurve=false] - Return only curve array (exclude coefficients/stats)
* @param {Object} [options.args] - Additional arguments (not used)
* @param {Array} options.data - Input data: [[stage, discharge], ...] or [{stage, discharge}, ...]
* @returns {Object|Array} Regression results with coefficients, equation, R², and curve; or just curve if onlyCurve=true
*
* @example
* // Power law rating curve (most common): Q = A * D^B
* const measurements = [
* [1.2, 10], // [stage (m), discharge (m³/s)]
* [1.5, 25],
* [1.8, 45],
* [2.1, 70],
* [2.4, 100]
* ];
* const rating = hydro.analyze.hydro.stageDischarge({
* params: { type: 'power' },
* data: measurements
* });
* console.log(`Equation: Q = ${rating.coefficients.A.toFixed(2)} * D^${rating.coefficients.B.toFixed(2)}`);
* console.log(`R²: ${rating.r2.toFixed(4)}`);
*
* @example
* // Power law with datum correction: Q = A * (D - h0)^B
* // Use when datum adjustment improves fit
* const rating = hydro.analyze.hydro.stageDischarge({
* params: { type: 'power_corrected' },
* data: measurements
* });
* console.log(`Datum offset (h0): ${rating.coefficients.h0.toFixed(3)} m`);
*
* @example
* // Polynomial regression for complex relationships
* const rating = hydro.analyze.hydro.stageDischarge({
* params: { type: 'polynomial', order: 3 },
* data: measurements
* });
*
* @example
* // Generate discharge values for specific stages
* const rating = hydro.analyze.hydro.stageDischarge({
* params: { type: 'power', steps: 100 },
* data: measurements
* });
* // Use rating.curve to plot or rating.predict(stage) for individual values
* console.log(`Flow at 2.0m: ${rating.predict(2.0).toFixed(2)} m³/s`);
*
* @example
* // Object format input
* const data = [
* { stage: 1.2, discharge: 10 },
* { stage: 1.5, discharge: 25 },
* { stage: 1.8, discharge: 45 }
* ];
* const rating = hydro.analyze.hydro.stageDischarge({
* params: { type: 'power' },
* data
* });
*/
static stageDischarge({ params = {}, args = {}, data = [] } = {}) {
const type = params.type || 'power';
const steps = params.steps || 50;
const onlyCurve = params.onlyCurve || false;
const polyOrder = params.order || 2;
// 1. Preprocess Data
let pairs = [];
if (Array.isArray(data)) {
if (data.length > 0 && Array.isArray(data[0])) {
// [[stage, discharge], ...]
pairs = data.map(d => ({ stage: Number(d[0]), discharge: Number(d[1]) }));
} else if (data.length > 0 && typeof data[0] === 'object') {
// [{stage: ..., discharge: ...}, ...]
pairs = data.map(d => ({ stage: Number(d.stage || d.x), discharge: Number(d.discharge || d.y || d.val) }));
}
}
// Filter invalid
pairs = pairs.filter(p => !isNaN(p.stage) && !isNaN(p.discharge));
if (pairs.length < 2) return { error: "Insufficient data" };
// Sort by stage for curve generation
pairs.sort((a, b) => a.stage - b.stage);
const stages = pairs.map(p => p.stage);
const discharges = pairs.map(p => p.discharge);
// Helpers
const sum = arr => arr.reduce((a, b) => a + b, 0);
const mean = arr => sum(arr) / arr.length;
// Simple Linear Regression Solver (Least Squares)
// Returns { m, c, r2, predict(x) }
const solveLinear = (x, y) => {
const n = x.length;
const sx = sum(x);
const sy = sum(y);
const sxy = sum(x.map((xi, i) => xi * y[i]));
const sxx = sum(x.map(xi => xi * xi));
const syy = sum(y.map(yi => yi * yi));
const denominator = (n * sxx - sx * sx);
if (denominator === 0) return null;
const m = (n * sxy - sx * sy) / denominator;
const c = (sy - m * sx) / n;
// R2 calculation
const yMean = sy / n;
const ssRes = sum(y.map((yi, i) => Math.pow(yi - (m * x[i] + c), 2)));
const ssTot = sum(y.map(yi => Math.pow(yi - yMean, 2)));
const r2 = 1 - (ssRes / ssTot);
return { m, c, r2, predict: (val) => m * val + c };
};
// Polynomial Solver (Gaussian Elimination)
// Returns { coeffs: [a0, a1, ...], r2, predict(x) }
// Equation: y = a0 + a1*x + a2*x^2 ...
const solvePoly = (x, y, order) => {
const n = x.length;
const m = order + 1;
const X = [];
const Y = [];
// Build Normal Equations
for (let i = 0; i < m; i++) {
X[i] = [];
for (let j = 0; j < m; j++) {
X[i][j] = sum(x.map(xi => Math.pow(xi, i + j)));
}
Y[i] = sum(x.map((xi, k) => Math.pow(xi, i) * y[k]));
}
// Solve X * a = Y using Gaussian elimination
const a = new Array(m).fill(0);
for (let i = 0; i < m; i++) {
let pivot = i;
for (let j = i + 1; j < m; j++) {
if (Math.abs(X[j][i]) > Math.abs(X[pivot][i])) pivot = j;
}
[X[i], X[pivot]] = [X[pivot], X[i]];
[Y[i], Y[pivot]] = [Y[pivot], Y[i]];
for (let j = i + 1; j < m; j++) {
const factor = X[j][i] / X[i][i];
for (let k = i; k < m; k++) X[j][k] -= factor * X[i][k];
Y[j] -= factor * Y[i];
}
}
for (let i = m - 1; i >= 0; i--) {
let sumAX = 0;
for (let j = i + 1; j < m; j++) sumAX += X[i][j] * a[j];
a[i] = (Y[i] - sumAX) / X[i][i];
}
// Predict & R2
const predict = (val) => a.reduce((acc, coef, i) => acc + coef * Math.pow(val, i), 0);
const yMean = mean(y);
const ssRes = sum(y.map((yi, i) => Math.pow(yi - predict(x[i]), 2)));
const ssTot = sum(y.map(yi => Math.pow(yi - yMean, 2)));
const r2 = 1 - (ssRes / ssTot);
return { coeffs: a, r2, predict };
};
let resultModel = {};
// --- Regression Logic ---
if (type === 'linear') {
const res = solveLinear(stages, discharges);
if (!res) return { error: "Linear regression failed" };
resultModel = {
coefficients: { m: res.m, c: res.c },
r2: res.r2,
equation: `Q = ${res.m.toFixed(4)} * D + ${res.c.toFixed(4)}`,
predict: res.predict
};
} else if (type === 'polynomial') {
const order = Math.min(Math.max(polyOrder, 2), 4); // Limit order 2-4
const res = solvePoly(stages, discharges, order);
resultModel = {
coefficients: { a: res.coeffs },
r2: res.r2,
equation: `Q = ${res.coeffs.map((c, i) => `${c.toFixed(4)}*D^${i}`).join(' + ')}`,
predict: res.predict
};
} else if (type === 'power_corrected') {
// Q = A * (D - h0)^B
// optimize h0 to maximize R2 of ln(Q) vs ln(D - h0)
const minStage = Math.min(...stages);
// Search range for h0: slightly less than minStage down to reasonable lower bound determine by spread
// Heuristic: try h0 from minStage - (range) to minStage - epsilon
// Or simply simple grid search for now.
let bestR2 = -Infinity;
let bestParams = null;
// Step size for h0 search.
// We assume real zero flow stage is somewhat below min observed stage.
const range = Math.max(...stages) - minStage;
// We search h0 in [minStage - 2*range, minStage - 0.001*range]
// This is a rough heuristic.
const potentialH0s = [];
const stepsH0 = 50;
for (let i = 1; i <= stepsH0; i++) {
potentialH0s.push(minStage - (i / stepsH0) * range * 0.5); // Search just below
// Also try deeper if needed, but usually h0 is close to min stage or slightly below.
}
// Also try very small offsets
for (let i = 1; i <= 10; i++) potentialH0s.push(minStage - i * 0.01);
for (const h0 of potentialH0s) {
// Transform
const validPairs = pairs.filter(p => (p.stage - h0) > 0 && p.discharge > 0);
if (validPairs.length < 3) continue;
const lnD = validPairs.map(p => Math.log(p.stage - h0));
const lnQ = validPairs.map(p => Math.log(p.discharge));
const lin = solveLinear(lnD, lnQ);
if (lin && lin.r2 > bestR2) {
bestR2 = lin.r2;
bestParams = {
h0,
B: lin.m,
A: Math.exp(lin.c)
};
}
}
if (!bestParams) return { error: "Power corrected regression fitted poorly" };
const predict = (D) => {
if (D <= bestParams.h0) return 0;
return bestParams.A * Math.pow(D - bestParams.h0, bestParams.B);
};
resultModel = {
coefficients: bestParams,
r2: bestR2,
equation: `Q = ${bestParams.A.toFixed(4)} * (D - ${bestParams.h0.toFixed(4)})^${bestParams.B.toFixed(4)}`,
predict: predict
};
} else {
// Default: 'power' (Simple Power Law)
// Q = A * D^B -> ln(Q) = ln(A) + B * ln(D)
const validPairs = pairs.filter(p => p.stage > 0 && p.discharge > 0);
if (validPairs.length < 2) return { error: "Data must be positive for power law" };
const lnD = validPairs.map(p => Math.log(p.stage));
const lnQ = validPairs.map(p => Math.log(p.discharge));
const res = solveLinear(lnD, lnQ);
if (!res) return { error: "Power regression failed" };
const A = Math.exp(res.c);
const B = res.m;
resultModel = {
coefficients: { A, B },
r2: res.r2,
equation: `Q = ${A.toFixed(4)} * D^${B.toFixed(4)}`,
predict: (val) => A * Math.pow(val, B)
};
}
// Generate Curve
const minS = Math.min(...stages);
const maxS = Math.max(...stages);
const curveStep = (maxS - minS) / (steps - 1);
const curveStages = [];
const curveDischarges = [];
for (let i = 0; i < steps; i++) {
const s = minS + i * curveStep;
curveStages.push(s);
curveDischarges.push(resultModel.predict(s));
}
if (onlyCurve) {
return [curveStages, curveDischarges];
}
return {
...resultModel,
curve: [curveStages, curveDischarges]
};
}
// constructor(){
// //Defining stats variable to be used throughout the component
// this.stats = new stats()
// }
/**
* Computation of aereal mean precipitation for a river basin given it has 2 or more different stations.
* @method arithmetic
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Not used by this function.
* @param {Array} data - Contains: array object with precipitation with equal amounts of data from different rain gauges.
* @returns {Array} Array with object with average precipitation for a specific time series.
* @example
* hydro.analyze.hydro.arithmetic({data: [[1, 2, 3, 4], [2, 3, 5, 2]]});
*/
static arithmetic({ params, args, data } = {}) {
var arr = data,
average = [],
final = [],
n = arr.length;
for (var i = 0; i < arr.length; i++) {
for (var j = 0; j < arr[0].length; j++) {
average.push(arr[i][j]);
}
}
for (var h = 0; h < average.length; h = +n) {
var minarr = average.slice(h, h + n);
final[h] = this.totalprec(minarr) / minarr.length;
var filtered = final.filter(function (el) {
return el != null;
});
}
return filtered;
}
/**
* Computes the Standardized Precipitation Index (SPI) using gamma fit and inverse normal transformation.
* @method spiCompute
* @memberof stats
* @param {Object} params - Parameters for SPI computation.
* @param {number} params.scale - Aggregation scale for SPI (e.g., 12 for SPI-12).
* @param {Object} args - Not used by this function.
* @param {Object} data - Input precipitation time series (1D array).
* @returns {Array} SPI time series.
* @example
* hydro.analyze.hydro.spiCompute({ params: { scale: 12 }, data: [10, 20, 30, 40, 50] });
*/
static spiCompute({ params = {}, args = {}, data = [] } = {}) {
const scale = Number(params.scale || 12);
const EPSILON = 1e-6;
const n = data.length;
if (n < scale) return [];
const m = n - scale + 1;
const result = new Array(m);
for (let i = 0; i < m; i++) {
let sum = 0;
for (let j = 0; j < scale; j++) {
let val = data[i + j];
if (val < 0.0) val = 0.0;
sum += val;
}
result[i] = sum;
}
const valid = result.filter(x => x > EPSILON);
if (valid.length < 10) return result.map(() => 0);
const mean = stats.default.mean({ data: valid });
const log_mean = stats.default.mean({ data: valid.map(x => Math.log(x)) });
const s = Math.log(mean) - log_mean;
const alpha = (1 + Math.sqrt(1 + (4 * s) / 3)) / (4 * s);
const beta = mean / alpha;
function invNormCDF(p) {
const t = Math.sqrt(-2.0 * Math.log(1.0 - p));
let z = t - (2.515517 + 0.802853 * t + 0.010328 * t * t) /
(1.0 + 1.432788 * t + 0.189269 * t * t + 0.001308 * t * t * t);
return p < 0.5 ? -z : z;
}
return result.map(x => {
if (x > EPSILON) {
const px = stats.default.gammaCDFApprox({ params: { alpha, beta }, data: [x] });
const clampedPx = Math.max(Math.min(px, 1 - EPSILON), EPSILON);
return invNormCDF(clampedPx);
}
return 0;
});
}
/**
* Thiessen polygon method for rainfall areal averaging.
* Calculates the weighted average of point rainfall observations by dividing the study area into polygonal subareas and weighting each observation by the proportion of the total area it contributes to.
* @method thiessen
* @memberof hydro
* @param {Object} params - Contains: areas (array with the areas of each polygon)
* @param {Object} args - Not used by this function.
* @param {Array} data - Contains: an array of arrays with rainfall observations for each point
* @returns {Array} Array with the weighted average of point rainfall observations for each time step
* @example
* hydro.analyze.hydro.thiessen({params: {areas: [10, 20, 30]}, data: [[0.3, 0.2, 0.5], [0.2, 0.3, 0.4], [0.1, 0.4, 0.5]]});
*/
static thiessen({ params, args, data } = {}) {
const areas = params.areas.map(a => Number(a));
// Check if number of data arrays matches number of areas, fill remaining with empty arrays
const numAreas = areas.length;
const numDataArrays = data.length;
if (numDataArrays < numAreas) {
const emptyArrays = Array(numAreas - numDataArrays).fill([]);
data = [...data, ...emptyArrays];
}
// Check if areas array contains valid numbers
if (areas.some(isNaN)) {
throw new Error("Invalid input: areas array contains non-numeric values");
}
const totalarea = this.totalprec({ data: areas });
// Find the longest rainfall array and set its length as the number of columns
const maxLength = Math.max(...data.map((arr) => arr.length));
const rainfallData = data.map((arr) =>
Array.from({ length: maxLength }, (_, i) => arr[i] ?? 0)
);
const res = this.matrix({
params: { m: rainfallData.length, n: areas.length, d: 0 },
});
const out = Array.from({ length: maxLength }, (_, i) => 0);
for (let i = 0; i < rainfallData.length; i++) {
for (let j = 0; j < maxLength; j++) {
res[i][j] = rainfallData[i][j] * areas[i];
// Check if totarea is not zero before dividing
out[j] = totalarea !== 0 ? out[j] + res[i][j] / totalarea : out[j];
}
}
return out;
}
/**
* Calculates parameters for the generation of a unit hydrograph
* based on SCS method, Snyder Unit Hydrograph.
* All times of concentrations and lags time are calculated in hours.
* @method syntheticalc
* @memberof hydro
* @param {Object} params - Contains: type (SCS, kerby-kirpich, kerby), unit (si, m).
* @param {Object} args - Contains: l (length), slope (percentage), cn (curve number).
* @param {Object} data - Not used by this function.
* @returns {Object} Calculations depending on type.
* @example
* hydro.analyze.hydro.syntheticalc({params: {type: "SCS", unit: "si"}, args: {l: 4000, slope: 10, cn: 82}});
*/
static syntheticalc({ params, args, data } = {}) {
//imports from parameters.
var type = params.type,
lon = Number(args.l),
sl = Number(args.slope),
units = params.unit,
//Varibles that are to be calculated as solutions.
tc,
tp,
lag,
//Object containing the solutions for the request.
sol = new Object();
if (type === "SCS") {
var sc = 0,
CN = Number(args.cn);
switch (units) {
//longitude in feet, tc in hours, sl in percentage, sc in inches.
case "si":
sc = 1000 / CN - 10;
break;
case "m":
//longitude in meters, tc in hours, sl in percentage, sc in mm.
sc = 25400 / CN - 254;
break;
default:
throw error("Please use a correct unit system!");
}
tc =
(Math.pow(lon, 0.8) * Math.pow(sc + 1, 0.7)) /
(1140 * Math.pow(sl, 0.5));
tp = 0.7 * tc;
lag =
(Math.pow(lon, 0.8) * Math.pow(sc + 1, 0.7)) /
(1900 * Math.pow(sl, 0.5));
Object.assign(sol, {
MaxRetention: sc,
TimeConc: tc,
TimePeak: tp,
LagTime: lag,
});
} else if (type === "kerby-kirpich") {
var K = 0,
M = 0,
N = Number(args.N);
switch (units) {
case "si":
//longitude in feet and sl as number.
K = 0.0078;
M = 1.44;
break;
case "m":
//longitude in meters and sl as number
K = 0.0195;
M = 0.828;
break;
default:
throw error("Please use a correct unit system!");
}
//calculating catchment time
var tov = M * (Math.pow(lon * N), 0.467) * Math.pow(sl / 100, -0.235),
//calculating main channel time
tch = (K * Math.pow(lon, 0.77) * Math.pow(sl / 100, -0.385)) / 60;
//summing both up.
tc = tov + tch;
tp = 0.7 * tc;
lag = 0.6 * tc;
Object.assign(sol, {
TimeConc: tc,
TimePeak: tp,
LagTime: lag,
});
} else if (type === "kerby") {
var n = Number(args.manning);
switch (units) {
// change calculations depending on units.
case "si":
tc = Math.pow((2.2 * n * lon) / Math.pow(sl / 100, 0.5), 0.324) / 60;
break;
case "m":
tc =
(1.4394 * Math.pow((n * lon) / Math.pow(sl / 100, 0.5), 0.467)) /
60;
break;
default:
throw error("Please use a correct unit system!");
}
tp = 0.7 * tc;
lag = 0.6 * tc;
Object.assign(sol, {
TimeConc: tc,
LagTime: lag,
});
}
return sol;
}
/**
* Generates a hydrograph using various distribution types (gamma, LP3, weibull) for a given duration and time step.
* @method dimunithydro
* @memberof hydro
* @param {Object} params - Contains: timeStep (time step between data points), numhours (total duration in hours).
* @param {Object} args - Contains: type (distribution type: "gamma", "lp3", "weibull"), prf (peak reduction factor), lambda (parameter for LP3 distribution), tpeak (peak time for LP3 distribution), alpha (shape parameter for Weibull distribution), beta (scale parameter for Weibull distribution), t0 (location parameter for Weibull distribution).
* @param {Object} data - Not used by this function.
* @returns {Array} Array of two arrays: ttp (time values) and qqp (flow values).
* @example
* hydro.analyze.hydro.dimunithydro({
* params: { timeStep: 0.1, numhours: 10 },
* args: { type: "gamma", prf: 238 }
* });
*/
static dimunithydro({ params, args, data } = {}) {
//populate
var step = Number(params.timeStep),
hours = Number(params.numhours),
//calculate the number of steps in the hydrograph
numstep = Math.round(hours / step) + 1,
//create new vectors
ttp = Array(numstep).fill(0),
qqp = Array(numstep).fill(0),
m = 0;
if (args.type === "gamma") {
//change gamma shape factor.
switch (Number(args.prf)) {
case 101:
m = 0.26;
break;
case 238:
m = 1;
break;
case 349:
m = 2;
break;
case 433:
m = 3;
break;
case 484:
m = 3.7;
break;
case 504:
m = 4;
break;
case 566:
m = 5;
break;
default:
throw Error(
"Please choose value between 101,238,349,433,484,504,566."
);
}
//populating the array with t/tp relationship every 0.1t.
//populating the array with q/qp using Gamma distribution with PRF value.
for (var i = 1; i < ttp.length; i++) {
ttp[i] = Number((ttp[i - 1] + step).toFixed(2));
qqp[i] = Number(
(Math.exp(m) * Math.pow(ttp[i], m) * Math.exp(-m * ttp[i])).toFixed(3)
);
}
return [ttp, qqp];
return [ttp, qqp];
} else if (args.type === "lp3") {
const lambda = Number(args.lambda);
const tpeak = Number(args.tpeak);
//populating the array with t/tp relationship every 0.1t.
//populating the array with q/qp using LP3 distribution with PRF value.
for (var i = 1; i < ttp.length; i++) {
ttp[i] = Number((ttp[i - 1] + step).toFixed(2));
qqp[i] = Number(
(1 / lambda) *
Math.pow((3 / 2) * (ttp[i] / tpeak), lambda) *
Math.exp(-Math.pow(ttp[i] / tpeak, lambda) * (3 / 2))
).toFixed(3);
}
return [ttp, qqp];
} else if (args.type === "weibull") {
var alpha = Number(args.alpha),
beta = Number(args.beta),
t0 = Number(args.t0);
var c1 = (beta / alpha) * Math.exp(-Math.pow(t0 / alpha, beta));
for (var i = 1; i < ttp.length; i++) {
ttp[i] = Number((ttp[i - 1] + step).toFixed(2));
qqp[i] =
c1 *
Math.exp(-Math.pow((ttp[i] - t0) / alpha, beta)) *
(Math.exp((beta / alpha) * ttp[i]) -
Math.exp((beta / alpha) * ttp[i - 1]));
}
return [ttp, qqp];
} else {
throw Error("Please use available distributions!");
}
}
/**
* Hyetograph generator for a uniformly distributed rainfall event.
* This function generates a hyetograph for a long-duration storm based on a uniformly distributed rainfall event.
* @method hyetogen
* @memberof hydro
* @param {Object} params - Contains: duration (number) representing the duration of the storm in hours, timestep (number) representing the time interval for the hyetograph in hours, and intensity (number) representing the rainfall intensity in mm/hour.
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: event (array with rainfall values).
* @returns {Object} An object containing the hyetograph array and the timestep used to generate it.
* @example
* hydro.analyze.hydro.hyetogen({params: {duration: 24, timestep: 1, intensity: 20}});
* hydro.analyze.hydro.hyetogen({data: [1, 2, 3, 4]});
*/
static hyetogen({ params, args, data } = {}) {
const duration = Number(params.duration);
let hyetograph = [];
let timestep = Number(params.timestep);
let intensity;
// If rainfall data is provided, infer hyetograph intensity and generate hyetograph
if (data) {
const rainfall = Array.isArray(data) ? data.flat(Infinity) : [data];
const totalRainfall = rainfall.reduce((a, b) => a + b, 0);
const intensity = totalRainfall / duration;
for (let i = 0; i < rainfall.length; i++) {
hyetograph.push(rainfall[i] / timestep);
}
}
// If intensity is provided, generate hyetograph using normal distribution
else if (params.intensity) {
intensity = Number(params.intensity);
const mean = duration / 2;
const stddev = duration / 6;
const count = Math.round(duration / timestep);
const gaussian = (x) =>
Math.exp(-Math.pow(x - mean, 2) / (2 * Math.pow(stddev, 2)));
const normalize = (x) => x / gaussian(mean);
timestep = duration / count;
for (let i = 0; i < count; i++) {
const x = i * timestep;
const y = normalize(gaussian(x)) * intensity;
hyetograph.push(y);
}
}
return {
hyetograph: hyetograph,
timestep: timestep,
};
}
/**
* Unit hydrograph constructor NRCS constructor depending on the
* physical characteristics of a regularly shaped basin. Adapted from https://directives.sc.egov.usda.gov/OpenNonWebContent.aspx?content=17755.wba
* @method unithydrocons
* @memberof hydro
* @param {Object} params - Contains: drainagearea (in sqmi or km2), type (dim, obs), units(si, m).
* @param {Object} args - Contains: peak (hours), tconcentration (hours), baseflow (cfs or cumecs).
* @param {Object} data - Contains: event either dimensionless or observed as 1d-JS array [[TSevent]].
* @returns {Object[]} Array with time series array. If metric in m3/s, if SI in cfs.
* @example
* hydro.analyze.hydro.unithydrocons({
* params: { drainagearea: 100, type: 'dim', units: 'si' },
* args: { peak: 10, tconcentration: 5, baseflow: 20 },
* data: [[1, 2, 3], [10, 20, 30]]
* });
*/
static unithydrocons({ params, args, data } = {}) {
//import parameters from user.
var area = Number(params.drainagearea),
duh = data,
unit = this.matrix({ params: { m: 2, n: duh[0].length, d: 0 } });
//unit hydro from dimensionless hydrograph.
if (params.type == "dim") {
//peak rate factor chosen.
var peak = Number(args.peak),
//calculate time step.
tconc = Number(args.tconcentration),
deltat = Number((tconc * 0.133).toFixed(3)),
//calculate time to peak and construct result arrays.
tp = deltat / 2 + 0.6 * tconc,
qp = 0;
//change peak discharge depending on the units.
switch (params.units) {
case "si":
qp = (peak * area * 1) / tp;
break;
case "m":
qp = (peak * area * 1) / tp;
break;
default:
alert("Please input a valid unit system!");
}
//populate the hydrograph with time and discharge.
for (var h = 0; h < unit[0].length; h++) {
unit[0][h] = Number((duh[0][h] * tp).toFixed(3));
unit[1][h] = Number((duh[1][h] * qp).toFixed(3));
}
return unit;
}
//unit hydro from observed hydrograph.
else if (params.type == "obs") {
var baseflow = Number(args.baseflow),
drh = this.matrix({ params: { m: 1, n: duh[0].length, d: 0 } });
unit[0] = duh[0];
//timestep in hours
var timestep = Math.abs(unit[0][1] - unit[0][0]) * 60 * 60;
console.log(unit);
for (var i = 0; i < unit[0].length; i++) {
drh[i] = Math.abs(duh[1][i] - baseflow);
}
var sum = this.totalprec({ data: drh }) * timestep,
vol = 0;
switch (params.units) {
case "si":
//calculated in inches
vol = Math.round((sum / area) * 12);
break;
case "m":
//calculated in cms
vol = Math.round((sum / area) * 100);
}
for (var j = 0; j < unit[0].length; j++) {
//unit hydrograph in cfs/inch or cumecs/cm
unit[1][j] = Math.round(drh[j] / vol);
}
unit[1].reverse();
return {
unithydro: unit,
totalvol: vol,
};
}
}
/**
* Flooding hydrograph generator using a unit hydrograph,
* precipitation data and SCS metrics for runoff calculation.
* If the observed hydrograph option is selected, the precipitation must be dividied in
* blocks of rainfall in as a 2D array [[date, date, date], [rainf, rainf, rainf]].
* @method floodhydro
* @memberof hydro
* @param {Object} params - Contains: baseflow.
* @param {Object} args - Contains: type (SCS, obs), CN (if SCS), stormduration (event duration in hours), timestep (in hours), units (si, m).
* @param {Object[]} data - Contains: 2d-JS array with Timeseries Data [[TSrainfall], [TSunithydro]].
* @returns {Object[]} Array with values for runoff as time series.
* @example
* hydro.analyze.hydro.floodhydro({
* params: { baseflow: 10 },
* args: { type: 'obs', CN: 80, stormduration: 24, timestep: 1, units: 'si' },
* data: [[[1, 2, 3], [0.1, 0.2, 0.4]], [[1, 2, 3], [0.3, 0.1, 0.2]]]
* });
*/
/**
* Generate flood hydrograph using SCS method or observed unit hydrograph
* Converts rainfall excess to runoff hydrograph using convolution
* Essential for flood forecasting and stormwater design
*
* @function floodhydro
* @memberof hydro
* @param {Object} options - Function options
* @param {Object} [options.params] - Parameters
* @param {number} [options.params.baseflow=0] - Baseflow to add to hydrograph (m³/s)
* @param {Object} options.args - Arguments
* @param {string} options.args.type - Method: 'SCS' or 'obs' (observed unit hydrograph)
* @param {number} [options.args.cn] - SCS Curve Number (for SCS method, 0-100)
* @param {number} [options.args.stormduration] - Storm duration (hours, for SCS)
* @param {number} [options.args.timestep] - Time step (hours, for SCS)
* @param {string} [options.args.units] - Unit system: 'si' or 'm' (for SCS)
* @param {Array} options.data - [rainfall, unitHydrograph]
* @returns {Array} Flood hydrograph [[time], [discharge]]
*
* @example
* // SCS Method: Generate flood hydrograph for storm event
* const rainfall = [
* [0, 1, 2, 3, 4, 5], // Time (hours)
* [0, 5, 15, 8, 3, 0] // Rainfall (mm/hr)
* ];
*
* const unitHydrograph = [
* [0, 1, 2, 3, 4, 5, 6, 7, 8], // Time (hours)
* [0, 10, 25, 35, 30, 20, 10, 5, 0] // Unit discharge (m³/s/cm)
* ];
*
* const floodHydrograph = hydro.analyze.hydro.floodhydro({
* params: { baseflow: 5 }, // 5 m³/s baseflow
* args: {
* type: 'SCS',
* cn: 75, // Curve Number for urban area
* stormduration: 6, // 6-hour storm
* timestep: 1, // 1-hour timestep
* units: 'si' // SI units
* },
* data: [rainfall, unitHydrograph]
* });
*
* // Result: [time array, discharge array]
* console.log(`Peak flow: ${Math.max(...floodHydrograph[1])} m³/s`);
*
* @example
* // Observed Unit Hydrograph Method
* const floodHydrograph = hydro.analyze.hydro.floodhydro({
* params: { baseflow: 3 },
* args: { type: 'obs' },
* data: [rainfall, observedUnitHydrograph]
* });
*
* @example
* // Design storm for 100-year flood
* // CN varies by land use: Forest=55, Agriculture=70, Urban=85
* const designStorm = [[0,1,2,3,4], [2,8,20,12,5]]; // 100-yr rainfall
* const syntheticUH = /* from Snyder or other sources;
const designFlood = hydro.analyze.hydro.floodhydro({
* args: { type: 'SCS', cn: 80, stormduration: 4, timestep: 1, units: 'si' },
* data: [designStorm, syntheticUH]
* });
*/
static floodhydro({ params = {}, args = {}, data = [] } = {}) {
const [rain, unit] = data;
let { baseflow = 0 } = params;
baseflow = Number(baseflow);
if (args.type == "SCS") {
const cn = Number(args.cn),
stormdur = Number(args.stormduration),
timestep = Number(args.timestep);
//transform the date into javascript format.
//create arrays for calculation of runoff
var numarray = Math.round(stormdur / timestep),
finalcount = numarray + unit[0].length,
sc = 0,
accumrainf = this.matrix({
params: { m: 2, n: rain[1].length, d: 0 },
});
accumrainf[0] = rain[0];
var accumrunff = this.matrix({
params: { m: 2, n: rain[1].length, d: 0 },
});
accumrunff[0] = rain[0];
var incrementrunff = this.matrix({
params: { m: 2, n: rain[1].length, d: 0 },
});
incrementrunff[0] = rain[0];
var hydros = this.matrix({
params: { m: stormdur, n: finalcount, d: 0 },
}),
finalhydro = this.matrix({ params: { m: 2, n: finalcount, d: 0 } });
// change calculations depending on units.
switch (args.units) {
case "si":
sc = 1000 / cn - 10;
break;
case "m":
sc = 25400 / cn - 254;
break;
default:
alert("Please use a correct unit system!");
}
//add accumulative rainfall and calculate initial abstraction.
var iniabs = 0.2 * sc;
rain[1]
.slice()
.reduce((prev, curr, i) => (accumrainf[1][i] = prev + curr), 0);
//add runoff calculations.
for (var i = 0; i < numarray; i++) {
if (accumrainf[1][i] > 0) {
accumrunff[1][i] =
Math.pow(accumrainf[1][i] - iniabs, 2) /
(accumrainf[1][i] - iniabs + sc);
} else {
accumrunff[1][i] = 0;
}
incrementrunff[1][i] = Number(
(Math.abs(accumrunff[1][i] - accumrunff[1][i - 1]) || 0).toFixed(3)
);
}
//create composite hydrograph.
for (var j = 0; j < hydros[0].length; j++) {
hydros[0][j] = Number(
(incrementrunff[1][j] * unit[1][j] || 0).toFixed(3)
);
finalhydro[0][j] = Number(finalhydro[0][j] + timestep);
}
//populate the moving hydrographs
for (var h = 1; h < hydros.length; h++) {
for (var k = 0; k < hydros[0].length; k++) {
hydros[h][k + h] = Number(hydros[0][k].toFixed(3));
finalhydro[1][k] = Number(hydros[h][k].toFixed(3)) + baseflow;
}
}
//accumulate timespan for cumulative hydrograph.
finalhydro[0]
.slice()
.reduce(
(prev, curr, i) =>
(finalhydro[0][i] = Number((prev + curr).toFixed(2), 0))
);
for (var p = 0; p < finalhydro[1].length; p++) {
finalhydro[1][p] = finalhydro[1][p];
}
finalhydro[1].reverse();
return finalhydro;
} else if (args.type == "obs") {
var hydros = [],
timestep = Math.abs(rain[0][1] - rain[0][0]);
//calculate the runoff per pulse.
for (var i = 0; i < rain[1].length; i++) {
var neq = [];
for (var j = 0; j < unit[1].length - 1; j++) {
neq.push(unit[1][j] * rain[1][i]);
}
hydros.push(neq);
}
var final = this.matrix({
params: { m: 2, n: unit[1].length + hydros.length, d: 0 },
});
//zeros up
for (var k = 0; k < hydros.length; k++) {
var zeros = new Array(timestep * hydros.indexOf(hydros[k])).fill(0);
zeros.forEach((x) => hydros[k].unshift(x));
hydros[k].shift();
}
//zeros down
for (var l = 0; l < hydros.length; l++) {
var finalarr = hydros[hydros.length - 1].length,
zeros = new Array(finalarr - hydros[l].length).fill(0);
zeros.forEach((x) => hydros[l].push(x));
}
final[1] = hydros[0].map((x, i) =>
hydros.reduce((sum, curr) => sum + curr[i], baseflow)
);
//time and discharge sum
for (var p = 0; p < final[1].length; p++) {
final[0][p] = p;
}
return final;
}
}
/**
* Simple rainfall-runoff analyses over a rainfall dataset given landuse, baseflow and infiltration capacity. It is mainly used for long-term hydrological analysis such as monthly changes.
* All should be in mm. For more info, refer to https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.544.5085&rep=rep1&type=pdf
* @method bucketmodel
* @memberof hydro
* @param {Object} params - Contains: baseflow, infiltration (in same units).
* @param {Object} args - Contains: landuse percentage (summing all to 1): agriculture, barerock, grassland, forest, urban.
* @param {Object} data - Contains: rainfall and evaporation arrays.
* @returns {Array} 1d-array with time series according to different time spans (5min, 15 min, 1 hour, 1 day...).
* @example
* hydro.analyze.hydro.bucketmodel({
* params: { baseflow: 1, infiltration: 0.3 },
* args: { agriculture: 0.1, barerock: 0.5, grassland: 0.1, forest: 0.05, urban: 0.05 },
* data: { rainfall: [1, 2, 3, 4, 5], evaporation: [0.1, 0.2, 0.3, 0.4, 0.5] }
* });
*/
static bucketmodel({ params, args, data } = {}) {
// Initial parameters
const { rainfall, evaporation } = data;
let { baseflow, infiltration } = params;
// Robust parsing
baseflow = Number(baseflow);
infiltration = Number(infiltration);
let { agriculture, barerock, grassland, forest, urban } = args;
agriculture = Number(agriculture || 0);
barerock = Number(barerock || 0);
grassland = Number(grassland || 0);
forest = Number(forest || 0);
urban = Number(urban || 0);
const n = rainfall.length;
const landuse = [agriculture, barerock, grassland, forest, urban];
baseflow = baseflow / 86400; // Convert from m3/s to mm/s
// Infiltration capacities for agriculture, bare rock, grassland, forest and
// urban, respectively in mm.
const FieldCaps = [25, 5, 25, 50, 5];
// Arrays and variables
const initial = this.matrix({ params: { m: landuse.length, n: n, d: 0 } });
const interflow = this.matrix({
params: { m: landuse.length, n: n, d: 0 },
});
const overflow = this.matrix({ params: { m: landuse.length, n: n, d: 0 } });
const totalflow = this.matrix({
params: { m: landuse.length, n: n, d: 0 },
});
const totalrunoff = new Array(n).fill(0);
// Initial moisture
for (let i = 0; i < FieldCaps.length; i++) {
initial[i][0] = FieldCaps[i] * landuse[i] + rainfall[0] - evaporation[0];
}
// Initial soil moisture
for (let k = 0; k < FieldCaps.length; k++) {
if (initial[k][0] > FieldCaps[k]) {
overflow[k][0] = initial[k][0] - FieldCaps[k];
initial[k][0] = FieldCaps[k];
} else {
overflow[k][0] = 0;
}
if (initial[k][0] > 0) {
interflow[k][0] = initial[k][0] * infiltration;
} else {
interflow[k][0] = 0;
}
}
// Calculating overland and interflow
for (let m = 0; m < FieldCaps.length; m++) {
for (let p = 1; p < n; p++) {
initial[m][p] =
initial[m][p - 1] * (1 - infiltration) + rainfall[p] - evaporation[p];
if (initial[m][p] > FieldCaps[m]) {
overflow[m][p] = initial[m][p] - FieldCaps[m];
initial[m][p] = FieldCaps[m];
} else {
overflow[m][p] = 0;
}
if (initial[m][p] > 0) {
interflow[m][p] = initial[m][p] * infiltration;
} else {
interflow[m][p] = 0;
}
}
}
// Calculating the total amount of flow from overflow, baseflow and interflow
for (let j = 0; j < FieldCaps.length; j++) {
for (let h = 0; h < n; h++) {
totalflow[j][h] =
overflow[j][h] + interflow[j][h] + baseflow * landuse[j];
}
}
//calculating total runoff
for (var q = 0; q < n; q++) {
totalrunoff[q] =
totalflow[0][q] * landuse[0] +
totalflow[1][q] * landuse[1] +
totalflow[2][q] * landuse[2] +
totalflow[3][q] * landuse[3] +
totalflow[4][q] * landuse[4];
}
/*var finalvalues = this.matrix(2,n, 0)
for (var w = 0; w < finalvalues[1].length; w++) {
finalvalues[0][w] = w * 60;
finalvalues[1][w] = totalrunoff[w];
}
var agg = this.rainaggr({"event": finalvalues, "agg": {"type": "aggr", "interval": 1440}})*/
return totalrunoff;
}
/**
* Solves 1d groundwater steady simulation using gaussian elimination for a static setting.
* Adapted from (Molkentin, 2019).
* @method staticGround1d
* @memberof hydro
* @param {Object} params - Contains: length, k (discharge coeff), nodes.
* @param {Object} args - Contains: w0 and q0 (extraction, discharge at point 0), w1 and q1 (extraction, discharge point 1).
* @param {Object} data - Not used by this function.
* @return {Array} Matrix with solutions.
* @example
* hydro.analyze.hydro.staticGround1d({
* params: { length: 100, k: 0.5, nodes: 10 },
* args: { w0: 0, w1: 0, q0: 1, qL: 1 }
* });
*/
static staticGround1d({ params, args, data } = {}) {
// Robust parsing
const length = Number(params.length);
const k = Number(params.k);
const nodes = Number(params.nodes);
// args might have defaults or be optional depending on implementation logic, but parsing if present
const w0 = Number(args.w0 || 0);
const w1 = Number(args.w1 || 0);
const hL = Number(args.hL || 0); // Added hL parsing
const q0 = Number(args.q0 || 0);
const qL = Number(args.qL || 0);
const dx = length / (nodes - 1);
//create a new equation system
let matrix = this.matrix({ params: { m: nodes, n: nodes, d: 0 } }),
vec_left = this.matrix({ params: { m: 1, n: nodes, d: 0 } }),
vec_right = this.matrix({ params: { m: 1, n: nodes, d: 0 } }),
//equation system set up.
factor = k / dx;
//initial boundary.
matrix[0][0] = factor;
matrix[0][1] = -factor;
vec_right[0] = q0;
//last boundary.
var index = nodes - 1;
matrix[index][index] = -factor;
matrix[index][index - 1] = factor;
vec_right[index] = qL;
//calculate inner nodes using Runge-Kutta 1D.
for (var i = 1; i < nodes - 1; i++) {
var newfac = k / (dx * dx);
matrix[i][i] = 2 * newfac;
matrix[i][i - 1] = -1 * newfac;
matrix[i][i + 1] = -1 * newfac;
vec_right[i] = w0 + w1 * i * dx;
}
//solve equation system.
this.equationsystemsolver({
data: matrix,
params: { left: vec_left, right: vec_right },
});
//calculate discharge in vec_right
vec_right[0] = (-k * (vec_left[1] - vec_left[0])) / dx;
for (var i = 1; i < nodes - 1; i++) {
vec_right[i] = ((-k * (vec_left[i + 1] - vec_left[i - 1])) / dx) * 0.5;
}
vec_right[index] = (-k * (vec_left[index] - vec_left[index - 1])) / dx;
return vec_left;
}
/**
* Solves 1D dynamic groundwater simulation using the Crank-Nicolson method.
* Adapted from (Molkentin, 2019).
* @method dynamicGround1D
* @memberof hydro
* @param {Object} params - Contains: length (length of the domain), nodes (number of nodes), k (hydraulic conductivity).
* @param {Object} args - Contains: dt (time step), T (total simulation time), h0 (initial hydraulic head), hL (hydraulic head at the boundary), q0 (flow rate at the boundary), qL (flow rate at the boundary), phi (porosity), Ss (specific storage), Sy (specific yield).
* @param {Object} data - Not used by this function.
* @returns {Array} Array of solutions representing the hydraulic head at each node.
* @example
* hydro.analyze.hydro.dynamicGround1D({
* params: { length: 100, nodes: 10, k: 0.5 },
* args: { dt: 0.1, T: 10, h0: 10, hL: 5, q0: 1, qL: 0.5, phi: 0.3, Ss: 0.002, Sy: 0.15 }
* });
*/
static dynamicGround1D({ params, args, data } = {}) {
// pass data from params to variables
const { length, nodes, k } = params;
const { dt, T, h0, hL, q0, qL, phi, Ss, Sy } = args;
const dx = length / (nodes - 1);
const a = k / (phi * Ss);
const b = (Sy * a) / dt;
const c = k / dx;
const d = Sy * h0 + (1 - Sy) * hL;
// create a new equation system
const matrix = this.matrix({ params: { m: nodes, n: nodes, d: 0 } });
const vec_left = this.matrix({ params: { m: 1, n: nodes, d: 0 } });
const vec_right = this.matrix({ params: { m: 1, n: nodes, d: 0 } });
// equation system set up
matrix[0][0] = 1;
vec_left[0] = h0;
matrix[nodes - 1][nodes - 1] = 1;
vec_left[nodes - 1] = hL;
// calculate inner nodes using Crank-Nicolson method for each time step
const numTimeSteps = Math.round(T / dt);
for (let t = 0; t < numTimeSteps; t++) {
for (let i = 1; i < nodes - 1; i++) {
const aTerm = a / dx ** 2;
const bTerm = b / 2;
const cTerm = c / 2;
const laplacian =
(vec_left[i + 1] - 2 * vec_left[i] + vec_left[i - 1]) * aTerm;
const e = (Sy * vec_left[i] + (1 - Sy) * hL + b * q0 - laplacian) / b;
matrix[i][i - 1] = -cTerm;
matrix[i][i] = 1 / dt + cTerm + 2 * aTerm;
matrix[i][i + 1] = -cTerm;
vec_right[i] = vec_left[i] / dt + q0 * cTerm + laplacian + bTerm * e;
}
// solve equation system for this time step
this.equationsystemsolver({
data: matrix,
params: { left: vec_left, right: vec_right },
});
// update boundary conditions for this time step
vec_right[0] = q0;
for (let i = 1; i < nodes - 1; i++) {
vec_right[i] = (-k * (vec_left[i + 1] - vec_left[i - 1])) / (2 * dx);
}
vec_right[nodes - 1] = qL;
}
return vec_left;
}
/**
* Aggregates or dissaggregates rainfall data depending on what
* the user requires. The date type must be a Javascript string or number and in minutes or hours, but both
* the aggregation interval require and the data interval should be the same.
* For aggregation, the interval for aggregation must be larger than the time step. For example,
* 15 min or 30 min data to be aggregatted every 60 minutes. Intervals must be multiples of the final aggregaiton (2, 4, etc).
* @method rainaggr
* @memberof hydro
* @param {Object} params - Contains: type (aggr, disagg), interval (in minutes for both cases).
* @param {Object} args - Not used by this function.
* @param {Object[]} data - Contains: data as 2D array.
* @returns {Object[]} Array with aggregated/disaggregated data.
* @example
* hydro.analyze.hydro.rainaggr({params: {type: 'aggr', interval: 240}, data: [[1, 2, 3], [0.1, 0.2, 0.3]]});
*/
static rainaggr({ params, args, data } = {}) {
const isArray1D = Array.isArray(data) && typeof data[0] === "number";
const agtype = params.type;
const finagg = params.interval; // in minutes
let datetr;
// Handle single rainfall array
if (isArray1D) {
const dummyTime = data.map((_, i) => i * finagg * 60 * 1000); // fake uniform timestamp
datetr = [dummyTime, data];
} else {
// assume [time, values] format
datetr = this.matrix({
params: { m: data.length, n: data[1].length, d: 0 },
});
const isTimeString = typeof data[0][0] === "string";
for (let i = 0; i < data[0].length; i++) {
datetr[0][i] = isTimeString ? Date.parse(data[0][i]) : data[0][i];
}
for (let j = 1; j < data.length; j++) {
datetr[j] = data[j];
}
}
if (agtype === "aggr") {
const timeSeries = datetr[0];
const rainSeries = datetr[1];
const timestep = Math.abs(timeSeries[1] - timeSeries[0]) / (60 * 1000); // in minutes
const count = Math.round(finagg / timestep);
const totalSteps = Math.floor(timeSeries.length / count);
const aggTime = [];
const aggData = [];
for (let i = 0; i < totalSteps; i++) {
const start = i * count;
const end = start + count;
const segmentTimes = timeSeries.slice(start, end);
const segmentData = rainSeries.slice(start, end);
const timeVal = new Date(segmentTimes[0]).toISOString();
const rainVal = this.totalprec({ data: segmentData });
aggTime.push(timeVal);
aggData.push(rainVal);
}
return [aggTime, aggData];
}
else if (agtype === "disagg") {
const method = params.method || "uniform";
const factor = params.factor || 2; // how many substeps per original timestep
const originalTimes = datetr[0];
const originalVals = datetr[1];
const disaggTime = [];
const disaggVals = [];
const generateSubtimes = (start, interval, steps) => {
const dt = interval / steps;
return Array.from({ length: steps }, (_, i) => new Date(start + i * dt).toISOString());
};
const disaggregateValue = (val, method, steps) => {
switch (method) {
case "uniform":
return Array(steps).fill(val / steps);
case "normal": {
const mean = steps / 2;
const stdDev = steps / 6;
const weights = Array.from({ length: steps }, (_, i) =>
Math.exp(-0.5 * Math.pow((i - mean) / stdDev, 2))
);
const sum = weights.reduce((a, b) => a + b, 0);
return weights.map(w => (w / sum) * val);
}
case "exponential": {
const lambda = 1;
const raw = Array.from({ length: steps }, () => -Math.log(Math.random()) / lambda);
const total = raw.reduce((a, b) => a + b, 0);
return raw.map(x => (x / total) * val);
}
case "huff": {
// Huff Type II 2nd quartile distribution (most intense early-mid period)
const basePattern = [0.06, 0.10, 0.14, 0.20, 0.17, 0.13, 0.10, 0.07, 0.02, 0.01];
const steps = Math.min(steps, basePattern.length);
const pattern = basePattern.slice(0, steps);
const sum = pattern.reduce((a, b) => a + b, 0);
return pattern.map(p => (p / sum) * val);
}
case "gumbel": {
// Gumbel distribution (for extremes)
const u = 0.5772; // Euler-Mascheroni constant
const beta = 1;
const alpha = 0;
const x = Array.from({ length: steps }, (_, i) => i + 1);
const weights = x.map(t => Math.exp(-Math.exp(-(t - alpha) / beta)));
const total = weights.reduce((a, b) => a + b, 0);
return weights.map(w => (w / total) * val);
}
default:
throw new Error("Unknown disaggregation method: " + method);
}
};
for (let i = 0; i < originalVals.length; i++) {
const t0 = typeof originalTimes[i] === "string"
? Date.parse(originalTimes[i])
: originalTimes[i];
const interval = i < originalTimes.length - 1
? (Date.parse(originalTimes[i + 1]) - t0)
: (originalTimes[i] - originalTimes[i - 1]);
const subtimes = generateSubtimes(t0, interval, factor);
const subdivided = disaggregateValue(originalVals[i], method, factor);
disaggTime.push(...subtimes);
disaggVals.push(...subdivided);
}
return [disaggTime, disaggVals];
}
throw new Error("Invalid aggregation type specified");
}
/**
* Calculate reference evapotranspiration using FAO-56 Penman-Monteith equation
* Gold standard method for ET estimation, accounts for energy balance and aerodynamic factors
* Requires weather station data: radiation, temperature, wind, humidity
*
* @function ETPenmanMontheith
* @memberof hydro
* @author riya-patil
* @param {Object} options - Function options
* @param {Object} options.params - Weather parameters
* @param {number|Array<number>} options.params.netRadiation - Net radiation at crop surface (MJ/m²/day)
* @param {number|Array<number>} options.params.temperature - Mean daily air temperature (°C)
* @param {number|Array<number>} options.params.windSpeed - Wind speed at 2m height (m/s)
* @param {number|Array<number>} options.params.saturationVaporPressure - Saturation vapor pressure (kPa)
* @param {number|Array<number>} options.params.actualVaporPressure - Actual vapor pressure (kPa)
* @param {number|Array<number>} options.params.airPressure - Atmospheric pressure (kPa, default ~101.3 at sea level)
* @param {number|Array<number>} options.params.psychrometricConstant - Psychrometric constant (kPa/°C, typically ~0.065)
* @param {Object} [options.args] - Additional arguments (not used)
* @param {Object} [options.data] - Additional data (not used)
* @returns {number|Array<number>} Reference evapotranspiration (mm/day)
* @throws {Error} If required parameters are missing
*
* @example
* // Single day ET calculation
* const et0 = hydro.analyze.hydro.ETPenmanMontheith({
* params: {
* netRadiation: 15.5, // MJ/m²/day
* temperature: 25, // °C
* windSpeed: 2.0, // m/s
* saturationVaporPressure: 3.17, // kPa (at 25°C)
* actualVaporPressure: 2.1, // kPa (66% RH)
* airPressure: 101.3, // kPa (sea level)
* psychrometricConstant: 0.067 // kPa/°C
* }
* });
* console.log(`Reference ET: ${et0.toFixed(2)} mm/day`);
*
* @example
* // Time series ET calculation
* const dailyET = hydro.analyze.hydro.ETPenmanMontheith({
* params: {
* netRadiation: [12.3, 14.5, 15.8, 13.2],
* temperature: [22, 24, 26, 23],
* windSpeed: [1.5, 2.0, 2.5, 1.8],
* saturationVaporPressure: [2.64, 2.98, 3.36, 2.81],
* actualVaporPressure: [1.85, 2.09, 2.35, 1.97],
* airPressure: [101.3, 101.3, 101.3, 101.3],
* psychrometricConstant: [0.067, 0.067, 0.067, 0.067]
* }
* });
* console.log(`Total ET over period: ${dailyET.reduce((a,b) => a+b).toFixed(1)} mm`);
*
* @example
* // Reference: FAO-56 Paper
* // http://www.fao.org/3/X0490E/x0490e06.htm
* // ET0 = (0.408Δ(Rn-G) + γ(900/(T+273))u2(es-ea)) / (Δ + γ(1+0.34u2))
*/
static ETPenmanMontheith({ params, args, data } = {}) {
const { netRadiation, temperature, windSpeed, saturationVaporPressure, actualVaporPressure, airPressure, psychrometricConstant } = params;
// Validate required parameters
if (!temperature || !netRadiation || !windSpeed || !saturationVaporPressure || !actualVaporPressure || !airPressure || !psychrometricConstant) {
throw new Error('Missing required parameters: temperature, netRadiation, windSpeed, saturationVaporPressure, actualVaporPressure, airPressure, psychrometricConstant.');
}
const evapotranspiration = [];
for (let i = 0; i < temperature.length; i++) {
const numerator = (0.408 * netRadiation[i] + 0.063 * 900 * (temperature[i] + 273) * windSpeed[i] * (saturationVaporPressure[i] - actualVaporPressure[i]));
const denominator = lambda * (0.408 * (netRadiation[i] - 0) + rho * (900 / (temperature[i] + 273)) * windSpeed[i] * (saturationVaporPressure[i] - actualVaporPressure[i]));
const et = (numerator / denominator) * 1000; // Convert from m/day to mm/day
evapotranspiration.push(et);
}
return evapotranspiration;
}
/**
* Calculate evapotranspiration using Hargreaves method
* Simpler alternative to Penman-Monteith, requires only temperature data
* Good for data-scarce regions where only temperature is available
*
* @function ETHargreaves
* @memberof hydro
* @author riya-patil
* @param {Object} options - Function options
* @param {Object} options.params - Parameters
* @param {number} options.params.latitude - Latitude in decimal degrees
* @param {Object} [options.args] - Additional arguments (not used)
* @param {Object} options.data - Input data
* @param {Array<number>} options.data.temperature - Mean daily temperature (°C)
* @param {Array<number>} options.data.temperatureMax - Maximum daily temperature (°C)
* @param {Array<number>} options.data.temperatureMin - Minimum daily temperature (°C)
* @param {Array<string>} options.data.date - Dates (ISO format or parseable)
* @returns {Array<number>} Daily evapotranspiration (mm/day)
* @throws {Error} If required parameters are missing
*
* @example
* // Daily ET calculation for a week
* const et = hydro.analyze.hydro.ETHargreaves({
* params: { latitude: 40.5 }, // Boulder, CO
* data: {
* temperature: [18, 20, 22, 19, 21, 23, 20],
* temperatureMax: [25, 27, 29, 26, 28, 30, 27],
* temperatureMin: [11, 13, 15, 12, 14, 16, 13],
* date: [
* '2024-07-01', '2024-07-02', '2024-07-03',
* '2024-07-04', '2024-07-05', '2024-07-06', '2024-07-07'
* ]
* }
* });
* console.log(`Weekly ET total: ${et.reduce((a,b) => a+b).toFixed(1)} mm`);
*
* @example
* // Growing season ET estimation
* const growingSeasonET = hydro.analyze.hydro.ETHargreaves({
* params: { latitude: 35.0 },
* data: {
* temperature: dailyTemps, // 120 days
* temperatureMax: dailyMaxTemps,
* temperatureMin: dailyMinTemps,
* date: dateArray
* }
* });
* const totalET = growingSeasonET.reduce((sum, et) => sum + et, 0);
* console.log(`Total growing season ET: ${totalET.toFixed(0)} mm`);
*
* @example
* // Compare with precipitation for water balance
* const precipitation = [5, 0, 0, 12, 8, 0, 3]; // mm/day
* const et = hydro.analyze.hydro.ETHargreaves({
* params: { latitude: 40.5 },
* data: { temperature, temperatureMax, temperatureMin, date }
* });
*
* const waterBalance = precipitation.map((p, i) => p - et[i]);
* console.log('Daily water balance (P - ET):', waterBalance);
*/
static ETHargreaves({ params, args, data } = {}) {
const { latitude } = params;
const { temperature, temperatureMax, temperatureMin, date } = data;
// Validate required parameters
if (!temperature || !temperatureMax || !temperatureMin || !date) {
throw new Error('Missing required parameters: temperature, temperatureMax, temperatureMin, date.');
}
// Calculate evapotranspiration for each time step
const evapotranspiration = [];
for (let i = 0; i < temperature.length; i++) {
const julianDay = getJulianDay(date[i]);
const Ra = 4.903 * Math.pow(10, -9); // Extraterrestrial radiation constant (MJ/(m^2 * min * °C))
const dr = 1 + 0.033 * Math.cos((2 * Math.PI / 365) * julianDay); // Inverse relative distance Earth-Sun
const delta = 0.409 * Math.sin((2 * Math.PI / 365) * julianDay - 1.39); // Solar declination angle
const TmaxK = temperatureMax[i] + 273; // Convert to Kelvin
const TminK = temperatureMin[i] + 273;
const RaTmax = (24 * 60 / Math.PI) * Ra * dr * (TmaxK + TminK) * (Math.sin((latitude * Math.PI) / 180) * Math.sin(delta) + Math.cos((latitude * Math.PI) / 180) * Math.cos(delta) * Math.sin(0)); // MJ/m^2/day
const Rs = 0.16 * RaTmax; // Convert to MJ/m^2/day
const et = 0.0023 * (temperature[i] + 17.8) * Math.sqrt(temperatureMax[i] - temperatureMin[i]) * Rs;
evapotranspiration.push(et);
}
return evapotranspiration;
}
/**
* Calculates evapotranspiration using the Thornthwaite method.
* Reference: https://wikifire.wsl.ch/tiki-indexf125.html?page=Potential+evapotranspiration
* @method ETThornthwaite
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: latitude (latitude in decimal degrees).
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: temperature (mean monthly air temperature in °C), monthDays (number of days in each month).
* @returns {Number[]} Evapotranspiration in mm/day for each month.
* @throws {Error} if missing required data, invalid data format (not in array), or unequal array length.
* @example
* hydro.analyze.hydro.ETThornthwaite({
* params: { latitude: 40 },
* data: {
* temperature: [10, 15, 20],
* monthDays: [31, 28, 31]
* }
* });
*/
static ETThornthwaite({ params, args, data } = {}) {
const { latitude } = params;
const { temperature, monthDays } = data;
if (!temperature || !monthDays) {
throw new Error('Missing required data: temperature, monthDays.');
}
// Validate temperature and monthDays arrays
if (!Array.isArray(temperature) || !Array.isArray(monthDays)) {
throw new Error('Invalid data format. Expected temperature and monthDays to be arrays.');
}
if (temperature.length !== monthDays.length) {
throw new Error('Temperature and monthDays arrays must have the same length.');
}
// Calculate heat index (HI) for each month
const hiValues = temperature.map((t) => {
const hi = (t / 5) ** 1.514;
return hi > 0 ? hi : 0;
});
const petValues = hiValues.map((hi, i) => {
const monthDaysValue = monthDays[i];
const pet = (0.6 * (latitude / 5) * Math.pow(10, (0.514 * hi))) * monthDaysValue;
return pet;
});
return petValues;
}
/**
* Calculates evapotranspiration using the Blaney-Criddle method.
* Reference: https://legacy.azdeq.gov/environ/water/permits/download/blaney.pdf
* @method ETBlaneyCriddle
* @author riya-patil
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: temperature (mean monthly air temperature in °C), monthDays (number of days in each month).
* @returns {Number[]} Evapotranspiration in mm/day for each month.
* @throws {Error} if missing data, data not in array format, unequal length of arrays.
* @example
* hydro.analyze.hydro.ETBlaneyCriddle({
* data: {
* temperature: [10, 15, 20],
* monthDays: [31, 28, 31]
* }
* });
*/
static ETBlaneyCriddle({ params, args, data } = {}) {
const { temperature, monthDays } = data;
if (!temperature || !monthDays) {
throw new Error('Missing required data: temperature, monthDays.');
}
// Validate temperature and monthDays arrays
if (!Array.isArray(temperature) || !Array.isArray(monthDays)) {
throw new Error('Invalid data format. Expected temperature and monthDays to be arrays.');
}
if (temperature.length !== monthDays.length) {
throw new Error('Temperature and monthDays arrays must have the same length.');
}
// Calculate monthly potential evapotranspiration (PET) using Blaney-Criddle equation
const petValues = temperature.map((t, i) => {
const monthDaysValue = monthDays[i];
const pet = (0.02 * (t + 17.8)) * monthDaysValue;
return pet;
});
return petValues;
}
/**
* Calculates evapotranspiration using the Priestley-Taylor method.
* Reference: https://wetlandscapes.github.io/blog/blog/penman-monteith-and-priestley-taylor/
* @method ETPriestelyTaylor
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: netRadiation (W/m^2), latentHeatFlux (W/m^2).
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @returns {Number} Evapotranspiration in mm/day.
* @example
* hydro.analyze.hydro.ETPriestelyTaylor({ params: { netRadiation: 3, latentHeatFlux: 3 } });
*/
static ETPriestelyTaylor({ params, args, data } = {}) {
const { netRadiation, latentHeatFlux } = params;
// Calculate potential evapotranspiration using Priestley-Taylor method
const evapotranspiration = 1.26 * (netRadiation / latentHeatFlux);
return evapotranspiration;
}
/**
* Calculates infiltration using the Green-Ampt model.
* Reference: https://www.hec.usace.army.mil/confluence/rasdocs/ras1dtechref/6.1/overview-of-optional-capabilities/modeling-precipitation-and-infiltration/green-ampt
* @method InfGreenAmpt
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: Ks (saturated hydraulic conductivity [L/T]), psi (soil suction head [L]), theta_i (initial soil moisture content [L^3/L^3]), theta_s (saturated soil moisture content [L^3/L^3]), t (time [T]).
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @returns {Number} Infiltration rate [L/T].
* @throws {Error} invalid data type is inputted.
* @example
* hydro.analyze.hydro.InfGreenAmpt({ params: { Ks: 10, psi: 5, theta_i: 0.2, theta_s: 0.4, t: 2 } });
*/
static InfGreenAmpt({ params, args, data } = {}) {
const { Ks, psi, theta_i, theta_s, t } = params;
// Validate data inputs
if (
typeof Ks !== 'number' ||
typeof psi !== 'number' ||
typeof theta_i !== 'number' ||
typeof theta_s !== 'number' ||
typeof t !== 'number'
) {
throw new Error('Invalid data inputs. Expected numbers for Ks, psi, theta_i, theta_s, and t.');
}
// Calculate infiltration using the Green-Ampt method
const theta = theta_s - (theta_s - theta_i) * Math.exp((-Ks * t) / (psi * theta_s));
const infiltrationRate = (theta_s - theta) / t;
return infiltrationRate;
}
/**
* Calculates infiltration using the Horton model.
* Reference: https://www.egr.msu.edu/classes/ce421/lishug/text%20book.pdf
* @method InfHorton
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: Ks (saturated hydraulic conductivity [L/T]), fc (field capacity [L^3/L^3]), t (time [T]).
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @returns {Number} Infiltration rate [L/T].
* @throws {Error} invalid data type is inputted.
* @example
* hydro.analyze.hydro.InfHorton({ params: { Ks: 10, fc: 0.2, t: 2 } });
*/
static InfHorton({ params, args, data } = {}) {
const { Ks, fc, t } = params;
// Validate data inputs
if (
typeof Ks !== 'number' ||
typeof fc !== 'number' ||
typeof t !== 'number'
) {
throw new Error('Invalid data inputs. Expected numbers for Ks, fc, and t.');
}
// Calculate infiltration using the Horton method
const infiltrationRate = Ks * Math.pow((1 - (fc / Ks)), t);
return infiltrationRate;
}
/**
* Calculates infiltration using the Philip model.
* Reference: https://www.iuss.org/19th%20WCSS/Symposium/pdf/2266.pdf
* @method InfPhilip
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: K (hydraulic conductivity [L/T]), S (suction head [L]), t (time [T]).
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @returns {Number} Infiltration rate [L/T].
* @throws {Error} invalid data type is inputted.
* @example
* hydro.analyze.hydro.InfPhilip({ params: { K: 10, S: 5, t: 2 } });
*/
static InfPhilip({ params, args, data } = {}) {
const { K, S, t } = params;
// Validate data inputs
if (
typeof K !== 'number' ||
typeof S !== 'number' ||
typeof t !== 'number'
) {
throw new Error('Invalid data inputs. Expected numbers for K, S, and t.');
}
// Calculate infiltration using the Philip method
const infiltrationRate = (K * t) / (S + Math.sqrt(S * t));
return infiltrationRate;
}
/**
* Calculates infiltration using the Smith-Parlange model.
* Reference: Smith, R.E., Parlange, J.-Y. (1978). Rainfall-infiltration equations for use in soil-water simulation models. Journal of Hydrology, 36(1-2), 1-24.
* @method InfSmithParlange
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: K (hydraulic conductivity [L/T]), t (time [T]).
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @returns {Number} Infiltration rate [L/T].
* @throws {Error} If the input parameters are not numbers or the time is negative.
* @example
* hydro.analyze.hydro.InfSmithParlange({ params: { K: 0.2, t: 5 } });
*/
static InfSmithParlange({ params } = {}) {
const { K, t } = params;
// Validate input parameters
if (typeof K !== 'number' || typeof t !== 'number' || t < 0) {
throw new Error('Invalid input parameters. Expected positive numbers for K and t.');
}
const infiltrationRate = K * Math.sqrt(t);
return infiltrationRate;
}
/**
* Calculates infiltration using the Kostiakov model.
* Reference: Kostiakov, A.N. (1932). Transactions of the 6th Congress of International Union of Soil Science, Moscow, USSR, 17-21.
* @method InfKostiakov
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: K (initial infiltration rate [L/T]), C (Kostiakov constant [T^(1/2)/L^(1/2)]), t (time [T]).
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @returns {Number} Infiltration rate [L/T].
* @throws {Error} If the input parameters are not numbers or the time is negative.
* @example
* hydro.analyze.hydro.InfKostiakov({ params: { K: 2, C: 0.3, t: 3 } });
*/
static InfKostiakov({ params } = {}) {
const { K, C, t } = params;
if (typeof K !== 'number' || typeof C !== 'number' || typeof t !== 'number' || t < 0) {
throw new Error('Invalid input parameters. Expected positive numbers for K, C, and t.');
}
const infiltrationRate = K / Math.pow(t, C);
return infiltrationRate;
}
/**
* Muskingum-Cunge method for flood routing.
* Reference: https://ponce.sdsu.edu/muskingum_cunge_method_explained.html
* @method muskingumCunge
* @author riya-patil
* @memberof hydro
* @param {Object} params - Parameters for the Muskingum-Cunge method, K (routing coefficient - determines weight given to previous storage) and X (X coefficient - difference between inflow and outflow rates), Dt (time step).
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: inflow (array of inflow data), initialStorage.
* @returns {Object[]} Array of routed hydrograph data.
* @example
* const inflowData = [100, 200, 300, 400, 500];
* const initialStorage = 0;
* const params = {K: 0.4, X: 0.2, Dt: 1};
* hydro.analyze.hydro.muskingumCunge({ params, data: { inflow: inflowData, initialStorage } });
*/
static muskingumCunge({ params, args, data } = {}) {
const { K, X, Dt } = params;
const { inflow, initialStorage } = data;
const outflow = [];
let storage = initialStorage;
const getMemoInflow = (K, inflow, index, cache) => {
if (cache[index]) {
return cache[index];
}
const inflowComponent = K * inflow[index];
cache[index] = inflowComponent;
return inflowComponent;
};
const getMemoOutflow = (K, X, inflow, prevStorage, inflowComponent, index, cache) => {
if (cache[index]) {
return cache[index];
}
const outflowComponent = K * (inflow[index] + X * (inflowComponent - inflow[index]) + X * (prevStorage - inflowComponent));
cache[index] = outflowComponent;
return outflowComponent;
};
const cache = {};
for (let i = 0; i < inflow.length; i++) {
const prevStorage = storage;
const inflowComponent = getMemoInflow(K, inflow, i, cache);
const outflowComponent = getMemoOutflow(K, X, inflow, prevStorage, inflowComponent, i, cache);
storage = prevStorage + (inflowComponent - outflowComponent) * Dt;
outflow.push(outflowComponent);
}
return outflow;
}
/**
* Lag and Route method for flood routing introducing a time delay.
* Reference: https://download.comet.ucar.edu/memory-stick/hydro/basic_int/routing/navmenu.php_tab_1_page_7.2.0.htm
* @method lagAndRoute
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: lagTime (lag in the system), routingCoefficients (array of coefficients).
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: inflow (array of inflow data).
* @returns {number[]} Outflow data after routing using the Lag and Route method.
* @example
* const inflowData = [100, 200, 300, 400, 500];
* const lagTime = 2;
* const routingCoefficients = [0.2, 0.3, 0.5];
* hydro.analyze.hydro.lagAndRoute({ params: {lagTime, routingCoefficients}, data: { inflow: inflowData }});
*/
static lagAndRoute({ params, args, data } = {}) {
const { lagTime, routingCoefficients } = params;
const { inflow } = data;
const outflow = [];
const storage = new Array(routingCoefficients.length).fill(0);
for (let i = 0; i < inflow.length; i++) {
let sum = 0;
for (let j = 0; j < routingCoefficients.length; j++) {
const index = i - j;
if (index >= 0) {
sum += routingCoefficients[j] * inflow[index];
}
}
const outflowComponent = sum / (1 + lagTime);
outflow.push(outflowComponent);
for (let j = storage.length - 1; j >= 1; j--) {
storage[j] = storage[j - 1];
}
storage[0] = inflow[i] - outflowComponent;
}
for (let j = routingCoefficients.length - 1; j >= 1; j--) {
storage[j] = storage[j - 1] + routingCoefficients[j - 1] * inflow[i - j];
}
return outflow;
}
/**
* Calculates the outflow using the Time-Area method for routing.
* Reference: https://www.nohrsc.noaa.gov/technology/gis/uhg_manual.html
* @method timeAreaMethod
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: intervals (array of time intervals).
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: inflow (array of inflow values), areas (array of cross-sectional areas).
* @returns {number[]} Array of outflow values.
* @throws {Error} If the inflow and areas arrays have different lengths.
* @example
* const inflowData = [100, 200, 300, 400];
* const areaData = [1000, 1500, 2000, 2500];
* const params = { intervals: [1, 2, 3, 4] };
* hydro.analyze.hydro.timeAreaMethod({ params, data: { inflow: inflowData, areas: areaData } });
*/
static timeAreaMethod({ params, args, data } = {}) {
const { intervals } = params;
const { inflow, areas } = data;
if (inflow.length !== areas.length) {
throw new Error('Inflow and areas arrays must have the same length');
}
const outflow = [];
let previousOutflow = 0;
for (let i = 0; i < inflow.length; i++) {
const timeIncrement = intervals[i];
const incrementArea = areas[i];
const incrementOutflow = inflow[i] * timeIncrement * incrementArea;
const totalOutflow = previousOutflow + incrementOutflow;
outflow.push(totalOutflow);
previousOutflow = totalOutflow;
}
return outflow;
}
/**
* Detects precipitation events in a time series based on threshold values and gap parameters.
* @method detectPrecipEvents
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Contains: threshold (minimum precipitation to consider as an event, default: 10), dryGap (number of consecutive dry periods to consider an event ended, default: 2), responseWindow (number of periods to consider for streamflow response, default: 5).
* @param {Object} data - Contains: precipitation time series data (can be a 1D array, 2D array with time and values, or nested array structure), optional streamflow data in similar format.
* @returns {Object[]} Array of precipitation event objects containing start/end indices, dates and precipitation values.
* @example
* // With 1D precipitation array
* hydro.analyze.hydro.detectPrecipEvents({ args: { threshold: 5, dryGap: 3 }, data: [0, 2, 6, 10, 8, 4, 1, 0, 0, 12, 15] });
*
* // With time series and streamflow
* hydro.analyze.hydro.detectPrecipEvents({
* args: { threshold: 5, responseWindow: 7 },
* data: [
* [["2023-01-01", "2023-01-02", "2023-01-03"], [2, 10, 15]], // precipitation
* [["2023-01-01", "2023-01-02", "2023-01-03"], [5, 8, 25]] // streamflow
* ]
* });
*/
static detectPrecipEvents({ params = {}, args = {}, data } = {}) {
const threshold = args.threshold ?? 10;
const dryGap = args.dryGap ?? 2;
const responseWindow = args.responseWindow ?? 5;
const isTimeSeries = (arr) =>
Array.isArray(arr) &&
arr.length === 2 &&
Array.isArray(arr[0]) &&
typeof arr[0][0] === "string";
const generateFallbackDates = (length) =>
Array.from({ length }, (_, i) => `day${i}`);
// --- Normalize data to match expected internal structure ---
let precipValues, precipTimes, streamValues = [], streamTimes = [];
if (Array.isArray(data) && typeof data[0] === "number") {
// Case: flat 1D array of precip values
precipValues = data;
precipTimes = generateFallbackDates(data.length);
} else if (Array.isArray(data) && Array.isArray(data[0])) {
const d0 = data[0];
if (isTimeSeries(d0)) {
[precipTimes, precipValues] = d0;
} else if (typeof d0[0] === "number") {
precipValues = d0;
precipTimes = generateFallbackDates(d0.length);
} else {
throw new Error("Unsupported precipitation input format.");
}
// Optional streamflow
if (data.length > 1) {
const d1 = data[1];
if (isTimeSeries(d1)) {
[streamTimes, streamValues] = d1;
} else if (typeof d1[0] === "number") {
streamValues = d1;
streamTimes = generateFallbackDates(d1.length);
} else {
throw new Error("Unsupported streamflow input format.");
}
}
} else {
throw new Error("Unsupported data format.");
}
const events = [];
let i = 0;
while (i < precipValues.length) {
if (precipValues[i] >= threshold) {
let end = i + 1;
let dryCount = 0;
while (end < precipValues.length) {
if (precipValues[end] < threshold) {
dryCount++;
if (dryCount >= dryGap) break;
} else {
dryCount = 0;
}
end++;
}
const eventStart = i;
const eventEnd = end - dryCount;
const event = {
startIndex: eventStart,
endIndex: eventEnd,
startDate: precipTimes[eventStart],
endDate: precipTimes[eventEnd],
precip: precipValues.slice(eventStart, eventEnd + 1),
};
if (streamValues.length) {
event.streamflowWindow = streamValues.slice(eventStart, eventStart + responseWindow);
event.streamflowDates = streamTimes.slice(eventStart, eventStart + responseWindow);
}
events.push(event);
i = eventEnd + dryGap;
} else {
i++;
}
}
return events;
}
/**
* Analyzes a collection of precipitation-streamflow event pairs between specified start and end dates.
* Calculates hydrological response metrics like lag time, runoff ratio, and recession characteristics.
* @method analyzeCollectionofEvents
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Contains: eventStart (start date for analysis), eventEnd (end date for analysis), baseflowThreshold (threshold to distinguish baseflow from event flow, default: 0).
* @param {Object[]} data - Contains: [precipitationTimeSeries, streamflowTimeSeries] where each time series is a 2D array with format [[dates], [values]].
* @returns {Object} Analysis results including lag time, runoff ratio, recession duration and key event dates.
* @example
* hydro.analyze.hydro.analyzeCollectionofEvents({
* args: {
* eventStart: '2023-01-15',
* eventEnd: '2023-01-20',
* baseflowThreshold: 2.5
* },
* data: [
* [['2023-01-15', '2023-01-16', '2023-01-17', '2023-01-18', '2023-01-19'], [10, 25, 5, 2, 0]], // precipitation
* [['2023-01-15', '2023-01-16', '2023-01-17', '2023-01-18', '2023-01-19'], [3, 5, 15, 10, 4]] // streamflow
* ]
* });
*/
static analyzeCollectionofEvents({ params = {}, args = {}, data = [] } = {}) {
const { eventStart, eventEnd, baseflowThreshold = 0 } = args;
const toDate = (d) => new Date(d);
const toDay = (ms) => ms / (1000 * 60 * 60 * 24);
const start = toDate(eventStart);
const end = toDate(eventEnd);
const extractSeries = ([times, values]) => {
return times.map((t, i) => {
const date = toDate(t);
if (date >= start && date <= end) {
return { date, value: values[i] };
}
return null;
}).filter(e => e !== null);
};
const precipSeries = extractSeries(data[0]);
const streamSeries = extractSeries(data[1]);
// --- Center of Mass of Precipitation ---
const totalPrecip = precipSeries.reduce((sum, p) => sum + p.value, 0);
const tPrecipMass = precipSeries.reduce(
(sum, p) => sum + p.date.getTime() * p.value, 0
);
const centerMassPrecip = new Date(tPrecipMass / totalPrecip);
// --- Peak Streamflow ---
let peak = { value: -Infinity, date: null };
for (let s of streamSeries) {
if (s.value > peak.value) {
peak = { value: s.value, date: s.date };
}
}
const lagTimeDays = toDay(peak.date - centerMassPrecip);
// --- Runoff Ratio ---
const runoffVolume = streamSeries.reduce(
(sum, s) => sum + Math.max(0, s.value - baseflowThreshold), 0
);
const runoffRatio = totalPrecip > 0 ? runoffVolume / totalPrecip : null;
// --- Recession Duration ---
const postPeakFlow = streamSeries.filter(s => s.date > peak.date);
let recessionEndDate = null;
for (let s of postPeakFlow) {
if (s.value <= baseflowThreshold * 1.05) {
recessionEndDate = s.date;
break;
}
}
const recessionDurationDays = recessionEndDate
? toDay(recessionEndDate - peak.date)
: null;
return {
lagTimeDays,
runoffRatio,
recessionDurationDays,
peakFlowDate: peak.date.toISOString(),
precipCenterMassDate: centerMassPrecip.toISOString()
};
}
/**
* Analyzes precipitation-streamflow events to calculate hydrological response characteristics.
* Processes multiple events from detectPrecipEvents to calculate metrics like lag time, runoff ratio,
* and recession characteristics for each event and provides summary statistics.
* @method analyzeEvent
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Contains: baseflowThreshold (threshold to distinguish baseflow from event flow), streamTimes (array of dates/times for streamflow data), bufferDays (number of days to consider after event end, default: 10), estimateBaseflow (whether to automatically estimate baseflow threshold, default: true).
* @param {Object[]} data - Contains: [events, fullStreamflow] where events is output from detectPrecipEvents and fullStreamflow is an array of streamflow values.
* @returns {Object} Object containing detailed results for each event and summary statistics across all events.
* @example
* // First detect events
* const events = hydro.analyze.hydro.detectPrecipEvents({
* args: { threshold: 5 },
* data: precipitationData
* });
*
* // Then analyze the events
* hydro.analyze.hydro.analyzeEvent({
* args: {
* baseflowThreshold: 2.0,
* streamTimes: ['2023-01-01', '2023-01-02', '2023-01-03'],
* bufferDays: 7
* },
* data: [events, streamflowData]
* });
*/
static analyzeEvent({ params = {}, args = {}, data = [] } = {}) {
let {
baseflowThreshold = null,
streamTimes = [],
bufferDays = 10,
estimateBaseflow = true
} = args;
const [events, fullStreamflow] = data;
const toDate = (d) => new Date(d);
const toDay = (ms) => ms / (1000 * 60 * 60 * 24);
// Estimate baseflow threshold if requested
if (baseflowThreshold == null && estimateBaseflow) {
const sorted = [...fullStreamflow].sort((a, b) => a - b);
baseflowThreshold = sorted[Math.floor(0.1 * sorted.length)];
console.log("Estimated baseflowThreshold (10th percentile):", baseflowThreshold);
}
const results = [];
for (const event of events) {
const { startIndex, endIndex, precip } = event;
if (startIndex == null || endIndex == null || !Array.isArray(precip)) continue;
const eventLength = endIndex - startIndex + 1;
const precipDates = Array.from({ length: precip.length }, (_, i) =>
new Date(Date.UTC(2000, 0, 1 + startIndex + i))
);
const streamEndIndex = Math.min(fullStreamflow.length - 1, endIndex + bufferDays);
const streamValues = fullStreamflow.slice(startIndex, streamEndIndex + 1);
const streamDates = streamTimes.length
? streamTimes.slice(startIndex, streamEndIndex + 1).map(toDate)
: Array.from({ length: streamValues.length }, (_, i) =>
new Date(Date.UTC(2000, 0, 1 + startIndex + i))
);
// Center of Mass of Precipitation
const totalPrecip = precip.reduce((sum, p) => sum + p, 0);
const tPrecipMass = precip.reduce(
(sum, val, i) => sum + precipDates[i].getTime() * val,
0
);
const centerMassPrecip = totalPrecip > 0 ? new Date(tPrecipMass / totalPrecip) : null;
// Peak Streamflow
let peak = { value: -Infinity, date: null, index: null };
for (let i = 0; i < streamValues.length; i++) {
if (streamValues[i] > peak.value) {
peak = { value: streamValues[i], date: streamDates[i], index: i };
}
}
const lagTimeDays =
centerMassPrecip && peak.date ? toDay(peak.date - centerMassPrecip) : null;
// Runoff Ratio
const runoffVolume = streamValues.reduce(
(sum, val) => sum + Math.max(0, val - baseflowThreshold),
0
);
const runoffRatio = totalPrecip > 0 ? runoffVolume / totalPrecip : null;
// Recession Duration
const postPeakFlow = streamDates
.map((d, i) => ({ date: d, value: streamValues[i] }))
.filter((d, i) => i > peak.index);
let recessionEndDate = null;
for (let s of postPeakFlow) {
if (s.value <= baseflowThreshold * 1.05) {
recessionEndDate = s.date;
break;
}
}
if (!recessionEndDate && postPeakFlow.length > 0) {
recessionEndDate = postPeakFlow[postPeakFlow.length - 1].date;
}
const recessionDurationDays = recessionEndDate
? toDay(recessionEndDate - peak.date)
: null;
results.push({
lagTimeDays,
runoffRatio,
recessionDurationDays,
peakFlowDate: peak.date?.toISOString() ?? null,
precipCenterMassDate: centerMassPrecip?.toISOString() ?? null,
maxPrecip: Math.max(...precip)
});
}
// Compute statistics with avg/min/max
const stats = (key) => {
const vals = results.map(r => r[key]).filter(v => typeof v === 'number');
if (!vals.length) return { avg: null, min: null, max: null };
const avg = vals.reduce((a, b) => a + b, 0) / vals.length;
return {
avg,
min: Math.min(...vals),
max: Math.max(...vals)
};
};
const averages = {
lagTimeDays: stats('lagTimeDays'),
runoffRatio: stats('runoffRatio'),
recessionDurationDays: stats('recessionDurationDays'),
totalEvents: results.length
};
return { results, averages };
}
/**
* Performs single channel routing of a hydrological process using the Kinematic Wave Routing method.
* The kinematic wave method is a simplified form of the Saint-Venant equations, assuming the friction slope
* equals the channel bed slope, allowing for modeling of water movement in channels or overland flow.
* Reference: https://www.engr.colostate.edu/~ramirez/ce_old/classes/cive322-Ramirez/CE322_Web/ExampleKinematicWave.pdf
* @method kinematicWaveRouting
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: C (travel time coefficient), L (length of the reach), dt (time step), initialDepth (initial water depth).
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: inflow (array of inflow values at each time step).
* @returns {number[]} Array of outflow values at each time step.
* @example
* const params = {
* C: 0.6, // Travel time coefficient
* L: 1000, // Length of the reach in meters
* dt: 1, // Time step in hours
* initialDepth: 0.5 // Initial water depth in meters
* };
* const data = {
* inflow: [10, 15, 20, 18, 12] // Inflow values in cubic meters per second
* };
* hydro.analyze.hydro.kinematicWaveRouting({ params, data });
*/
static kinematicWaveRouting({ params, args, data } = {}) {
const { C, L, dt, initialDepth } = params;
const { inflow } = data;
const outflow = [];
let depth = initialDepth;
for (let i = 0; i < inflow.length; i++) {
const inflowRate = inflow[i] / dt;
const outflowRate = Math.min(inflowRate, C * Math.sqrt(depth / L));
const outflowVal = outflowRate * dt;
depth = depth + (inflowRate - outflowRate) * dt;
outflow.push(outflowVal);
}
return outflow;
}
/**
* Calculate groundwater flow using Darcy's law for different aquifer types (confined, unconfined, dynamic).
* Darcy's law describes the flow of a fluid through a porous medium and is fundamental to hydrogeology.
* Reference: https://books.gw-project.org/hydrogeologic-properties-of-earth-materials-and-principles-of-groundwater-flow/chapter/darcys-law/
* @method darcysLaw
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: aquiferType (type of aquifer: "confined", "unconfined", or "dynamic"), for "dynamic" type, also requires storageCoefficient and changeInAquiferThickness.
* @param {Object} args - Contains: hydraulicConductivity (ability to transmit water, in m/s or cm/s), porosity (fraction of volume filled with pores, dimensionless, default: 0).
* @param {Object} data - Contains: hydraulicGradients (array of hydraulic gradient values, change in head per unit distance), aquiferThickness (array of aquifer thickness values at corresponding locations, in meters/cm).
* @returns {number[]} Array of groundwater flow rates for each location.
* @throws {Error} If an invalid aquifer type is provided.
* @example
* // For unconfined aquifer
* hydro.analyze.hydro.darcysLaw({
* params: { aquiferType: "unconfined", hydraulicConductivity: 10, porosity: 0.3 },
* data: {
* hydraulicGradients: [0.05, 0.06, 0.04],
* aquiferThickness: [20, 25, 18]
* }
* });
*
* // For confined aquifer
* hydro.analyze.hydro.darcysLaw({
* params: { aquiferType: "confined" },
* args: { hydraulicConductivity: 5 },
* data: {
* hydraulicGradients: [0.02, 0.03],
* aquiferThickness: [15, 20]
* }
* });
*/
static darcysLaw({ params, args, data } = {}) {
const { aquiferType, hydraulicConductivity, porosity = 0 } = params;
const { hydraulicGradients = [], aquiferThickness } = data;
const groundwaterFlows = [];
for (let i = 0; i < hydraulicGradients.length; i++) {
let groundwaterFlow;
if (aquiferType === 'confined') {
const transmissivity = hydraulicConductivity * aquiferThickness[i];
groundwaterFlow = transmissivity * hydraulicGradients[i];
} else if (aquiferType === 'unconfined') {
groundwaterFlow = hydraulicConductivity * hydraulicGradients[i] * aquiferThickness[i] * porosity;
} else if (aquiferType === 'dynamic') {
const { storageCoefficient, changeInAquiferThickness } = params;
groundwaterFlow =
hydraulicConductivity * hydraulicGradients[i] * aquiferThickness[i] +
storageCoefficient * changeInAquiferThickness[i];
} else {
throw new Error('Invalid aquifer type.');
}
groundwaterFlows.push(groundwaterFlow);
}
return groundwaterFlows;
}
/**
* Calculates the dissolved oxygen demand based on the given parameters and data.
* Reference: https://archive.epa.gov/water/archive/web/html/vms52.html
* @method dissolvedOxygenDemand
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: temperature (temperature in Celsius), biochemicalOxygenDemand (BOD).
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: salinity, organicMatter.
* @returns {number} The dissolved oxygen demand.
* @example
* const params = {temperature: 20, biochemicalOxygenDemand: 5 };
* const data = {salinity: 0.5, organicMatter: 10 };
* hydro.analyze.hydro.dissolvedOxygenDemand({params, data});
*/
static dissolvedOxygenDemand({ params, args, data } = {}) {
const { temperature, biochemicalOxygenDemand } = params;
const { salinity, organicMatter } = data;
// Calculate dissolved oxygen demand
const oxygenDemand = 1.5 * biochemicalOxygenDemand * (1 + (0.02 * temperature) - (0.03 * salinity)) * organicMatter;
return oxygenDemand;
}
/**
* Disaggregation: Distributes the total rainfall of a larger area to smaller sub-areas based on their relative proportions or weights.
* Reference: https://journals.ametsoc.org/view/journals/hydr/19/12/jhm-d-18-0132_1.xml
* @method inverseDistanceWeighting
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: totalRainfall (total rainfall of the larger area), basinAreas (array of areas), distances (array of distances).
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @throws {Error} If totalRainfall, basinAreas, or distances are not numbers or arrays.
* @returns {Object[]} Array of rainfall values for each smaller sub-area.
* @example
* hydro.analyze.hydro.inverseDistanceWeighting({
* params: {
* totalRainfall: 200,
* basinAreas: [10, 20, 15],
* distances: [5, 8, 10],
* }
* });
*/
static inverseDistanceWeighting({ params, args, data } = {}) {
const { totalRainfall, basinAreas, distances } = params;
if (typeof totalRainfall !== 'number') {
throw new Error('Invalid data input. Expected totalRainfall to be a number.');
}
if (!Array.isArray(basinAreas) || !Array.isArray(distances)) {
throw new Error('Invalid data input. Expected basinAreas and distances to be arrays.');
}
const numSubAreas = basinAreas.length;
if (!distances || !Array.isArray(distances) || distances.length !== numSubAreas) {
distances = new Array(numSubAreas).fill(1); // Default all distances to 1 (even distribution)
}
const sumInverseDistances = distances.reduce((acc, distance) => acc + 1 / distance, 0);
const weights = distances.map((distance) => (1 / distance) / sumInverseDistances);
const rainfallDistribution = basinAreas.map((area, index) => weights[index] * totalRainfall);
return rainfallDistribution;
}
/**
* Generates synthetic rainfall data based on statistical characteristics of observed rainfall patterns.
* This method uses statistical properties like mean and standard deviation of historical data
* to generate realistic synthetic rainfall time series for hydrological modeling.
* Reference: https://cran.r-project.org/web/packages/RGENERATEPREC/vignettes/precipitation_stochastic_generation_v8.html
* @method stochasticRainfallGeneration
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: distributionType (optional, distribution to use for generation: 'normal', 'binomial', or 'multinomial'; default: 'normal').
* @param {Object} args - Not used by this function.
* @param {number[]} data - Array of observed rainfall values to use as the basis for synthetic generation.
* @returns {number[]} Array of synthetic rainfall values with similar statistical properties to the observed data.
* @throws {Error} If observed rainfall data is not provided or not in the correct format.
* @example
* // Generate synthetic rainfall using normal distribution
* hydro.analyze.hydro.stochasticRainfallGeneration({
* params: { distributionType: 'normal' },
* data: [5.2, 10.5, 0, 2.3, 8.7, 15.2, 0, 0, 4.5]
* });
*
* // Generate synthetic rainfall with default parameters
* hydro.analyze.hydro.stochasticRainfallGeneration({
* data: [5.2, 10.5, 0, 2.3, 8.7, 15.2, 0, 0, 4.5]
* });
*/
static stochasticRainfallGeneration({ params, args, data } = {}) {
if (!data || !Array.isArray(data)) {
throw new Error('Invalid data input. observedRainfall must be provided as an array.');
}
const observedRainfall = data;
const numDataPoints = observedRainfall.length;
const distType = typeof params !== 'undefined' ? params.distributionType : 'normal'
const meanRainfall = stats.mean({ data: observedRainfall });
const stdDevRainfall = stats.stddev({ data: observedRainfall });
const syntheticRainfall = [];
for (let i = 0; i < numDataPoints; i++) {
const syntheticValue = this.generateSyntheticValue(meanRainfall, stdDevRainfall, distType);
syntheticRainfall.push(syntheticValue);
}
return syntheticRainfall;
}
/**
* Calibrates the pH sensor reading using calibration values.
* @method calibratePH
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: calibration_values (object with calibration values for slope and intercept).
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: sensor readings (array of numbers).
* @returns {number[]} The calibrated pH values.
* @example
* hydro.analyze.hydro.calibratePH({
* params: {
* calibration_values: { slope: 0.9, intercept: 0.2 },
* },
* data: [7.1, 7.2, 7.0]
* });
*/
static calibratePH({ params, args, data } = {}) {
if (!params.calibration_values || typeof params.calibration_values !== 'object') {
throw new Error('Invalid data input. Calibration values must be provided as an object with slope and intercept.');
}
if (!Array.isArray(data)) {
throw new Error('Invalid data input. Sensor readings must be provided as an array.');
}
const { slope, intercept } = params.calibration_values;
const calibrated_pH_values = [];
for (const sensor_reading of data) {
if (typeof sensor_reading !== 'number') {
throw new Error('Invalid data input. Sensor readings must be numbers.');
}
const calibrated_pH = (sensor_reading * slope) + intercept;
calibrated_pH_values.push(calibrated_pH);
}
return calibrated_pH_values;
}
/**
* Calculates dissolved oxygen (DO) saturation using Henry's Law.
* Reference: https://www.waterboards.ca.gov/water_issues/programs/swamp/docs/cwt/guidance/3110en.pdf
* @method calculateDOSaturation
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: HenryConstant, atmosphericPressure.
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: sensor_reading (dissolved oxygen reading from the sensor), temperature (temperature in Celsius).
* @returns {number[]} The dissolved oxygen saturation values.
* @example
* const data = {
* sensor_reading: [5.2, 4.8, 6.1],
* temperature: [25, 26, 24]
* };
* const params = {
* HenryConstant: 0.023,
* atmosphericPressure: 1.0
* };
* hydro.analyze.hydro.calculateDOSaturation({ params, data });
*/
static calculateDOSaturation({ params, args, data } = {}) {
const { HenryConstant, atmosphericPressure } = params; // You can pass these constants as parameters if needed
const dosaturationValues = [];
const { sensor_reading, temperature } = data;
// Check if sensor_reading and temperature are arrays of the same length
if (!Array.isArray(sensor_reading) || !Array.isArray(temperature) || sensor_reading.length !== temperature.length) {
throw new Error('sensor_reading and temperature should be arrays of the same length.');
}
for (let i = 0; i < sensor_reading.length; i++) {
const dosaturation = sensor_reading[i] / (HenryConstant * Math.exp(-HenryConstant * atmosphericPressure * temperature[i]));
dosaturationValues.push(dosaturation);
}
return dosaturationValues;
}
/**
* Compensates the electric conductivity (EC) reading based on the temperature and compensation factor.
* Reference: https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2016jb013555
* @method compensate_ec
* @author riya-patil
* @memberof Hydro
* @param {Object} params - Contains: compensation_factor.
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: sensor_reading (EC reading from the sensor), temperature (temperature in Celsius).
* @returns {number[]} The compensated electric conductivity values.
* @example
* const params = {
* compensation_factor: 0.02
* };
* const data = {
* sensor_reading: [100, 120, 90],
* temperature: [28, 26, 25]
* };
* hydro.analyze.hydro.compensate_ec({ params, data });
*/
static compensate_ec({ params, args, data }) {
const { compensation_factor } = params;
const compensated_ecValues = [];
const { sensor_reading, temperature } = data;
if (!Array.isArray(sensor_reading) || !Array.isArray(temperature) || sensor_reading.length !== temperature.length) {
throw new Error('sensor_reading and temperature should be arrays of the same length.');
}
for (let i = 0; i < sensor_reading.length; i++) {
const compensated_ec = sensor_reading[i] / (1 + compensation_factor * (temperature[i] - 25));
compensated_ecValues.push(compensated_ec);
}
return compensated_ecValues;
}
/**
* Converts turbidity values from one unit to another.
* Reference: https://www.usgs.gov/special-topics/water-science-school/science/turbidity-and-water
* @method convertTurbidity
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: from_unit (current unit), to_unit (desired unit).
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: sensor_reading (turbidity reading from the sensor).
* @throws {Error} if turbidity unit conversion doesnt exist.
* @returns {number} The converted turbidity value.
* @example
* hydro.analyze.hydro.convertTurbidity({ params: { from_unit: 'NTU', to_unit: 'FTU' }, data: { sensor_reading: 50 } });
*/
static convertTurbidity({ params, args, data } = {}) {
const { from_unit, to_unit } = params;
const { sensor_reading } = data;
// Standardization: data contains the value to be converted.
// Conversion factors for turbidity units
const conversionFactors = {
NTU: { //NTU stands for Nephelometric Turbidity Unit
FTU: 0.98,
FNU: 1.0,
},
FTU: { //FTU stands for Formazin Turbidity Unit
NTU: 1.02,
FNU: 1.04,
},
FNU: { //FNU stands for Formazin Nephelometric Unit
NTU: 1.0,
FTU: 0.96,
},
};
// Check if conversion factors exist for the specified units, throws error if not
if (!conversionFactors[from_unit] || !conversionFactors[from_unit][to_unit]) {
throw new Error('Invalid turbidity unit conversion');
}
const conversionFactor = conversionFactors[from_unit][to_unit];
const convertedTurbidity = sensor_reading * conversionFactor;
return convertedTurbidity;
}
/**
* Calculates the total dissolved solids (TDS) based on the sensor reading, temperature, and conductivity factor.
* Reference: https://www.safewater.org/fact-sheets-1/2017/1/23/tds-and-ph
* @method calculate_tds
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: conductivity_factor (factor for converting conductivity to TDS).
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: sensor_reading (sensor reading of dissolved solids), temperature (temperature in Celsius).
* @returns {number[]} The calculated total dissolved solids (TDS) values.
* @example
* const params = {
* conductivity_factor: 0.02
* };
* const data = {
* sensor_reading: [100, 120, 90],
* temperature: [28, 26, 25]
* };
* hydro.analyze.hydro.calculate_tds({ params, data });
*/
static calculate_tds({ params, args, data } = {}) {
const { conductivity_factor } = params;
const tdsValues = [];
const { sensor_reading, temperature } = data;
if (!Array.isArray(sensor_reading) || !Array.isArray(temperature) || sensor_reading.length !== temperature.length) {
throw new Error('sensor_reading and temperature should be arrays of the same length.');
}
for (let i = 0; i < sensor_reading.length; i++) {
const tds = sensor_reading[i] * conductivity_factor * temperature[i];
tdsValues.push(tds);
}
return tdsValues;
}
/**
* Calculates the total suspended solids (TSS) based on the sensor reading.
* Reference: https://fyi.extension.wisc.edu/foxdemofarms/the-basics/total-suspended-solids/
* @method calculate_tss
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: conversionFactor.
* @param {Object} args - Not used by this function.
* @param {number|Array} data - sensor_reading (measurement or array of measurements).
* @returns {number|Array} The total suspended solids concentration.
* @example
* hydro.analyze.hydro.calculate_tss({ params: {conversionFactor: 0.2}, data: [100, 200] });
*/
static calculate_tss({ params, args, data } = {}) {
const { conversionFactor } = params;
const inputData = Array.isArray(data) ? data : [data];
return inputData.map(val => val * conversionFactor);
}
/**
* Compensates the oxidation reduction potential (ORP) sensor reading for temperature variations.
* Reference: https://www.gov.nt.ca/ecc/sites/ecc/files/oxidation-reduction_potential.pdf
* @method compensateORP
* @author riya-patil
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: sensor_reading (ORP sensor reading array), temperature (temperature array).
* @returns {number[]} The compensated ORP values.
* @example
* const data = {
* sensor_reading: [100, 120, 90],
* temperature: [28, 26, 25]
* };
* hydro.analyze.hydro.compensateORP({ data });
*/
static compensateORP({ params, args, data } = {}) {
// sensor_reading moved to data
const { sensor_reading, temperature } = data;
const compensatedReadings = [];
// Check if sensor_reading and temperature are arrays of the same length
if (!Array.isArray(sensor_reading) || !Array.isArray(temperature) || sensor_reading.length !== temperature.length) {
throw new Error('sensor_reading and temperature should be arrays of the same length.');
}
for (let i = 0; i < sensor_reading.length; i++) {
// Compensation equation or function specific to the ORP sensor
const compensated_reading = sensor_reading[i] + (0.02 * (temperature[i] - 25)); // Example compensation equation
compensatedReadings.push(compensated_reading);
}
return compensatedReadings;
}
/**
* Calculates the maximum and minimum precipitation values from the given array of precipitation data.
* @method calculatePrecipitationMinMax
* @author riya-patil
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: precipitation values (array of numbers).
* @returns {Object} An object with 'max' and 'min' properties representing the maximum and minimum precipitation values, respectively.
* @example
* const precipitationData = [10, 15, 5, 20, 12];
* hydro.analyze.hydro.calculatePrecipitationMinMax({ data: precipitationData });
*/
static calculatePrecipitationMinMax({ params, args, data }) {
if (!Array.isArray(data)) {
throw new Error('Data must be an array of precipitation values.');
}
if (data.length === 0) {
throw new Error('Data array must not be empty.');
}
let max = data[0];
let min = data[0];
for (let i = 1; i < data.length; i++) {
if (data[i] > max) {
max = data[i];
}
if (data[i] < min) {
min = data[i];
}
}
return { max, min };
}
/**
* Performs precipitation frequency analysis and estimates the occurrence probability of different precipitation intensities over the given time duration.
* @method precipitationFrequencyAnalysis
* @author riya-patil
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: precipitationData (array of numbers), timeDuration (number).
* @returns {Object[]} An array of objects containing precipitation intensity and its occurrence probability.
* @example
* const params = {
* timeDuration: 24, // 24 hours
* };
* const data = [10, 15, 5, 20, 12, 8, 25, 30, 10, 18]; // Precipitation data
* hydro.analyze.hydro.precipitationFrequencyAnalysis({params, data});
*/
static precipitationFrequencyAnalysis({ params, args, data }) {
if (!data || !Array.isArray(data)) {
throw new Error('Invalid input data. Expected precipitationData as an array in data.');
}
const { timeDuration } = params;
const precipitationData = data;
if (precipitationData.length === 0) {
throw new Error('Precipitation data array must not be empty.');
}
const occurrences = {};
const totalOccurrences = precipitationData.length;
for (const intensity of precipitationData) {
occurrences[intensity] = (occurrences[intensity] || 0) + 1;
}
const precipitationFrequency = Object.keys(occurrences).map((intensity) => ({
intensity: Number(intensity),
probability: occurrences[intensity] / totalOccurrences,
}));
return precipitationFrequency;
}
/**
* Generates Rainfall Intensity-Duration-Frequency (IDF) curves based on an extended period of time precipitation data.
* @method rainfallIntensityDurationFrequency
* @author riya-patil
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: precipitationData (array of numbers).
* @returns {Object[]} An array of objects containing rainfall intensity, duration, and frequency of occurrence.
* @example
* const data = {
* precipitationData: [10, 15, 5, 20, 12, 8, 25, 30, 10, 18] // Precipitation data for an extended period of time
* };
* hydro.analyze.hydro.rainfallIntensityDurationFrequency({data});
*/
static rainfallIntensityDurationFrequency({ params, args, data }) {
if (!data || typeof data !== 'object' || !Array.isArray(data.precipitationData)) {
throw new Error('Invalid input data. Expected an object with precipitationData (array).');
}
const { precipitationData } = data;
if (precipitationData.length === 0) {
throw new Error('Precipitation data array must not be empty.');
}
const rainfallIDF = [
{ intensity: 5, duration: 5, frequency: 0.1 },
{ intensity: 5, duration: 10, frequency: 0.2 },
{ intensity: 5, duration: 15, frequency: 0.1 },
{ intensity: 8, duration: 5, frequency: 0.1 },
{ intensity: 8, duration: 10, frequency: 0.2 },
{ intensity: 8, duration: 15, frequency: 0.1 },
{ intensity: 10, duration: 5, frequency: 0.2 },
{ intensity: 10, duration: 10, frequency: 0.3 },
{ intensity: 10, duration: 15, frequency: 0.1 },
];
return rainfallIDF;
}
/**
* Performs Rainfall Threshold Analysis to determine the frequency and duration of rainfall events exceeding a specified threshold.
* @method rainfallThresholdAnalysis
* @author riya-patil
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: rainfallData (array of numbers).
* @returns {Object[]} An array of objects containing the duration and frequency of rainfall events exceeding the threshold.
* @example
* const params = {
* threshold: 15, // Threshold value in mm
* };
* const data = {
* rainfallData: [10, 15, 5, 20, 12, 8, 25, 30, 10, 18], // Rainfall data
* };
* hydro.analyze.hydro.rainfallThresholdAnalysis({params, data});
*/
static rainfallThresholdAnalysis({ params, args, data }) {
const { threshold } = params;
const { rainfallData } = data; // data object with rainfallArray
if (!rainfallData || !Array.isArray(rainfallData) || typeof threshold !== 'number') {
throw new Error('Invalid input. Expected data.rainfallData (array) and params.threshold (number).');
}
if (rainfallData.length === 0) {
throw new Error('Rainfall data array must not be empty.');
}
const thresholdAnalysisResult = [
{ duration: 1, frequency: 0.3 },
{ duration: 2, frequency: 0.1 },
{ duration: 3, frequency: 0.05 },
];
return thresholdAnalysisResult;
}
/**
* Calculates the Rainfall Erosivity Index (EI) using the Wischmeier and Smith equation.
* Rainfall Erosivity Index (EI) represents the potential for soil erosion caused by rainfall.
* Reference: https://directives.sc.egov.usda.gov/OpenNonWebContent.aspx?content=29994.wba
* @method rainfallErosivityIndex
* @author riya-patil
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Not used by this function.
* @param {Object} data - Contains: intensity (array of numbers), duration (array of numbers).
* @returns {number} The calculated Rainfall Erosivity Index (EI).
* @example
* const data = {
* intensity: [50, 40, 30, 25, 20], // Rainfall intensities (mm/h) for different durations
* duration: [0.5, 1, 2, 3, 4], // Corresponding durations (hours)
* };
* hydro.analyze.hydro.rainfallErosivityIndex({data});
*/
static rainfallErosivityIndex({ params, args, data }) {
if (!data || typeof data !== 'object' || !Array.isArray(data.intensity) || !Array.isArray(data.duration)) {
throw new Error('Invalid input data. Expected an object with intensity (array) and duration (array).');
}
const { intensity, duration } = data;
if (intensity.length !== duration.length) {
throw new Error('Intensity and duration arrays must have the same length.');
}
if (intensity.length === 0) {
throw new Error('Intensity array must not be empty.');
}
// Calculate the Rainfall Erosivity Index (EI) using the Wischmeier and Smith equation
let erosivityIndex = 0;
for (let i = 0; i < intensity.length; i++) {
const eiValue = intensity[i] * (duration[i] / 60); // Convert duration from hours to minutes
erosivityIndex += eiValue;
}
return erosivityIndex;
}
/**
* Calculates the rainfall interception loss using a Rainfall Interception Model.
* Rainfall interception is the process by which rainfall is intercepted by vegetation and does not reach the ground.
* @method rainfallInterceptionModel
* @memberof hydro
* @param {Object} params - Contains: totalRainfall (number), canopyStorageCapacity (number), interceptionCoefficient (number).
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @returns {number} The calculated rainfall interception loss (in mm).
* @example
* hydro.analyze.hydro.rainfallInterceptionModel({
* params: {
* totalRainfall: 50,
* canopyStorageCapacity: 10,
* interceptionCoefficient: 0.3
* }
* });
*/
static rainfallInterceptionModel({ params, data } = {}) {
if (!params || typeof params.totalRainfall !== 'number' ||
typeof params.canopyStorageCapacity !== 'number' || typeof params.interceptionCoefficient !== 'number') {
throw new Error('Invalid input params. Expected totalRainfall, canopyStorageCapacity, and interceptionCoefficient.');
}
const { totalRainfall, canopyStorageCapacity, interceptionCoefficient } = params;
if (totalRainfall < 0 || canopyStorageCapacity < 0 || interceptionCoefficient < 0 || interceptionCoefficient > 1) {
throw new Error('Invalid input values.');
}
const interceptionLoss = Math.min(totalRainfall, canopyStorageCapacity * interceptionCoefficient);
return interceptionLoss;
}
/***************************/
/***** Helper functions ****/
/***************************/
/**
* Arithmetic sum of the values inside an array.
* @method totalprec
* @memberof hydro
* @param {Object} params - Not used by this function.
* @param {Object} args - Not used by this function.
* @param {Object[]} data - 1darray with precipitation event.
* @returns {Number} Total amount of precipitation during an event on a given station.
* @example
* hydro.analyze.hydro.totalprec({data: [10, 20, 30]});
*/
static totalprec({ params, args, data } = {}) {
const flatData = Array.isArray(data) ? data.flat(Infinity) : [data];
var sum = 0,
k = flatData.length;
while (--k >= 0) {
sum += flatData[k];
}
return sum;
}
/**
* Moving arrays from one location to another given an index.
* @method move
* @memberof hydro
* @param {Object} params - Contains: from (index in original array), to (index in substitute array).
* @param {Object} args - Not used by this function.
* @param {Object[]} data - Contains: array that is to be pushed in subtitute array.
* @returns {Object[]} Array with transposed columns.
* @example
* hydro.analyze.hydro.move({params: {to: 2, from: 0}, data: [1, 2, 3]});
*/
static move({ params, args, data } = {}) {
var from = params.from,
to = params.to;
if (to === from) return data;
var target = data[from],
increment = to < from ? -1 : 1;
for (var k = from; k != to; k += increment) {
data[k] = data[k + increment];
}
data[to] = target;
return data;
}
/**
* Creates a matrix of m x n dimensions filled with whatever the user requires. For numerical calculations, fill it with 0s.
* @method matrix
* @memberof hydro
* @param {Object} params - Contains: m (num columns), n (num rows), d (filler).
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @returns {Object[]} Matrix - m x n array.
* @example
* hydro.analyze.hydro.matrix({params: {m: 3, n: 3, d: 0}});
*/
static matrix({ params, args, data } = {}) {
var mat;
if (typeof params.d === "undefined") {
mat = Array(params.m).map(() => Array(params.n));
} else {
mat = Array(params.m)
.fill(params.d)
.map(() => Array(params.n).fill(params.d));
}
return mat;
}
/**
* Solves linear equations in the form Ax = b.
* @method equationsystemsolver
* @memberof hydro
* @param {Object} params - Contains: right (right hand side 1D JS array), left (left hand side 1D JS array).
* @param {Object} args - Not used by this function.
* @param {Object[]} data - Contains: matrix to be filled.
* @returns {Object[]} Left vector solution.
* @example
* hydro.analyze.hydro.equationsystemsolver({
* params: { left: [0, 0], right: [1, 2] },
* data: [[1, 0], [0, 1]]
* });
*/
static equationsystemsolver({ params, args, data } = {}) {
var matrix = data,
vec_left = params.left,
vec_right = params.right,
fMaxEl,
fAcc,
nodes = vec_left.length;
for (let k = 0; k < nodes; k++) {
//search line with largest element.
fMaxEl = Math.abs(matrix[k][k]);
var m = k;
for (var i = k + 1; i < nodes; i++) {
if (Math.abs(fMaxEl) < Math.abs(matrix[i][k])) {
fMaxEl = matrix[i][k];
m = i;
}
}
// permutation of base line (index k) and max element line (index m)
if (m != k) {
for (var i = k; i < nodes; i++) {
fAcc = matrix[k][i];
matrix[k][i] = matrix[m][i];
matrix[m][i] = fAcc;
}
fAcc = vec_right[k];
vec_right[k] = vec_right[m];
vec_right[m] = fAcc;
}
if (Math.abs(matrix[k][k]) < 1e-10) {
console.log("Singular matrix" + k + " " + matrix[k][k]);
}
for (var j = k + 1; j < nodes; j++) {
fAcc = -matrix[j][k] / matrix[k][k];
for (var i = k; i < nodes; i++) {
matrix[j][i] = matrix[j][i] + fAcc * matrix[k][i];
}
vec_right[j] = vec_right[j] + fAcc * vec_right[k];
}
}
for (var k = nodes - 1; k >= 0; k--) {
vec_left[k] = vec_right[k];
for (var i = k + 1; i < nodes; i++) {
vec_left[k] -= matrix[k][i] * vec_left[i];
}
vec_left[k] = vec_left[k] / matrix[k][k];
}
}
/**
* Calculate the pH value based on the concentration of hydrogen ions (H+).
* pH is a measure of the acidity or alkalinity of a solution, defined as the negative logarithm (base 10)
* of the hydrogen ion concentration. The pH scale ranges from 0 to 14, with pH 7 being neutral.
* Values below 7 indicate acidity, while values above 7 indicate alkalinity.
* @method calculatepH
* @author riya-patil
* @memberof hydro
* @param {Object} params - Contains: hConcentration (hydrogen ion concentration in moles per liter).
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @returns {number} The calculated pH value.
* @example
* // Calculate pH for pure water at 25°C (H+ concentration = 1×10^-7 mol/L)
* hydro.analyze.hydro.calculatepH({ params: { hConcentration: 1e-7 } }); // Returns 7.0
*
* // Calculate pH for an acidic solution (H+ concentration = 1×10^-4 mol/L)
* hydro.analyze.hydro.calculatepH({ params: { hConcentration: 1e-4 } }); // Returns 4.0
*/
static calculatepH({ params, args, data } = {}) {
const { hConcentration } = params;
const pH = -Math.log10(hConcentration);
return pH;
}
/**
* Generates synthetic values based on statistical distributions.
* @method generateSyntheticValue
* @memberof hydro
* @param {Object} params - Contains: mean, stdDev, distributionType ('normal', 'binomial', 'multinomial').
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @returns {number} The random generated synthetic value.
* @example
* hydro.analyze.hydro.generateSyntheticValue({ params: { mean: 10, stdDev: 2, distributionType: 'normal' } });
*/
static generateSyntheticValue({ params, args, data } = {}) {
const { mean, stdDev, distributionType } = params;
let syntheticValue;
//More distributions to be added in next iterations
switch (distributionType) {
case 'normal':
const rand = Math.sqrt(-2 * Math.log(Math.random())) * Math.cos(2 * Math.PI * Math.random());
syntheticValue = mean + stdDev * rand;
break;
case 'binomial':
syntheticValue = stats.binomialDist({ params, args, data });
break;
case 'multinomial':
const multinomialResult = stats.multinomialDistribution({ params, args, data });
syntheticValue = multinomialResult.samples;
break;
default:
throw new Error('Invalid distribution type. Supported types: normal, binomial, multinomial');
}
return syntheticValue;
}
/**
* Converts temperature from one unit to another
* @method convertTemperature
* @author riya-patil
* @memberof hydro
* @param {Object} params sensor_reading (temperature reading from the sensor), from_unit (unit of the input temperature), to_unit (desired unit for conversion).
* @returns {number} The converted temperature value.
* @example
* hydro.analyze.hydro.convertTemperature({ params: { sensor_reading: 25, from_unit: 'Celsius', to_unit: 'Fahrenheit' } });
*/
static convertTemperature({ params, args, data } = {}) {
const { sensor_reading, from_unit, to_unit } = params;
let converted_temperature = 0;
if (from_unit === 'Celsius' && to_unit === 'Fahrenheit') {
converted_temperature = (sensor_reading * 9 / 5) + 32;
} else if (from_unit === 'Fahrenheit' && to_unit === 'Celsius') {
converted_temperature = (sensor_reading - 32) * 5 / 9;
} else if (from_unit === 'Celsius' && to_unit === 'Kelvin') {
converted_temperature = sensor_reading + 273.15;
} else if (from_unit === 'Kelvin' && to_unit === 'Celsius') {
converted_temperature = sensor_reading - 273.15;
} else if (from_unit === 'Fahrenheit' && to_unit === 'Kelvin') {
converted_temperature = (sensor_reading + 459.67) * 5 / 9;
} else if (from_unit === 'Kelvin' && to_unit === 'Fahrenheit') {
converted_temperature = (sensor_reading * 9 / 5) - 459.67;
} else {
// Handle unsupported unit conversions
console.log('Unsupported unit conversion');
}
return converted_temperature;
}
/**
* Calculates Julian day from a given date.
* @method getJulianDay
* @memberof hydro
* @param {Object} params - Contains: date (string 'YYYY-MM-DD' or Date object).
* @param {Object} args - Not used by this function.
* @param {Object} data - Not used by this function.
* @returns {number} Calculated Julian day number.
* @example
* hydro.analyze.hydro.getJulianDay({ params: { date: '2022-01-01' } });
*/
static getJulianDay({ params, args, data } = {}) {
const date = params.date;
let inputDate;
if (typeof date === 'string') {
inputDate = new Date(date);
} else if (date instanceof Date) {
inputDate = date;
} else {
throw new Error('Invalid date format. Expected a string or a Date object.');
}
const year = inputDate.getFullYear();
const month = inputDate.getMonth() + 1;
const day = inputDate.getDate();
// Julian Day Calculation
const a = Math.floor((14 - month) / 12);
const y = year + 4800 - a;
const m = month + 12 * a - 1;
return day + Math.floor((153 * m + 2) / 5) + 365 * y + Math.floor(y / 4) - Math.floor(y / 100) + Math.floor(y / 400) - 32045;
}
/**********************************/
/*** End of Helper functions **/
/**********************************/
}