Mercurial > dive4elements > river
changeset 343:f165c7d5d6db
Implementation of the "Ruecksprungkorrektur" to be done in "W fuer angepassten Abflusslaengschnitt".
flys-artifacts/trunk@1743 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Fri, 22 Apr 2011 09:57:27 +0000 (2011-04-22) |
parents | f72c63713099 |
children | 79401797f4e1 |
files | flys-artifacts/ChangeLog flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java flys-artifacts/src/main/java/de/intevation/flys/artifacts/states/RiverSelect.java |
diffstat | 3 files changed, 219 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/flys-artifacts/ChangeLog Thu Apr 21 08:49:58 2011 +0000 +++ b/flys-artifacts/ChangeLog Fri Apr 22 09:57:27 2011 +0000 @@ -1,3 +1,34 @@ +2011-04-22 Sascha L. Teichmann <sascha.teichmann@intevation.de> + + * src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java: + + Implementation of the "Ruecksprungkorrektur" to be done in + "W fuer angepassten Abflusslaengschnitt". + + All tests show the expected results. In some corner cases the + algorithm described in the "Anwenderhandbuch" chapter 3.3.4.3 "Korrektur" + has some definition shortcomings: + + a - What should happend when you cannot find point 2 because + you cannot step back one quarter from point 3 because there + is no data there any more (river too short in this direction)? + The implemented algorithm raises point 3' only to an + according factor. E.g. If you can step back the whole quarter + distance the elevation is the full quarter. If you can + step back only the half of the quarter the elevation is + only an eighth. + + b - If the water heights between point 2 and 3 are constant then + the algorithm will produce a spline interpolation that + lowers those values. Is this intended? + + For real data the back jumps are expected to be more in the middle + of the distance ranges so the corner cases are maybe not so + important. + + * src/main/java/de/intevation/flys/artifacts/states/RiverSelect.java: + Removed superfluous import. + 2011-04-21 Ingo Weinzierl <ingo@intevation.de> * src/main/java/de/intevation/flys/themes/ThemeFactory.java: Removed debug
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/math/BackJumpCorrector.java Fri Apr 22 09:57:27 2011 +0000 @@ -0,0 +1,188 @@ +package de.intevation.flys.artifacts.math; + +import java.util.ArrayList; +import java.util.List; + +import java.io.Serializable; + +import org.apache.commons.math.analysis.interpolation.SplineInterpolator; + +import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; + +import org.apache.commons.math.ArgumentOutsideDomainException; + +import org.apache.log4j.Logger; + +public class BackJumpCorrector +implements Serializable +{ + private static Logger log = Logger.getLogger(BackJumpCorrector.class); + + protected ArrayList<Double> backjumps; + + protected double [] corrected; + + public BackJumpCorrector() { + backjumps = new ArrayList<Double>(); + } + + public boolean hasBackJumps() { + return !backjumps.isEmpty(); + } + + public List<Double> getBackJumps() { + return backjumps; + } + + public double [] getCorrected() { + return corrected; + } + + public boolean doCorrection(double [] km, double [] ws) { + int N = km.length; + + if (N != ws.length) { + throw new IllegalArgumentException("km.length != ws.length"); + } + + if (N < 2) { + return false; + } + + SplineInterpolator interpolator = null; + + for (int i = 1; i < N; ++i) { + if (ws[i] <= ws[i-1]) { + // no back jump + continue; + } + backjumps.add(km[i]); + + if (corrected == null) { + // lazy cloning + ws = corrected = (double [])ws.clone(); + } + + int back = i-2; + double distance = Math.abs(km[i] - km[i-1]); + double rest = 0.0; + double ikm = 0.0; + + while (back > -1) { + if (ws[back] < ws[i]) { // under water + // continue scanning backwards + distance += Math.abs(km[back+1] - km[back]); + --back; + continue; + } + if (ws[back] > ws[i]) { // over water + // need to calculate intersection + log.debug("intersection case"); + double m = (km[back+1]-km[back])/(ws[back+1]-ws[back]); + double b = km[back]-ws[back]*m; + ikm = m*ws[i] + b; + distance += Math.abs(ikm - km[back+1]); + rest = Math.abs(km[back] - ikm); + } + else { + // equals height at ws[back] + log.debug("exact height match"); + distance += Math.abs(km[back+1] - km[back]); + ikm = km[back]; + } + break; + } + + if (back < 0) { + log.debug("run over left border"); + // fill with same as ws[i] + for (int j = 0; j < i; ++j) { + ws[j] = ws[i]; + } + continue; + } + + double quarterDistance = 0.25*distance; + + // Now further seek back for the max height + + double restDistance = Math.max(0.0, quarterDistance - rest); + --back; + + double mkmw = ws[i]; + double mkm = km[0]; + + while (back > -1) { + double d = Math.abs(km[back+1] - km[back]); + restDistance -= d; + if (restDistance > 0.0) { + --back; + continue; + } + if (restDistance < 0.0) { + // need to calculate intersection + if (km[back] == km[back+1]) { // should not happen + mkm = km[back]; + mkmw = 0.5*(ws[back] + ws[back+1]); + } + else { + double m = (ws[back+1]-ws[back])/(km[back+1]-km[back]); + double b = ws[back] - km[back]*m; + mkm = km[back] + restDistance; + mkmw = m*mkm + b; + } + } + else { + // exact match + mkm = km[back]; + mkmw = ws[back]; + } + break; + } + + double factor = back >= 0 && Math.abs(restDistance) < 1e-4 + ? 1.0 + : 1.0 - Math.min(1, Math.max(0, restDistance/quarterDistance)); + + double ikmw = factor*0.25*(mkmw-ws[i]) + ws[i]; + + double end = ikm + quarterDistance*factor; + + double [] x = { mkm, ikm, end }; + double [] y = { mkmw, ikmw, ws[i] }; + + if (interpolator == null) { + interpolator = new SplineInterpolator(); + } + + PolynomialSplineFunction spline = interpolator.interpolate(x, y); + + try { + if (log.isDebugEnabled()) { + log.debug("spline points:"); + for (int j = 0; j < x.length; ++j) { + log.debug(x[j] + " " + y[j] + " " + spline.value(x[j])); + } + } + + for (back = Math.max(back, 0); + back < i && km[back] < end; + ++back + ) { + // to 3/4 point fill with spline values + ws[back] = spline.value(km[back]); + } + while (back < i) { + // fill the rest with ws[i] + ws[back++] = ws[i]; + } + } + catch (ArgumentOutsideDomainException aode) { + log.error("spline interpolation failed", aode); + } + } // for all km + + return !backjumps.isEmpty(); + } +} +// vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/states/RiverSelect.java Thu Apr 21 08:49:58 2011 +0000 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/states/RiverSelect.java Fri Apr 22 09:57:27 2011 +0000 @@ -18,7 +18,6 @@ import de.intevation.flys.model.River; -import de.intevation.flys.artifacts.FLYSArtifact; import de.intevation.flys.artifacts.model.RiverFactory; import de.intevation.flys.artifacts.resources.Resources;