Mercurial > dive4elements > river

comparison artifacts/src/main/java/org/dive4elements/river/artifacts/sinfo/flowdepth/BedQualityD50KmValueFinder.java @ 8917:86395ab8ebc3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
SINFO Flowdepth: d50 aggregation by median instead of arithmetic mean
author mschaefer
date Mon, 26 Feb 2018 18:56:05 +0100
parents 37ff7f435912
children
comparison
equal deleted inserted replaced
8911:37ff7f435912 8917:86395ab8ebc3
8 * documentation coming with Dive4Elements River for details. 8 * documentation coming with Dive4Elements River for details.
9 */ 9 */
10 10
11 package org.dive4elements.river.artifacts.sinfo.flowdepth; 11 package org.dive4elements.river.artifacts.sinfo.flowdepth;
12 12
13 import java.util.Date;
14 import java.util.List; 13 import java.util.List;
15 14
16 import org.apache.commons.lang.math.DoubleRange; 15 import org.apache.commons.lang.math.DoubleRange;
17 import org.apache.commons.math.ArgumentOutsideDomainException; 16 import org.apache.commons.math.ArgumentOutsideDomainException;
18 import org.apache.commons.math.analysis.interpolation.LinearInterpolator; 17 import org.apache.commons.math.analysis.interpolation.LinearInterpolator;
20 import org.apache.log4j.Logger; 19 import org.apache.log4j.Logger;
21 import org.dive4elements.river.artifacts.math.Utils; 20 import org.dive4elements.river.artifacts.math.Utils;
22 import org.dive4elements.river.artifacts.model.DateRange; 21 import org.dive4elements.river.artifacts.model.DateRange;
23 import org.dive4elements.river.backend.SedDBSessionHolder; 22 import org.dive4elements.river.backend.SedDBSessionHolder;
24 import org.dive4elements.river.model.River; 23 import org.dive4elements.river.model.River;
25 import org.dive4elements.river.utils.DoubleUtil;
26 import org.hibernate.SQLQuery; 24 import org.hibernate.SQLQuery;
27 import org.hibernate.Session; 25 import org.hibernate.Session;
28 import org.hibernate.type.StandardBasicTypes; 26 import org.hibernate.type.StandardBasicTypes;
29 27
30 import gnu.trove.TDoubleArrayList; 28 import gnu.trove.TDoubleArrayList;
37 * @author Matthias Schäfer 35 * @author Matthias Schäfer
38 * 36 *
39 */ 37 */
40 public class BedQualityD50KmValueFinder { 38 public class BedQualityD50KmValueFinder {
41 39
42 /***** INNER CLASSES *****/
43
44 /**
45 * A bed measurements aggregate with its d50 characteristic grain diameter
46 */
47 private class D50Measurement {
48 private double km;
49 public double getKm() {
50 return km;
51 }
52 private Date mindate;
53 public Date getMinDate() {
54 return mindate;
55 }
56 private Date maxdate;
57 public Date getMaxDate() {
58 return maxdate;
59 }
60 private int cnt;
61 public int getCnt() {
62 return cnt;
63 }
64 private double mindepth;
65 public double getMinDepth() {
66 return mindepth;
67 }
68 private double maxdepth;
69 public double getMaxDepth() {
70 return maxdepth;
71 }
72 private double d50;
73 /**
74 * D50 in m
75 */
76 public double getD50() {
77 return d50;
78 }
79 /**
80 * Parameter constructor
81 */
82 public D50Measurement(double km, Date mindate, Date maxdate, int cnt, double mindepth, double maxdepth, double d50mm) {
83 this.km = km;
84 this.mindate = mindate;
85 this.maxdate = maxdate;
86 this.cnt = cnt;
87 this.mindepth = mindepth;
88 this.maxdepth = maxdepth;
89 this.d50 = d50mm / 1000;
90 }
91
92 /**
93 * Query result row constructor
94 */
95 public D50Measurement(Object[] tuple, String[] aliases) {
96 km = 0;
97 mindate = null;
98 maxdate = null;
99 cnt = 0;
100 mindepth = Double.NaN;
101 maxdepth = Double.NaN;
102 d50 = Double.NaN;
103 for (int i = 0; i < tuple.length; ++i) {
104 if (tuple[i] == null)
105 continue;
106 switch (aliases[i]) {
107 case "km":
108 km = ((Number) tuple[i]).doubleValue();
109 break;
110 case "mindate":
111 mindate = (Date) tuple[i];
112 break;
113 case "maxdate":
114 maxdate = (Date) tuple[i];
115 break;
116 case "cnt":
117 cnt = ((Number) tuple[i]).intValue();
118 break;
119 case "mindepth":
120 mindepth = ((Number) tuple[i]).doubleValue();
121 break;
122 case "maxdepth":
123 maxdepth = ((Number) tuple[i]).doubleValue();
124 break;
125 case "d50":
126 d50 = ((Number) tuple[i]).doubleValue() / 1000; // mm to m
127 break;
128 default:
129 break;
130 }
131 }
132 }
133 }
134
135 /***** FIELDS *****/ 40 /***** FIELDS *****/
136 41
137 /** 42 /**
138 * Private log to use here. 43 * Private log to use here.
139 */ 44 */
140 private static Logger log = Logger.getLogger(BedQualityD50KmValueFinder.class); 45 private static Logger log = Logger.getLogger(BedQualityD50KmValueFinder.class);
141 46
142 /** 47 /**
143 * Query that aggregates by km for a km range and a time period all sub layer bed measurements with their d50<br /> 48 * Query selecting all sub layer bed measurements with their d50 for a km range and a time period<br />
144 * <br /> 49 * <br />
145 * A km may have bed measurements for multiple dates, multiple distances from the river bank, and multiple depth layers. 50 * A km may have bed measurements for multiple dates, multiple distances from the river bank, and multiple depth layers.
146 * The query filters by km range, time period and layer (sub layer: below bed to max. 50 cm depth). 51 * The query filters by km range, time period and layer (sub layer: below bed to max. 50 cm depth).<br />
147 * Those measurements are then grouped by km, and the D50 aggregated as average value. 52 *
53 * If PostgreSQL would support a median aggregate function like Oracle does the aggregation could be placed into this query.
148 */ 54 */
149 private static final String SQL_BED_D50_SUBLAYER_MEASUREMENT = 55 private static final String SQL_BED_D50_SUBLAYER_MEASUREMENT =
150 "SELECT t.km, MIN(t.datum) AS mindate, MAX(t.datum) AS maxdate, COUNT(*) AS cnt," 56 "SELECT t.km, t.datum, p.tiefevon, p.tiefebis, a.d50 AS d50"
151 + " MIN(p.tiefevon) AS mindepth, MAX(p.tiefebis) AS maxdepth, AVG(a.d50) AS d50" 57 + " FROM sohltest t INNER JOIN station s ON t.stationid = s.stationid"
152 + " FROM sohltest t INNER JOIN station s ON t.stationid = s.stationid" 58 + " INNER JOIN gewaesser g ON s.gewaesserid = g.gewaesserid"
153 + " INNER JOIN gewaesser g ON s.gewaesserid = g.gewaesserid" 59 + " INNER JOIN sohlprobe p ON t.sohltestid = p.sohltestid"
154 + " INNER JOIN sohlprobe p ON t.sohltestid = p.sohltestid" 60 + " INNER JOIN siebanalyse a ON p.sohlprobeid = a.sohlprobeid"
155 + " INNER JOIN siebanalyse a ON p.sohlprobeid = a.sohlprobeid" 61 + " WHERE (g.name = :name) AND (s.km BETWEEN :fromkm - 0.0001 AND :tokm + 0.0001)"
156 + " WHERE (g.name = :name) AND (s.km BETWEEN :fromkm - 0.0001 AND :tokm + 0.0001)" 62 + " AND (p.tiefevon > 0.0) AND (p.tiefebis <= 0.5)"
157 + " AND (p.tiefevon > 0.0) AND (p.tiefebis <= 0.5)" 63 + " AND (t.datum BETWEEN :fromdate AND :todate)"
158 + " AND (t.datum BETWEEN :fromdate AND :todate)" 64 + " ORDER BY t.km ASC, a.d50 ASC";
159 + " GROUP BY t.km"
160 + " ORDER BY t.km";
161 private static final String[] SQL_BED_D50_SELECT_ALIAS = {"km", "mindate", "maxdate", "cnt", "mindepth", "maxdepth", "d50"};
162 65
163 /** 66 /**
164 * Real linear interpolator for kms and d50 values 67 * Real linear interpolator for kms and d50 values (m)
165 */ 68 */
166 private PolynomialSplineFunction interpolator; 69 private PolynomialSplineFunction interpolator;
167 70
168 71
169 /***** METHODS *****/ 72 /***** METHODS *****/
170 73
171 /** 74 /**
172 * Returns the d50 value interpolated according to a km 75 * Returns the d50 value interpolated according to a km
173 * @return d50 (mm) of the km, or NaN 76 * @return d50 (m) of the km, or NaN
174 */ 77 */
175 public double findD50(double km) throws ArgumentOutsideDomainException { 78 public double findD50(double km) throws ArgumentOutsideDomainException {
176 return interpolator.value(km); 79 return interpolator.value(km);
177 /* ohne interpolation: 80 /* ohne interpolation:
178 if ((kms == null) || (kms.size() == 0)) 81 if ((kms == null) || (kms.size() == 0))
195 */ 98 */
196 public boolean loadValues(final River river, final DoubleRange kmRange, final DateRange dateRange) { 99 public boolean loadValues(final River river, final DoubleRange kmRange, final DateRange dateRange) {
197 log.debug(String.format("loadValues km %.3f - %.3f %tF - %tF", kmRange.getMinimumDouble(), kmRange.getMaximumDouble(), dateRange.getFrom(), dateRange.getTo())); 100 log.debug(String.format("loadValues km %.3f - %.3f %tF - %tF", kmRange.getMinimumDouble(), kmRange.getMaximumDouble(), dateRange.getFrom(), dateRange.getTo()));
198 Session session = SedDBSessionHolder.HOLDER.get(); 101 Session session = SedDBSessionHolder.HOLDER.get();
199 SQLQuery sqlQuery = session.createSQLQuery(SQL_BED_D50_SUBLAYER_MEASUREMENT) 102 SQLQuery sqlQuery = session.createSQLQuery(SQL_BED_D50_SUBLAYER_MEASUREMENT)
200 .addScalar("km", StandardBasicTypes.DOUBLE) 103 .addScalar("km", StandardBasicTypes.DOUBLE)
201 .addScalar("mindate", StandardBasicTypes.DATE) 104 .addScalar("datum", StandardBasicTypes.DATE)
202 .addScalar("maxdate", StandardBasicTypes.DATE) 105 .addScalar("tiefevon", StandardBasicTypes.DOUBLE)
203 .addScalar("cnt", StandardBasicTypes.INTEGER) 106 .addScalar("tiefebis", StandardBasicTypes.DOUBLE)
204 .addScalar("mindepth", StandardBasicTypes.DOUBLE) 107 .addScalar("d50", StandardBasicTypes.DOUBLE);
205 .addScalar("maxdepth", StandardBasicTypes.DOUBLE)
206 .addScalar("d50", StandardBasicTypes.DOUBLE);
207 String seddbRiver = river.nameForSeddb(); 108 String seddbRiver = river.nameForSeddb();
208 sqlQuery.setString("name", seddbRiver); 109 sqlQuery.setString("name", seddbRiver);
209 sqlQuery.setDouble("fromkm", kmRange.getMinimumDouble()); 110 sqlQuery.setDouble("fromkm", kmRange.getMinimumDouble());
210 sqlQuery.setDouble("tokm", kmRange.getMaximumDouble()); 111 sqlQuery.setDouble("tokm", kmRange.getMaximumDouble());
211 sqlQuery.setDate("fromdate", dateRange.getFrom()); 112 sqlQuery.setDate("fromdate", dateRange.getFrom());
212 sqlQuery.setDate("todate", dateRange.getTo()); 113 sqlQuery.setDate("todate", dateRange.getTo());
213 @SuppressWarnings("unchecked") 114 @SuppressWarnings("unchecked")
214 final List<Object[]> rows = sqlQuery.list(); 115 final List<Object[]> rows = sqlQuery.list();
215 final double[] kms = new double[rows.size()]; 116 final TDoubleArrayList kms = new TDoubleArrayList();
216 final double[] values = new double[rows.size()]; 117 final TDoubleArrayList values = new TDoubleArrayList();
217 D50Measurement measurement; 118 final TDoubleArrayList kmd50s = new TDoubleArrayList();
218 int i = -1; 119 for (int i = 0; i <= rows.size() - 1; i++) {
219 for (Object[] row : rows) { 120 kmd50s.add((double) rows.get(i)[4]);
220 measurement = new D50Measurement(row, SQL_BED_D50_SELECT_ALIAS); 121 if (((i == rows.size() - 1) || !Utils.epsilonEquals((double) rows.get(i)[0], (double) rows.get(i+1)[0], 0.0001))) {
221 i++; 122 int k = kmd50s.size() / 2;
222 kms[i] = measurement.getKm(); 123 values.add(((k + k < kmd50s.size()) ? kmd50s.get(k) : (kmd50s.get(k-1) + kmd50s.get(k)) / 2) / 1000);
223 values[i] = measurement.getD50(); 124 kms.add((double) rows.get(i)[0]);
224 log.debug(String.format("loadValues km %.3f d50(mm) %.1f count %d", kms[i], values[i], measurement.getCnt())); 125 log.debug(String.format("loadValues km %.3f d50(mm) %.1f count %d", kms.get(kms.size()-1), values.get(values.size()-1), kmd50s.size()));
126 kmd50s.clear();
127 }
225 } 128 }
226 try { 129 try {
227 interpolator = new LinearInterpolator().interpolate(kms, values); 130 interpolator = new LinearInterpolator().interpolate(kms.toNativeArray(), values.toNativeArray());
228 return true; 131 return true;
229 } catch (Exception e) { 132 } catch (Exception e) {
230 interpolator = null; 133 interpolator = null;
231 return false; 134 return false;
232 } 135 }
233 } 136 }
234
235 } 137 }
http://dive4elements.wald.intevation.org