teichmann@5863: /* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde
teichmann@5863:  * Software engineering by Intevation GmbH
teichmann@5863:  *
teichmann@5994:  * This file is Free Software under the GNU AGPL (>=v3)
teichmann@5863:  * and comes with ABSOLUTELY NO WARRANTY! Check out the
teichmann@5994:  * documentation coming with Dive4Elements River for details.
teichmann@5863:  */
teichmann@5863: 
teichmann@5831: package org.dive4elements.river.artifacts.model.minfo;
teichmann@5778: 
felix@6382: import gnu.trove.TDoubleArrayList;
felix@6382: 
teichmann@5831: import org.dive4elements.artifactdatabase.state.Facet;
teichmann@5778: 
teichmann@5831: import org.dive4elements.artifacts.Artifact;
teichmann@5831: import org.dive4elements.artifacts.CallContext;
teichmann@5778: 
teichmann@5867: import org.dive4elements.river.artifacts.D4EArtifact;
teichmann@5778: 
teichmann@5831: import org.dive4elements.river.artifacts.model.CalculationResult;
teichmann@5831: import org.dive4elements.river.artifacts.model.DataFacet;
teichmann@5831: import org.dive4elements.river.artifacts.model.FacetTypes;
teichmann@5831: 
teichmann@5831: import org.dive4elements.river.artifacts.states.DefaultState.ComputeType;
teichmann@5831: 
teichmann@5831: import org.dive4elements.river.model.MeasurementStation;
teichmann@5831: 
teichmann@5865: import org.dive4elements.river.utils.RiverUtils;
rrenkert@4372: 
felix@5645: import java.util.ArrayList;
felix@5653: import java.util.Collections;
felix@5645: import java.util.List;
felix@5653: import java.util.Map;
felix@6940: import java.util.TreeSet;
felix@5653: import java.util.TreeMap;
felix@5645: 
teichmann@5778: import org.apache.log4j.Logger;
felix@5645: 
rrenkert@4372: 
felix@5647: /** Facet to access various sediment loads. */
rrenkert@4372: public class SedimentLoadFacet
rrenkert@4372: extends DataFacet
rrenkert@4372: {
felix@5645:     /** Very own logger. */
rrenkert@4372:     private static Logger logger = Logger.getLogger(SedimentLoadFacet.class);
rrenkert@4372: 
felix@5647:     /** Used as tolerance value when fetching measurement stations. */
felix@5645:     private static double EPSILON = 1e-5;
felix@5645: 
felix@5647: 
rrenkert@4372:     public SedimentLoadFacet() {
rrenkert@4372:     }
rrenkert@4372: 
rrenkert@4372:     public SedimentLoadFacet(int idx, String name, String description,
rrenkert@4372:         ComputeType type, String stateId, String hash) {
rrenkert@4372:         super(idx, name, description, type, hash, stateId);
rrenkert@7894:         this.metaData.put("X", "chart.longitudinal.section.xaxis.label");
rrenkert@7894:         this.metaData.put("Y", "");
rrenkert@4372:     }
rrenkert@4372: 
felix@6679:     @Override
rrenkert@4372:     public Object getData(Artifact artifact, CallContext context) {
rrenkert@4372:         logger.debug("Get data for sediment load at index: " + index);
rrenkert@4372: 
teichmann@5867:         D4EArtifact flys = (D4EArtifact) artifact;
rrenkert@4372: 
rrenkert@4372:         CalculationResult res = (CalculationResult) flys.compute(context, hash,
rrenkert@4372:             stateId, type, false);
rrenkert@4372: 
rrenkert@4372:         Object[] data =
rrenkert@4372:             (SedimentLoadResult[]) res.getData(); // TODO CAST TO SPECIFIC CLASS
rrenkert@4372: 
teichmann@5865:         List<MeasurementStation> allStations = RiverUtils.getRiver(flys).getMeasurementStations();
felix@5645:         SedimentLoadResult result = data != null && data.length > index ? (SedimentLoadResult)data[index] : null;
felix@5645:         if (result == null) {
felix@5645:             return null;
felix@5645:         }
felix@5645: 
felix@6382:         // These complicated calculations were necessary because
felix@6382:         // SedimentLoad/Fraction did not contain the ranges of the given
felix@6382:         // values. Since this changed, the code is somewhat obsolete, but stable.
felix@6382:         // For an example of easier calculation, see the "total" part below.
felix@6382: 
felix@5653:         List<Double> sortedStarts = new ArrayList<Double>();
felix@5645:         // Filter stations according to type.
felix@5645:         List<MeasurementStation> stations = new ArrayList<MeasurementStation>();
felix@5645:         for (MeasurementStation station: allStations) {
felix@5645:             if (station.getRange() == null || station.getMeasurementType() == null) {
felix@5645:                 continue;
felix@5645:             }
felix@6412:             if (this.getName().contains("susp_sediment")
felix@6412:                 && station.getMeasurementType().equals("Schwebstoff")) {
felix@5645:                 stations.add(station);
felix@5653:                 sortedStarts.add(station.getStation());
felix@5653:             }
felix@6412:             else if (!this.getName().contains("susp_sediment")
felix@6412:                 && station.getMeasurementType().equals("Geschiebe")) {
felix@5645:                 stations.add(station);
felix@5653:                 sortedStarts.add(station.getStation());
felix@5653:             }
felix@5645:         }
felix@5653:         Collections.sort(sortedStarts);
felix@5645: 
felix@6382:         // Handle sediment load differently, as it respects already
felix@6382:         // the ranges that were added to SedimentLoad/Fraction.
felix@7501:         if (getName().equals(FacetTypes.SEDIMENT_LOAD_TA_TOTAL)
felix@7501:            ||getName().equals(FacetTypes.SEDIMENT_LOAD_M3A_TOTAL)) {
teichmann@8024:             SedimentLoadLSData load = result.getLoad();
felix@6382:             TDoubleArrayList xPos = new TDoubleArrayList();
felix@6382:             TDoubleArrayList yPos = new TDoubleArrayList();
felix@6382:             double lastX = -1d;
felix@6940:             for (double km: new TreeSet<Double>(load.getKms())) {
felix@6382:                 SedimentLoadFraction fraction = load.getFraction(km);
felix@6382:                 if (fraction.getTotal() != 0) {
felix@6382:                     if (Math.abs(lastX-km) >= EPSILON) {
felix@6382:                        xPos.add(Double.NaN);
felix@6382:                        yPos.add(Double.NaN);
felix@6382:                     }
felix@6382:                     xPos.add(km);
felix@6382:                     yPos.add(fraction.getTotal());
felix@6382:                     xPos.add(fraction.getTotalRange().getEnd());
felix@6382:                     yPos.add(fraction.getTotal());
felix@6382:                     lastX = fraction.getTotalRange().getEnd();
felix@6382:                 }
felix@6382:             }
felix@6382:             return new double[][] {xPos.toNativeArray(), yPos.toNativeArray()};
felix@6382:         }
felix@6382: 
felix@6382:         // Access data according to type (except total - see above).
felix@5648:         double[][] sd = getLoadData(result);
felix@5645: 
felix@5653:         // Sort by km.
felix@6340:         TreeMap<Double, Double> sortedKmLoad = new TreeMap<Double,Double>();
felix@5653: 
felix@5645:         double[] km   = sd[0];
felix@5645:         double[] load = sd[1];
felix@5645: 
felix@6342:         // Build map of km->load, but directly exclude the ones which do
felix@6342:         // not match against a measurements station ranges start.
felix@5653:         for (int i = 0 ; i < km.length; i++) {
felix@6342:             for (MeasurementStation station: stations) {
felix@6342:                 if (Math.abs(station.getStation() - km[i]) <= EPSILON) {
felix@6342:                     sortedKmLoad.put(km[i], load[i]);
felix@6342:                     continue;
felix@6342:                 }
felix@6342:             }
felix@5653:         }
felix@5653: 
felix@6341:         // [0] -> x, [1] -> y
felix@5645:         double[][] values = new double[2][];
felix@6342:         values[0] = new double[sortedKmLoad.size()*3];
felix@6342:         values[1] = new double[sortedKmLoad.size()*3];
felix@5645: 
felix@5645:         // Find station via its station (km).
felix@5653:         // TODO use a binarySearch instead of linear absdiff approach
felix@5653:         int i = 0;
felix@6340:         for (Map.Entry<Double, Double> kmLoad: sortedKmLoad.entrySet()) {
felix@5645:             boolean matchFound = false;
teichmann@7254:             for (int k = 0, S = stations.size(); k < S; k++) {
felix@6343:                 MeasurementStation station = stations.get(k);
felix@6343:                 if (Math.abs(station.getStation() - kmLoad.getKey()) < EPSILON) {
felix@6343:                     // Value has been taken at measurement station.
felix@6343:                     values[0][i*3]   = station.getRange().getA().doubleValue() + EPSILON;
felix@6343:                     values[1][i*3]   = kmLoad.getValue();
felix@6343:                     double endValue = 0d;
felix@6343:                     // Valid until next measurements stations begin of range,
felix@6343:                     // or end of current range if last value.
teichmann@7254:                     if (k+2 <= S) {
felix@6343:                         endValue = stations.get(k+1).getRange().getA().doubleValue();
felix@6343:                     }
felix@6343:                     else {
felix@6343:                         endValue = station.getRange().getB().doubleValue();
felix@6343:                     }
felix@6343:                     values[0][i*3+1] = endValue;
felix@6340:                     values[1][i*3+1] = kmLoad.getValue();
felix@6343:                     values[0][i*3+2] = endValue;
felix@6340:                     values[1][i*3+2] = kmLoad.getValue();
felix@5645:                     matchFound = true;
felix@5645:                 }
felix@5645:             }
felix@5653:             // Store points without match for later assessment.
felix@5645:             if (!matchFound) {
felix@6340:                 logger.warn("measurement without station ("+kmLoad.getKey()+")!");
felix@5653:             }
felix@5653:             i++;
felix@5653:         }
felix@5653: 
felix@5653:         for (int x = 0; x < values[0].length-1; x++) {
felix@5653:             // Introduce gaps where no data in measurement station.
felix@5653:             if (Math.abs(values[0][x+1] - values[0][x]) > 3*EPSILON
felix@5653:                     && values[1][x+1] != values[1][x]) {
felix@5653:                 values[0][x] = Double.NaN;
felix@5653:                 values[1][x] = Double.NaN;
felix@5645:             }
felix@5645:         }
felix@5653: 
felix@5645:         return values;
rrenkert@4372:     }
rrenkert@4372: 
felix@5645: 
felix@5648:     /** Get data according to type of facet. */
felix@5648:     private double[][] getLoadData(SedimentLoadResult result) {
felix@7501:         String name = getName();
felix@7501:         if (FacetTypes.IS.SEDIMENT_LOAD_SAND(name))
felix@5648:             return result.getSandData();
felix@7501:         else if (FacetTypes.IS.SEDIMENT_LOAD_COARSE(name))
felix@5648:             return result.getCoarseData();
felix@7501:         else if (FacetTypes.IS.SEDIMENT_LOAD_FINEMIDDLE(name))
felix@5648:             return result.getFineMiddleData();
felix@7501:         else if (FacetTypes.IS.SEDIMENT_LOAD_SUSP_SAND(name))
felix@5648:             return result.getSuspSandData();
felix@7501:         else if (FacetTypes.IS.SEDIMENT_LOAD_SUSP_SAND_BED(name))
felix@5648:             return result.getSuspSandBedData();
felix@7501:         else if (FacetTypes.IS.SEDIMENT_LOAD_SUSP_SEDIMENT(name))
felix@5648:             return result.getSuspSedimentData();
felix@7501:         else if (FacetTypes.IS.SEDIMENT_LOAD_TOTAL_LOAD(name))
felix@5648:             return result.getTotalLoadData();
felix@7501:         else if (FacetTypes.IS.SEDIMENT_LOAD_TOTAL(name))
felix@5648:             return result.getTotalData();
felix@5648:         else {
felix@7501:             logger.error("SedimentLoadFacet " + name + " cannot determine data type.");
felix@5648:             return null;
felix@5648:         }
felix@5648:     }
felix@5648: 
rrenkert@4372:     /** Copy deeply. */
rrenkert@4372:     @Override
rrenkert@4372:     public Facet deepCopy() {
rrenkert@4372:         SedimentLoadFacet copy = new SedimentLoadFacet();
rrenkert@4372:         copy.set(this);
rrenkert@4372:         copy.type = type;
rrenkert@4372:         copy.hash = hash;
rrenkert@4372:         copy.stateId = stateId;
rrenkert@4372:         return copy;
rrenkert@4372:     }
rrenkert@4372: }
rrenkert@4372: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :