comparison artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flowdepth/FlowVelocityModelKmValueFinder.java @ 8898:89f3c5462a16

Implemented S-INFO Flowdepth TKH calculation
author mschaefer
date Thu, 22 Feb 2018 12:03:31 +0100
parents
children
comparison
equal deleted inserted replaced
8891:f431aec10d2c 8898:89f3c5462a16
1 /* Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
2 * Software engineering by
3 * Björnsen Beratende Ingenieure GmbH
4 * Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
5 *
6 * This file is Free Software under the GNU AGPL (>=v3)
7 * and comes with ABSOLUTELY NO WARRANTY! Check out the
8 * documentation coming with Dive4Elements River for details.
9 */
10
11 package org.dive4elements.river.artifacts.sinfo.flowdepth;
12
13 import java.util.ArrayList;
14 import java.util.List;
15
16 import org.apache.commons.lang.math.DoubleRange;
17 import org.apache.log4j.Logger;
18 import org.hibernate.SQLQuery;
19 import org.hibernate.Session;
20 import org.hibernate.type.StandardBasicTypes;
21
22 import gnu.trove.TDoubleArrayList;
23
24 import org.dive4elements.river.artifacts.math.Linear;
25 import org.dive4elements.river.artifacts.math.Utils;
26 import org.dive4elements.river.backend.SessionHolder;
27 import org.dive4elements.river.model.River;
28
29 /**
30 * Searchable sorted km array with parallel FlowVelocityKmModelValues array and linear interpolation for km and the model values between the array elements.<br />
31 * {@link loadValues} loads all the model values for a given km range of a river.<br />
32 * {@link findKmQValues} then searches a km in the values table or the nearest including km interval, resp.
33 * The v and tau values for a given discharge are either found directly or also interpolated linearly.<br />
34 *
35 * (Created based on a copy of FlowVelocityMeasurementFactory.)
36 *
37 * @author Matthias Schäfer
38 *
39 */
40 public class FlowVelocityModelKmValueFinder
41 {
42 /***** FIELDS *****/
43
44 /**
45 * Private log to use here.
46 */
47 private static Logger log = Logger.getLogger(FlowVelocityModelKmValueFinder.class);
48
49 /**
50 * Query for a range of stations of a river with all their q, main-v and tau values.<br />
51 * (Might be several 10000 rows if many stations and large q range)
52 */
53 private static final String SQL_SELECT_ALL =
54 "SELECT fvmv.station AS station, fvmv.q AS q, fvmv.main_channel AS vmain, fvmv.shear_stress AS tau"
55 + " FROM (discharge_zone dz INNER JOIN flow_velocity_model fvm ON dz.id = fvm.discharge_zone_id)"
56 + " INNER JOIN flow_velocity_model_values fvmv ON fvm.id = fvmv.flow_velocity_model_id"
57 + " WHERE (dz.river_id = :river_id) AND (fvmv.station BETWEEN :kmfrom - 0.0001 AND :kmto + 0.0001)"
58 /* + " WHERE (dz.river_id = :river_id) AND (fvmv.q BETWEEN :qmin AND :qmax)" */
59 + " ORDER BY fvmv.station ASC, fvmv.q ASC";
60
61 /**
62 * Query for a river's max km below a limit with all its q, main-v and tau values.
63 */
64 private static final String SQL_SELECT_KMLOWER =
65 "SELECT fvmv.station AS station, fvmv.q AS q, fvmv.main_channel AS vmain, fvmv.shear_stress AS tau"
66 + " FROM flow_velocity_model_values fvmv"
67 + " INNER JOIN (SELECT MAX(fvmvi.station) AS kmmax"
68 + " FROM(discharge_zone dz INNER JOIN flow_velocity_model fvm ON dz.id = fvm.discharge_zone_id)"
69 + " INNER JOIN flow_velocity_model_values fvmvi ON fvm.id = fvmvi.flow_velocity_model_id"
70 + " WHERE (dz.river_id = :river_id) AND (fvmvi.station < :kmfrom - 0.0001)) finf ON fvmv.station = finf.kmmax"
71 + " ORDER BY fvmv.q ASC";
72
73 /**
74 * Query for a river's min km above a limit with all its q, main-v and tau values.
75 */
76 private static final String SQL_SELECT_KMUPPER =
77 "SELECT fvmv.station AS station, fvmv.q AS q, fvmv.main_channel AS vmain, fvmv.shear_stress AS tau"
78 + " FROM flow_velocity_model_values fvmv"
79 + " INNER JOIN (SELECT MIN(fvmvi.station) AS kmmin"
80 + " FROM(discharge_zone dz INNER JOIN flow_velocity_model fvm ON dz.id = fvm.discharge_zone_id)"
81 + " INNER JOIN flow_velocity_model_values fvmvi ON fvm.id = fvmvi.flow_velocity_model_id"
82 + " WHERE (dz.river_id = :river_id) AND (fvmvi.station > :kmto + 0.0001)) fsup ON fvmv.station = fsup.kmmin"
83 + " ORDER BY fvmv.q ASC";
84
85 /**
86 * Query to select all km-q-v-tau of a river that are the q maxima below a q limit
87 */
88 private static final String SQL_SELECT_QLOWER =
89 "SELECT fvmv.station AS station, fvmv.q AS q, fvmv.main_channel AS vmain, fvmv.shear_stress AS tau"
90 + " FROM flow_velocity_model_values fvmv"
91 + " INNER JOIN (SELECT fv2.station, MAX(fv2.q) AS q"
92 + " FROM (discharge_zone dz INNER JOIN flow_velocity_model fvm ON dz.id = fvm.discharge_zone_id)"
93 + " INNER JOIN flow_velocity_model_values fv2 ON fvm.id = fv2.flow_velocity_model_id"
94 + " WHERE (dz.river_id = :river_id) AND (fv2.q < :qlim) GROUP BY fv2.station) qx"
95 + " ON (fvmv.station=qx.station) AND (fvmv.q=qx.q)"
96 + " ORDER BY fvmv.station ASC";
97
98 /**
99 * Query to select all km-q-v-tau of a river that are the q minima above a q limit
100 */
101 private static final String SQL_SELECT_QUPPER =
102 "SELECT fvmv.station AS station, fvmv.q AS q, fvmv.main_channel AS vmain, fvmv.shear_stress AS tau"
103 + " FROM flow_velocity_model_values fvmv"
104 + " INNER JOIN (SELECT fv2.station, MIN(fv2.q) AS q"
105 + " FROM (discharge_zone dz INNER JOIN flow_velocity_model fvm ON dz.id = fvm.discharge_zone_id)"
106 + " INNER JOIN flow_velocity_model_values fv2 ON fvm.id = fv2.flow_velocity_model_id"
107 + " WHERE (dz.river_id = :river_id) AND (fv2.q > :qlim) GROUP BY fv2.station) qx"
108 + " ON (fvmv.station=qx.station) AND (fvmv.q=qx.q)"
109 + " ORDER BY fvmv.station ASC";
110
111 /**
112 * Kms of the loaded river range
113 */
114 private TDoubleArrayList kms;
115
116 /**
117 * For each km in kms a list of q-v-tau-tupels
118 */
119 private List<FlowVelocityKmModelValues> values;
120
121 /**
122 * Searched km of the last findKmValue
123 */
124 private double findKm;
125
126 /**
127 * kms and values index of the interval start found by the last findKmValue
128 */
129 private int leftIndexFound = -1;
130
131 /**
132 * kms and values index of the interval end found by the last findKmValue
133 */
134 private int rightIndexFound = -1;
135
136 /**
137 * Q of the last findKmQValues
138 */
139 private double findQ;
140
141 /***** METHODS *****/
142
143 /**
144 * Discharge of the last {@link findKmQValue}
145 */
146 public double getFindQ() {
147 return findQ;
148 }
149
150 /**
151 * Velocity of the last {@link findKmQValues}
152 */
153 public double getFindVmainFound() {
154 if (leftIndexFound < 0)
155 return Double.NaN;
156 else if (leftIndexFound == rightIndexFound)
157 return getLeftValues().getVmainFound();
158 else
159 return Linear.linear(findKm, getLeftValues().getKm(), getRightValues().getKm(), getLeftValues().getVmainFound(), getRightValues().getVmainFound());
160 }
161
162 /**
163 * Shear stress tau of the last {@link findKmQValues}
164 */
165 public double getFindTauFound() {
166 if (leftIndexFound < 0)
167 return Double.NaN;
168 else if (leftIndexFound == rightIndexFound)
169 return getLeftValues().getTauFound();
170 else
171 return Linear.linear(findKm, getLeftValues().getKm(), getRightValues().getKm(), getLeftValues().getTauFound(), getRightValues().getTauFound());
172 }
173
174 /**
175 * Whether the discharge has been interpolated in the last {@link findKmQValues}
176 */
177 public boolean getFindIsQInterpolated() {
178 return (getLeftValues() != null) && (getLeftValues().getIsInterpolated() || getRightValues().getIsInterpolated());
179 }
180
181 /**
182 * Queries a range of a river's kms with all their q-v-tau values.
183 * @return Whether the load has been successful
184 */
185 @SuppressWarnings("unchecked")
186 public boolean loadValues(River river, DoubleRange kmRange, DoubleRange qRange) {
187 // DB session
188 log.debug(String.format("loadValues km %.3f - %.3f / q %.1f - %.1f", kmRange.getMinimumDouble(), kmRange.getMaximumDouble(), qRange.getMinimumDouble(), qRange.getMaximumDouble()));
189 kms = new TDoubleArrayList();
190 values = new ArrayList<FlowVelocityKmModelValues>();
191 boolean isDemoValuesCorrection = river.getName().equalsIgnoreCase("beispielfluss");
192 final Session session = SessionHolder.HOLDER.get();
193
194 // Select km infimum
195 SQLQuery sqlQuery = session.createSQLQuery(SQL_SELECT_KMLOWER)
196 .addScalar("station", StandardBasicTypes.DOUBLE)
197 .addScalar("q", StandardBasicTypes.DOUBLE)
198 .addScalar("vmain", StandardBasicTypes.DOUBLE)
199 .addScalar("tau", StandardBasicTypes.DOUBLE);
200 sqlQuery.setParameter("river_id", river.getId());
201 sqlQuery.setParameter("kmfrom", kmRange.getMinimumDouble());
202 addKms(sqlQuery.list(), isDemoValuesCorrection);
203
204 // Select km range
205 sqlQuery = session.createSQLQuery(SQL_SELECT_ALL)
206 .addScalar("station", StandardBasicTypes.DOUBLE)
207 .addScalar("q", StandardBasicTypes.DOUBLE)
208 .addScalar("vmain", StandardBasicTypes.DOUBLE)
209 .addScalar("tau", StandardBasicTypes.DOUBLE);
210 sqlQuery.setParameter("river_id", river.getId());
211 sqlQuery.setParameter("kmfrom", kmRange.getMinimumDouble());
212 sqlQuery.setParameter("kmto", kmRange.getMaximumDouble());
213 //sqlQuery.setParameter("qmin", qRange.getMinimumDouble());
214 //sqlQuery.setParameter("qmax", qRange.getMaximumDouble());
215 int kmcount = kms.size();
216 int rowcount = addKms(sqlQuery.list(), isDemoValuesCorrection);
217 kmcount = kms.size() - kmcount;
218
219 // Select km supremum
220 sqlQuery = session.createSQLQuery(SQL_SELECT_KMUPPER)
221 .addScalar("station", StandardBasicTypes.DOUBLE)
222 .addScalar("q", StandardBasicTypes.DOUBLE)
223 .addScalar("vmain", StandardBasicTypes.DOUBLE)
224 .addScalar("tau", StandardBasicTypes.DOUBLE);
225 sqlQuery.setParameter("river_id", river.getId());
226 sqlQuery.setParameter("kmto", kmRange.getMaximumDouble());
227 int supcnt = addKms(sqlQuery.list(), isDemoValuesCorrection);
228
229 // Add copy of last km for search of max km value
230 if ((supcnt == 0) && (values.size() >= 1)) {
231 kms.add(kms.getQuick(kms.size()) + 0.0001);
232 values.add(new FlowVelocityKmModelValues(kms.getQuick(kms.size()-1), values.get(values.size()-1)));
233 }
234
235 // log.debug
236 if (values.size() - 1 >= 0) {
237 log.debug(String.format("loadValues %d: km %.3f - %d values", 0, values.get(0).getKm(), values.get(0).size()));
238 if (values.size() - 1 >= 1) {
239 log.debug(String.format("loadValues %d: km %.3f - %d values", 1, values.get(1).getKm(), values.get(1).size()));
240 if (values.size() - 1 >= 2)
241 log.debug("loadValues ...");
242 if (values.size() - 2 >= 3)
243 log.debug(String.format("loadValues %d: km %.3f - %d values", values.size()-2, values.get(values.size()-2).getKm(), values.get(values.size()-2).size()));
244 if (values.size() - 1 >= 3)
245 log.debug(String.format("loadValues %d: km %.3f - %d values", values.size()-1, values.get(values.size()-1).getKm(), values.get(values.size()-1).size()));
246 }
247 }
248 log.debug(String.format("loadValues %d kms, %d values loaded", kmcount, rowcount));
249 return (kms.size() >= 1);
250 }
251
252 /**
253 * Adds the km-q-v-tau values of a query result row to the last km of the list, or a new one resp.
254 * @return Number of rows
255 */
256 private int addKms(List<Object[]> rows, boolean isDemoValuesCorrection) {
257 for (Object[] row : rows) {
258 if ((kms.size() == 0) || !Utils.epsilonEquals(kms.get(kms.size()-1), (double) row[0], 0.0001)) {
259 kms.add((double) row[0]);
260 values.add(new FlowVelocityKmModelValues(kms.get(kms.size()-1)));
261 }
262 if (isDemoValuesCorrection)
263 // "Verfremdung" der v-Werte etwas korrigieren (Originalwerte wurden mit Zufallswert zwischen 10 und 20 multipliziert)
264 values.get(values.size()-1).addValues((double) row[1], ((double) row[2]) / 10, (double) row[3]);
265 else
266 values.get(values.size()-1).addValues((double) row[1], (double) row[2], (double) row[3]);
267 }
268 return rows.size();
269 }
270
271 /**
272 * Searches a km and finds or interpolates the velocity and shear stress values for a discharge<br />
273 * The values may be got via {@link getVmainFound} etc.
274 * @return Whether values have been found
275 */
276 public boolean findKmQValues(double km, double q) {
277 findQ = q;
278 if (!searchKm(km))
279 return false;
280 if (leftIndexFound == rightIndexFound) {
281 // Exact km match
282 final double qfound = getLeftValues().findQ(q);
283 log.debug(String.format("findKmQValues km %.3f q %.0f = %.0f (%d)", km, q, qfound, leftIndexFound));
284 return !Double.isNaN(qfound);
285 }
286 else {
287 final double[] qfound = {getLeftValues().findQ(q), getRightValues().findQ(q)};
288 log.debug(String.format("findKmQValues km %.3f q %.0f = %.0f (%d, %.3f) - %.0f (%d, %.3f)", km, q, qfound[0], leftIndexFound,
289 getLeftValues().getKm(), qfound[1], rightIndexFound, getRightValues().getKm()));
290 return !Double.isNaN(qfound[0]) && !Double.isNaN(qfound[1]);
291 }
292 }
293
294 /**
295 * Searches a km
296 * @return Whether the km was within the supported range
297 */
298 private boolean searchKm(double km) {
299 findKm = km;
300 leftIndexFound = -1;
301 rightIndexFound = -1;
302 if ((kms == null) || (kms.size() == 0))
303 return false;
304 int i = kms.binarySearch(km);
305 if (i >= 0) {
306 // Exact km match
307 leftIndexFound = i;
308 rightIndexFound = i;
309 return true;
310 }
311 else {
312 // Out of range or within km interval
313 if (i < 0)
314 i = -i - 1;
315 if ((i <= 0) || (i >= kms.size()))
316 return false;
317 leftIndexFound = i - 1;
318 rightIndexFound = i;
319 return true;
320 }
321 }
322
323 private FlowVelocityKmModelValues getLeftValues() {
324 return values.get(leftIndexFound);
325 }
326 private FlowVelocityKmModelValues getRightValues() {
327 return values.get(rightIndexFound);
328 }
329
330 }

http://dive4elements.wald.intevation.org