Mercurial > dive4elements > river

comparison artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/BedQualityCalculation.java @ 5838:5aa05a7a34b7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rename modules to more fitting names.
author Sascha L. Teichmann <teichmann@intevation.de>
date Thu, 25 Apr 2013 15:23:37 +0200
parents flys-artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/BedQualityCalculation.java@bd047b71ab37
children 4897a58c8746
comparison
equal deleted inserted replaced
5837:d9901a08d0a6 5838:5aa05a7a34b7
1 package org.dive4elements.river.artifacts.model.minfo;
2
3 import gnu.trove.TDoubleArrayList;
4
5 import java.util.LinkedList;
6 import java.util.List;
7 import java.util.Map;
8
9 import org.apache.log4j.Logger;
10
11 import org.dive4elements.river.artifacts.access.BedQualityAccess;
12 import org.dive4elements.river.artifacts.model.Calculation;
13 import org.dive4elements.river.artifacts.model.CalculationResult;
14 import org.dive4elements.river.artifacts.model.DateRange;
15 import org.dive4elements.river.backend.SedDBSessionHolder;
16
17
18 public class BedQualityCalculation extends Calculation {
19
20 private static final Logger logger = Logger
21 .getLogger(BedQualityCalculation.class);
22
23 protected String river;
24 protected double from;
25 protected double to;
26 protected List<String> bedDiameter;
27 protected List<String> bedloadDiameter;
28 protected List<DateRange> ranges;
29
30 public BedQualityCalculation() {
31 }
32
33 public CalculationResult calculate(BedQualityAccess access) {
34 logger.info("BedQualityCalculation.calculate");
35
36 String river = access.getRiver();
37 Double from = access.getFrom();
38 Double to = access.getTo();
39 List<String> bedDiameter = access.getBedDiameter();
40 List<String> bedloadDiameter = access.getBedloadDiameter();
41 List<DateRange> ranges = access.getDateRanges();
42
43 if (river == null) {
44 // TODO: i18n
45 addProblem("minfo.missing.river");
46 }
47
48 if (from == null) {
49 // TODO: i18n
50 addProblem("minfo.missing.from");
51 }
52
53 if (to == null) {
54 // TODO: i18n
55 addProblem("minfo.missing.to");
56 }
57
58 if (ranges == null) {
59 // TODO: i18n
60 addProblem("minfo.missing.periods");
61 }
62
63 if (!hasProblems()) {
64 this.river = river;
65 this.from = from;
66 this.to = to;
67 this.ranges = ranges;
68 this.bedDiameter = bedDiameter;
69 this.bedloadDiameter = bedloadDiameter;
70
71 SedDBSessionHolder.acquire();
72 try {
73 return internalCalculate();
74 }
75 finally {
76 SedDBSessionHolder.release();
77 }
78 }
79
80 return new CalculationResult();
81 }
82
83 protected CalculationResult internalCalculate() {
84
85 List<BedQualityResult> results = new LinkedList<BedQualityResult>();
86 // Calculate for all time periods.
87 for (DateRange dr : ranges) {
88 QualityMeasurements loadMeasurements =
89 QualityMeasurementFactory.getBedloadMeasurements(
90 river,
91 from,
92 to,
93 dr.getFrom(),
94 dr.getTo());
95 QualityMeasurements bedMeasurements =
96 QualityMeasurementFactory.getBedMeasurements(
97 river,
98 from,
99 to,
100 dr.getFrom(),
101 dr.getTo());
102 BedQualityResult result = new BedQualityResult();
103 result.setDateRange(dr);
104 if (bedDiameter != null) {
105 result.add(calculateBedParameter(bedMeasurements, dr));
106 for (String bd : bedDiameter) {
107 BedDiameterResult bedResult =
108 calculateBed(bedMeasurements, bd, dr);
109
110 // Avoid adding empty result sets.
111 if (!bedResult.isEmpty()) {
112 result.add(bedResult);
113 }
114 }
115 }
116 if (bedloadDiameter != null) {
117 for (String bld : bedloadDiameter) {
118 BedloadDiameterResult loadResult =
119 calculateBedload(loadMeasurements, bld, dr);
120 result.add(loadResult);
121 }
122 }
123 results.add(result);
124 }
125
126 return new CalculationResult(
127 results.toArray(new BedQualityResult[results.size()]), this);
128 }
129
130 private BedParametersResult calculateBedParameter(
131 QualityMeasurements qm,
132 DateRange dr
133 ) {
134 List<Double> kms = qm.getKms();
135 QualityMeasurements capFiltered = filterCapMeasurements(qm);
136 QualityMeasurements subFiltered = filterSubMeasurements(qm);
137 TDoubleArrayList location = new TDoubleArrayList();
138 TDoubleArrayList porosityCap = new TDoubleArrayList();
139 TDoubleArrayList porositySub = new TDoubleArrayList();
140 TDoubleArrayList densityCap = new TDoubleArrayList();
141 TDoubleArrayList densitySub = new TDoubleArrayList();
142
143 for(double km : kms) {
144 double[] pCap = calculatePorosity(capFiltered, km);
145 double[] pSub = calculatePorosity(subFiltered, km);
146 double[] dCap = calculateDensity(capFiltered, pCap);
147 double[] dSub = calculateDensity(subFiltered, pSub);
148
149 double pCapRes = 0d;
150 double pSubRes = 0d;
151 double dCapRes = 0d;
152 double dSubRes = 0d;
153 for (int i = 0; i < pCap.length; i++) {
154 pCapRes += pCap[i];
155 dCapRes += dCap[i];
156 }
157 for (int i = 0; i < pSub.length; i++) {
158 pSubRes += pSub[i];
159 dSubRes += dSub[i];
160 }
161 location.add(km);
162 porosityCap.add((pCapRes / pCap.length) * 100 );
163 porositySub.add((pSubRes / pSub.length) * 100);
164 densityCap.add((dCapRes / dCap.length) / 1000);
165 densitySub.add((dSubRes / dSub.length) / 1000);
166
167 }
168
169 return new BedParametersResult(
170 location,
171 porosityCap,
172 porositySub,
173 densityCap,
174 densitySub);
175 }
176
177 protected BedDiameterResult calculateBed(
178 QualityMeasurements qm,
179 String diameter,
180 DateRange range
181 ) {
182 List<Double> kms = qm.getKms();
183 TDoubleArrayList location = new TDoubleArrayList();
184 TDoubleArrayList avDiameterCap = new TDoubleArrayList();
185 TDoubleArrayList avDiameterSub = new TDoubleArrayList();
186 for (double km : kms) {
187 //Filter cap and sub measurements.
188 QualityMeasurements capFiltered = filterCapMeasurements(qm);
189 QualityMeasurements subFiltered = filterSubMeasurements(qm);
190
191 List<QualityMeasurement> cm = capFiltered.getMeasurements(km);
192 List<QualityMeasurement> sm = subFiltered.getMeasurements(km);
193
194 double avCap = calculateAverage(cm, diameter);
195 double avSub = calculateAverage(sm, diameter);
196 location.add(km);
197 avDiameterCap.add(avCap * 1000);// bring to mm.
198 avDiameterSub.add(avSub * 1000);
199 }
200 return new BedDiameterResult(
201 diameter,
202 avDiameterCap,
203 avDiameterSub,
204 location);
205 }
206
207 private double[] calculateDensity(
208 QualityMeasurements capFiltered,
209 double[] porosity
210 ) {
211 double[] density = new double[porosity.length];
212 for (int i = 0; i < porosity.length; i++) {
213 density[i] = (1 - porosity[i]) * 2650;
214 }
215 return density;
216 }
217
218 private double[] calculatePorosity(
219 QualityMeasurements capFiltered,
220 double km
221 ) {
222 List<QualityMeasurement> list = capFiltered.getMeasurements(km);
223 double[] results = new double[list.size()];
224 int i = 0;
225 for (QualityMeasurement qm : list) {
226 double deviation = calculateDeviation(qm);
227 double p = calculateP(qm);
228 double porosity = 0.353 - 0.068 * deviation + 0.146 * p;
229 results[i] = porosity;
230 i++;
231 }
232
233 return results;
234 }
235
236 protected BedloadDiameterResult calculateBedload(
237 QualityMeasurements qm,
238 String diameter,
239 DateRange range
240 ) {
241 List<Double> kms = qm.getKms();
242 TDoubleArrayList location = new TDoubleArrayList();
243 TDoubleArrayList avDiameter = new TDoubleArrayList();
244 for (double km : kms) {
245 List<QualityMeasurement> measurements = qm.getMeasurements(km);
246 double mid = calculateAverage(measurements, diameter);
247 location.add(km);
248 avDiameter.add(mid);
249 }
250 return new BedloadDiameterResult(
251 diameter,
252 avDiameter,
253 location,
254 range);
255 }
256
257 protected double calculateAverage(
258 List<QualityMeasurement> list,
259 String diameter
260 ) {
261 double av = 0;
262 for (QualityMeasurement qm : list) {
263 av += qm.getDiameter(diameter);
264 }
265 return av/list.size();
266 }
267
268 protected QualityMeasurements filterCapMeasurements(
269 QualityMeasurements qms
270 ) {
271 List<QualityMeasurement> result = new LinkedList<QualityMeasurement>();
272 for (QualityMeasurement qm : qms.getMeasurements()) {
273 if (qm.getDepth1() == 0d && qm.getDepth2() <= 0.3) {
274 result.add(qm);
275 }
276 }
277 return new QualityMeasurements(result);
278 }
279
280 protected QualityMeasurements filterSubMeasurements(
281 QualityMeasurements qms
282 ) {
283 List<QualityMeasurement> result = new LinkedList<QualityMeasurement>();
284 for (QualityMeasurement qm : qms.getMeasurements()) {
285 if (qm.getDepth1() > 0d && qm.getDepth2() <= 0.5) {
286 result.add(qm);
287 }
288 }
289 return new QualityMeasurements(result);
290 }
291
292 public double calculateDeviation(QualityMeasurement qm) {
293 Map<String, Double> dm = qm.getAllDiameter();
294 double phiM = 0;
295 double[] phis = new double[dm.size()];
296 double[] ps = new double[dm.size()];
297 int i = 0;
298 for (String key : dm.keySet()) {
299 double d = dm.get(key);
300 double phi = -Math.log(d)/Math.log(2);
301 phis[i] = phi;
302 double p = calculateWeight(qm, key);
303 ps[i] = p;
304 phiM += phi * p;
305 i++;
306 }
307
308 double sig = 0d;
309 for (i = 0; i < dm.size(); i++) {
310 sig += ps[i] * Math.exp(phis[i] - phiM);
311 }
312 double deviation = Math.sqrt(sig);
313 return deviation;
314 }
315
316 protected double calculateP(QualityMeasurement qm) {
317 return calculateWeight(qm, "dmin");
318 }
319
320 public double calculateWeight(QualityMeasurement qm, String diameter) {
321 Map<String, Double> dm = qm.getAllDiameter();
322 double value = qm.getDiameter(diameter);
323
324 double sum = 0d;
325 for (Double d : dm.values()) {
326 sum =+ d.doubleValue();
327 }
328 double weight = sum/100*value;
329 return weight;
330 }
331 }
332 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org