# HG changeset patch
# User Raimund Renkert <raimund.renkert@intevation.de>
# Date 1347632442 0
# Node ID 728ecd2afa200ec653ea87b213c4d07842f51259
# Parent  00aafe1fedd767cb42c544a3d3831fba1a89b1db
Implemented bed quality calculation in minfo module.



flys-artifacts/trunk@5474 c6561f87-3c4e-4783-a992-168aeb5c3f6f

diff -r 00aafe1fedd7 -r 728ecd2afa20 flys-artifacts/ChangeLog
--- a/flys-artifacts/ChangeLog	Fri Sep 14 14:14:46 2012 +0000
+++ b/flys-artifacts/ChangeLog	Fri Sep 14 14:20:42 2012 +0000
@@ -1,3 +1,15 @@
+2012-09-14  Raimund Renkert <raimund.renkert@intevation.de>
+
+	* src/main/java/de/intevation/flys/artifacts/model/minfo/QualityMeasurementFactory.java:
+	  Fixed SQL-statement.
+
+	* src/main/java/de/intevation/flys/artifacts/model/minfo/QualityMeasurements.java:
+	  Added getter for all kms.
+
+	* src/main/java/de/intevation/flys/artifacts/model/minfo/BedQualityCalculation.java:
+	  Implemented the bed quality calculation. There are still some fixes to do,
+	  e.g. extract a single result object for porosity and density.
+
 2012-09-14  Raimund Renkert <raimund.renkert@intevation.de>
 
 	* src/main/java/de/intevation/flys/artifacts/access/BedQualityAccess.java:
diff -r 00aafe1fedd7 -r 728ecd2afa20 flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/BedQualityCalculation.java
--- 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;
     }
 }
diff -r 00aafe1fedd7 -r 728ecd2afa20 flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/QualityMeasurementFactory.java
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/QualityMeasurementFactory.java	Fri Sep 14 14:14:46 2012 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/QualityMeasurementFactory.java	Fri Sep 14 14:20:42 2012 +0000
@@ -4,6 +4,7 @@
 import java.util.HashMap;
 import java.util.Map;
 
+import org.apache.log4j.Logger;
 import org.hibernate.SQLQuery;
 import org.hibernate.Session;
 import org.hibernate.transform.BasicTransformerAdapter;
@@ -14,11 +15,13 @@
 
 public class QualityMeasurementFactory {
 
+    private static Logger logger = Logger.getLogger(QualityMeasurementFactory.class);
+
     private static final String SQL_BED_MEASUREMENT =
         "SELECT st.km       as km," +
-        "       st.datum    as date," +
-        "       sp.tiefevon as depth1" +
-        "       sp.tiefebis as depth2" +
+        "       st.datum    as datum," +
+        "       sp.tiefevon as depth1," +
+        "       sp.tiefebis as depth2," +
         "       sa.d10      as d10," +
         "       sa.d16      as d16," +
         "       sa.d20      as d20," +
@@ -34,19 +37,21 @@
         "       sa.d90      as d90," +
         "       sa.dmin     as dmin," +
         "       sa.dmax     as dmax " +
-        "FROM sohltest st" +
-        "    JOIN station sn ON sn.stationid = st.stationid" +
-        "    JOIN gewaesser gw ON gw.gewaesserid = sn.gewaesserid" +
-        "    JOIN sohlprobe sp ON sp.sohltestid = st.sohltestid" +
+        "FROM sohltest st " +
+        "    JOIN station sn ON sn.stationid = st.stationid " +
+        "    JOIN gewaesser gw ON gw.gewaesserid = sn.gewaesserid " +
+        "    JOIN sohlprobe sp ON sp.sohltestid = st.sohltestid " +
         "    JOIN siebanalyse sa ON sa.sohlprobeid = sp.sohlprobeid " +
         "WHERE gw.name = :name AND " +
         "      st.km IS NOT NULL AND " +
-        "      st.km BETWEEN :from - 0.001 AND :to + 0.001 AND" +
+        "      sp.tiefevon IS NOT NULL AND " +
+        "      sp.tiefebis IS NOT NULL AND " +
+        "      st.km BETWEEN :from - 0.001 AND :to + 0.001 AND " +
         "      st.datum BETWEEN :start AND :end";
 
     private static final String SQL_BEDLOAD_MEASUREMENT =
         "SELECT m.km    as km," +
-        "       m.datum as date," +
+        "       m.datum as datum," +
         "       m.d10   as d10," +
         "       m.d16   as d16," +
         "       m.d20   as d20," +
@@ -83,14 +88,14 @@
             Map<String, Double> map = new HashMap<String, Double>();
             double km = 0;
             Date d = null;
-            double depth1 = 0;
-            double depth2 = 0;
+            double depth1 = Double.NaN;
+            double depth2 = Double.NaN;
             for (int i = 0; i < tuple.length; ++i) {
                 if (tuple[i] != null) {
                     if (aliases[i].equals("km")) {
                         km = ((Number) tuple[i]).doubleValue();
                     }
-                    else if (aliases[i].equals("date")) {
+                    else if (aliases[i].equals("datum")) {
                         d = (Date) tuple[i];
                     }
                     else if (aliases[i].equals("depth1")) {
@@ -100,7 +105,7 @@
                         depth2 = ((Number) tuple[i]).doubleValue();
                     }
                     else {
-                        map.put(aliases[i], (Double) tuple[i]);
+                        map.put(aliases[i], ((Double) tuple[i])/1000);
                     }
                 }
             }
@@ -122,7 +127,7 @@
     ) {
         SQLQuery query = session.createSQLQuery(statement)
             .addScalar("km", StandardBasicTypes.DOUBLE)
-            .addScalar("date", StandardBasicTypes.DATE)
+            .addScalar("datum", StandardBasicTypes.DATE)
             .addScalar("d10", StandardBasicTypes.DOUBLE)
             .addScalar("d16", StandardBasicTypes.DOUBLE)
             .addScalar("d20", StandardBasicTypes.DOUBLE)
@@ -168,7 +173,7 @@
             SQL_BED_MEASUREMENT);
         }
         finally {
-            session.close();
+            //session.close();
         }
     }
 
diff -r 00aafe1fedd7 -r 728ecd2afa20 flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/QualityMeasurements.java
--- a/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/QualityMeasurements.java	Fri Sep 14 14:14:46 2012 +0000
+++ b/flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/minfo/QualityMeasurements.java	Fri Sep 14 14:20:42 2012 +0000
@@ -3,8 +3,10 @@
 import java.util.LinkedList;
 import java.util.List;
 
+import org.apache.log4j.Logger;
+
 public class QualityMeasurements {
-
+    private static Logger logger = Logger.getLogger(QualityMeasurements.class);
     private List<QualityMeasurement> measurements;
 
     public QualityMeasurements() {
@@ -28,6 +30,16 @@
         return res;
     }
 
+    public List<Double> getKms() {
+        List<Double> result = new LinkedList<Double>();
+        for (QualityMeasurement qm : measurements) {
+            if (result.indexOf(qm.getKm()) < 0) {
+                result.add(qm.getKm());
+            }
+        }
+        return result;
+    }
+
     public void setMeasurements(List<QualityMeasurement> list) {
         this.measurements = list;
     }