Mercurial > dive4elements > river
diff flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/BedQualityCalculation.java @ 3769:728ecd2afa20
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 |
line wrap: on
line diff
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/BedQualityCalculation.java Fri Sep 14 14:14:46 2012 +0000 +++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/BedQualityCalculation.java Fri Sep 14 14:20:42 2012 +0000 @@ -1,7 +1,10 @@ package de.intevation.flys.artifacts.model.minfo; +import gnu.trove.TDoubleArrayList; + import java.util.LinkedList; import java.util.List; +import java.util.Map; import org.apache.log4j.Logger; @@ -20,6 +23,8 @@ protected String river; protected double from; protected double to; + protected List<String> bedDiameter; + protected List<String> bedloadDiameter; protected List<DateRange> ranges; public BedQualityCalculation() { @@ -31,6 +36,8 @@ String river = access.getRiver(); Double from = access.getFrom(); Double to = access.getTo(); + List<String> bedDiameter = access.getBedDiameter(); + List<String> bedloadDiameter = access.getBedloadDiameter(); List<DateRange> ranges = access.getDateRanges(); if (river == null) { @@ -58,6 +65,8 @@ this.from = from; this.to = to; this.ranges = ranges; + this.bedDiameter = bedDiameter; + this.bedloadDiameter = bedloadDiameter; SedDBSessionHolder.acquire(); try { @@ -76,28 +85,224 @@ List<BedQualityResult> results = new LinkedList<BedQualityResult>(); // Calculate for all time periods. for (DateRange dr : ranges) { - QualityMeasurements bedMeasurements = QualityMeasurementFactory - .getBedMeasurements(river, from, to, dr.getFrom(), dr.getTo()); - QualityMeasurements loadMeasurements = QualityMeasurementFactory - .getBedMeasurements(river, from, to, dr.getFrom(), dr.getTo()); + QualityMeasurements loadMeasurements = + QualityMeasurementFactory.getBedMeasurements( + river, + from, + to, + dr.getFrom(), + dr.getTo()); + QualityMeasurements bedMeasurements = + QualityMeasurementFactory.getBedMeasurements( + river, + from, + to, + dr.getFrom(), + dr.getTo()); + + for (String bd : bedDiameter) { + BedQualityResult bedResult = + calculateBed(bedMeasurements, bd, dr); - BedQualityResult bedResult = calculateBed(bedMeasurements); - BedQualityResult loadResult = calculateBedload(loadMeasurements); - results.add(bedResult); - results.add(loadResult); + //Avoid adding empty result sets. + if (!bedResult.isEmpty()) { + results.add(bedResult); + } + } + for (String bld : bedloadDiameter) { + BedQualityResult loadResult = + calculateBedload(loadMeasurements, bld, dr); + results.add(loadResult); + } } return new CalculationResult( results.toArray(new BedQualityResult[results.size()]), this); } - protected BedQualityResult calculateBed(QualityMeasurements qm) { - // TODO - return new BedQualityResult(); + protected BedQualityResult calculateBed( + QualityMeasurements qm, + String diameter, + DateRange range + ) { + List<Double> kms = qm.getKms(); + TDoubleArrayList location = new TDoubleArrayList(); + TDoubleArrayList avDiameterCap = new TDoubleArrayList(); + TDoubleArrayList avDiameterSub = new TDoubleArrayList(); + TDoubleArrayList porosityCap = new TDoubleArrayList(); + TDoubleArrayList porositySub = new TDoubleArrayList(); + TDoubleArrayList densityCap = new TDoubleArrayList(); + TDoubleArrayList densitySub = new TDoubleArrayList(); + for (double km : kms) { + //Filter cap and sub measurements. + QualityMeasurements capFiltered = filterCapMeasurements(qm); + QualityMeasurements subFiltered = filterSubMeasurements(qm); + + List<QualityMeasurement> cm = capFiltered.getMeasurements(km); + List<QualityMeasurement> sm = subFiltered.getMeasurements(km); + + double avCap = calculateAverage(cm, diameter); + double avSub = calculateAverage(sm, diameter); + double[] pCap = calculatePorosity(capFiltered, km, diameter); + double[] pSub = calculatePorosity(subFiltered, km, diameter); + double[] dCap = calculateDensity(capFiltered, pCap); + double[] dSub = calculateDensity(subFiltered, pSub); + + double pCapRes = 0d; + double pSubRes = 0d; + double dCapRes = 0d; + double dSubRes = 0d; + for (int i = 0; i < pCap.length; i++) { + pCapRes += pCap[i]; + dCapRes += dCap[i]; + } + for (int i = 0; i < pSub.length; i++) { + pSubRes += pSub[i]; + dSubRes += dSub[i]; + } + location.add(km); + avDiameterCap.add(avCap); + avDiameterSub.add(avSub); + porosityCap.add((pCapRes / pCap.length) * 100 ); + porositySub.add((pSubRes / pSub.length) * 100); + densityCap.add((dCapRes / dCap.length) / 1000); + densitySub.add((dSubRes / dSub.length) / 1000); + } + return new BedBedQualityResult( + diameter, + avDiameterCap, + avDiameterSub, + location, + range, + porosityCap, + porositySub, + densityCap, + densitySub); } - protected BedQualityResult calculateBedload(QualityMeasurements qm) { - // TODO - return new BedQualityResult(); + private double[] calculateDensity( + QualityMeasurements capFiltered, + double[] porosity + ) { + double[] density = new double[porosity.length]; + for (int i = 0; i < porosity.length; i++) { + density[i] = (1 - porosity[i]) * 2650; + } + return density; + } + + private double[] calculatePorosity( + QualityMeasurements capFiltered, + double km, + String diameter + ) { + List<QualityMeasurement> list = capFiltered.getMeasurements(km); + double[] results = new double[list.size()]; + int i = 0; + for (QualityMeasurement qm : list) { + double deviation = calculateDeviation(qm); + double p = calculateP(qm); + double porosity = 0.353 - 0.068 * deviation + 0.146 * p; + results[i] = porosity; + i++; + } + + return results; + } + + protected BedQualityResult calculateBedload( + QualityMeasurements qm, + String diameter, + DateRange range + ) { + List<Double> kms = qm.getKms(); + TDoubleArrayList location = new TDoubleArrayList(); + TDoubleArrayList avDiameter = new TDoubleArrayList(); + for (double km : kms) { + List<QualityMeasurement> measurements = qm.getMeasurements(km); + double mid = calculateAverage(measurements, diameter); + location.add(km); + avDiameter.add(mid); + } + return new BedLoadBedQualityResult( + diameter, + avDiameter, + location, + range); + } + + protected double calculateAverage( + List<QualityMeasurement> list, + String diameter + ) { + double av = 0; + for (QualityMeasurement qm : list) { + av += qm.getDiameter(diameter); + } + return av/list.size(); + } + + protected QualityMeasurements filterCapMeasurements( + QualityMeasurements qms + ) { + List<QualityMeasurement> result = new LinkedList<QualityMeasurement>(); + for (QualityMeasurement qm : qms.getMeasurements()) { + if (qm.getDepth1() == 0d && qm.getDepth2() <= 0.3) { + result.add(qm); + } + } + return new QualityMeasurements(result); + } + + protected QualityMeasurements filterSubMeasurements( + QualityMeasurements qms + ) { + List<QualityMeasurement> result = new LinkedList<QualityMeasurement>(); + for (QualityMeasurement qm : qms.getMeasurements()) { + if (qm.getDepth1() > 0d && qm.getDepth2() <= 0.5) { + result.add(qm); + } + } + return new QualityMeasurements(result); + } + + public double calculateDeviation(QualityMeasurement qm) { + Map<String, Double> dm = qm.getAllDiameter(); + double phiM = 0; + double[] phis = new double[dm.size()]; + double[] ps = new double[dm.size()]; + int i = 0; + for (String key : dm.keySet()) { + double d = dm.get(key); + double phi = -Math.log(d)/Math.log(2); + phis[i] = phi; + double p = calculateWeight(qm, key); + ps[i] = p; + phiM += phi * p; + i++; + } + + double sig = 0d; + for (i = 0; i < dm.size(); i++) { + sig += ps[i] * Math.exp(phis[i] - phiM); + } + double deviation = Math.sqrt(sig); + return deviation; + } + + protected double calculateP(QualityMeasurement qm) { + return calculateWeight(qm, "dmin"); + } + + public double calculateWeight(QualityMeasurement qm, String diameter) { + Map<String, Double> dm = qm.getAllDiameter(); + double value = qm.getDiameter(diameter); + + double sum = 0d; + for (Double d : dm.values()) { + sum =+ d.doubleValue(); + } + double weight = sum/100*value; + return weight; } }