Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/BedQualityCalculation.java @ 3769:728ecd2afa20

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