# HG changeset patch # User Sascha L. Teichmann # Date 1303309611 0 # Node ID 4509ba8fae684389f39cb1aa37c1f942ee0a54e1 # Parent cf84f0f926e9bafe9cb598fd89f8db9cdef3742f Implementation of "Abflusskurve/Abflusstafel" calculation. flys-artifacts/trunk@1738 c6561f87-3c4e-4783-a992-168aeb5c3f6f diff -r cf84f0f926e9 -r 4509ba8fae68 flys-artifacts/ChangeLog --- a/flys-artifacts/ChangeLog Wed Apr 20 13:13:19 2011 +0000 +++ b/flys-artifacts/ChangeLog Wed Apr 20 14:26:51 2011 +0000 @@ -1,3 +1,17 @@ +2011-04-20 Sascha L. Teichmann + + * src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java: + Implementation of "Abflusskurve/Abflusstafel" calculation. + + Added method interpolateWQ() which takes an km and results in a + tuple of two double arrays containing the w/q values interpolated + between the surrounding w/q values of the table. + w values are interpolated linear, q values with a cubic spline. + + Drawing w over q gives you the discharge table at the given km. + + !!! This code needs testing !!! + 2011-04-20 Sascha L. Teichmann * pom.xml: Added dependency to Apache Commons Math 2.2 (Apache License 2.0) @@ -15,7 +29,7 @@ 2011-04-19 Sascha L. Teichmann * src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java: - add a method interpolateW() method which takes an array of + add a method interpolateW() which takes an array of q values and returns an equal sized array of w values. This is essentially the "Wasserstand/Wasserspiegellagen" calculation of desktop FLYS. diff -r cf84f0f926e9 -r 4509ba8fae68 flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java --- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java Wed Apr 20 13:13:19 2011 +0000 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java Wed Apr 20 14:26:51 2011 +0000 @@ -15,6 +15,14 @@ import java.util.Collections; import java.util.Iterator; +import org.apache.log4j.Logger; + +import org.apache.commons.math.analysis.interpolation.SplineInterpolator; + +import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; + +import org.apache.commons.math.ArgumentOutsideDomainException; + import org.hibernate.Session; import org.hibernate.Query; import org.hibernate.SQLQuery; @@ -24,6 +32,7 @@ public class WstValueTable implements Serializable { + private static Logger log = Logger.getLogger(WstValueTable.class); // TODO: put this into a property file public static final String SQL_POS_WQ = @@ -31,6 +40,8 @@ " FROM wst_value_table" + " WHERE wst_id = :wst_id"; + public static final double DEFAULT_STEP_WIDTH = 0.01; + public static class Column implements Serializable { @@ -188,6 +199,134 @@ if (d > 0.0) return +1; return 0; } + + public double [][] cloneWQs() { + return new double [][] { + (double [])ws.clone(), + (double [])ws.clone() }; + } + + public double [][] interpolateWQ(Row other, double km, double stepWidth) { + + int W1 = ascendingWs(); + int W2 = other.ascendingWs(); + + int W = Math.min(W1, W2); + + if (W < 1) { + return new double[2][0]; + } + + double factor = factor(km, this.km, other.km); + + double minW = weight(factor, ws[0], other.ws[0]); + double maxW = weight(factor, ws[W-1], other.ws[W-1]); + + if (minW > maxW) { + double t = minW; + minW = maxW; + maxW = t; + } + double [] x = new double[W]; + double [] y = new double[W]; + + for (int i = 0; i < W; ++i) { + x[i] = weight(factor, ws[i], other.ws[i]); + y[i] = weight(factor, qs[i], other.qs[i]); + } + + SplineInterpolator interpolator = new SplineInterpolator(); + PolynomialSplineFunction spline = interpolator.interpolate(x, y); + + double [] outWs = + new double[(int)Math.ceil((maxW - minW)/stepWidth)]; + double [] outQs = + new double[outWs.length]; + + try { + double w = minW; + for (int i = 0; i < outWs.length; ++i, w += stepWidth) { + outQs[i] = spline.value(outWs[i] = w); + } + } + catch (ArgumentOutsideDomainException aode) { + log.error("Spline interpolation failed.", aode); + } + + return new double [][] { outWs, outQs }; + + } + + public double [][] interpolateWQ(double stepWidth) { + int W = ascendingWs(); // ignore back jumps + + if (W < 1) { + return new double[2][0]; + } + + double [] x = new double[W]; + double [] y = new double[W]; + + for (int i = 0; i < W; ++i) { + x[i] = ws[i]; + y[i] = qs[i]; + } + + SplineInterpolator interpolator = new SplineInterpolator(); + + PolynomialSplineFunction spline = interpolator.interpolate(x, y); + + double minW = ws[0]; + double maxW = ws[W-1]; + + double [] outWs = + new double[(int)Math.ceil((maxW - minW)/stepWidth)]; + double [] outQs = + new double[outWs.length]; + + try { + double w = minW; + for (int i = 0; i < outWs.length; ++i, w += stepWidth) { + outQs[i] = spline.value(outWs[i] = w); + } + } + catch (ArgumentOutsideDomainException aode) { + log.error("Spline interpolation failed.", aode); + } + + return new double [][] { outWs, outQs }; + } + + public int ascendingWs() { + if (ws.length < 2) { + return ws.length; + } + + int idx = 1; + + for (; idx < ws.length; ++idx) { + if (Double.isNaN(ws[idx]) || ws[idx] < ws[idx-1]) { + return idx; + } + } + + return idx; + } + + public double [][] weightWQs(Row other, double km) { + int W = Math.min(ws.length, other.ws.length); + double factor = factor(km, this.km, other.km); + + double [] outWs = new double[W]; + double [] outQs = new double[W]; + + for (int i = 0; i < W; ++i) { + outWs[i] = weight(factor, ws[i], other.ws[i]); + outQs[i] = weight(factor, qs[i], other.qs[i]); + } + + return new double [][] { outWs, outQs }; + } } // class Row @@ -234,6 +373,40 @@ return ws; } + + public double [][] interpolateWQ(double km) { + return interpolateWQ(km, DEFAULT_STEP_WIDTH, true); + } + + public double [][] interpolateWQ( + double km, + double stepWidth, + boolean checkAscending + ) { + int rowIndex = Collections.binarySearch(rows, new Row(km)); + + if (rowIndex >= 0) { // direct row match + Row row = rows.get(rowIndex); + return checkAscending + ? row.interpolateWQ(stepWidth) + : row.cloneWQs(); + } + + rowIndex = -rowIndex -1; + + if (rowIndex < 1 || rowIndex >= rows.size()) { + // do not extrapolate + return new double[2][0]; + } + + Row r1 = rows.get(rowIndex-1); + Row r2 = rows.get(rowIndex); + + return checkAscending + ? r1.interpolateWQ(r2, km, stepWidth) + : r1.weightWQs(r2, km); + } + public static WstValueTable getTable(River river) { return getTable(river, 0); }