Mercurial > dive4elements > river
view artifacts/src/main/java/org/dive4elements/river/artifacts/model/minfo/SedimentLoadCalculation.java @ 7027:21064459dc5d
Doc.
author | Felix Wolfsteller <felix.wolfsteller@intevation.de> |
---|---|
date | Mon, 16 Sep 2013 17:29:16 +0200 |
parents | a18a43764793 |
children | 7ee382b3f839 |
line wrap: on
line source
/* Copyright (C) 2011, 2012, 2013 by Bundesanstalt für Gewässerkunde * Software engineering by Intevation GmbH * * This file is Free Software under the GNU AGPL (>=v3) * and comes with ABSOLUTELY NO WARRANTY! Check out the * documentation coming with Dive4Elements River for details. */ package org.dive4elements.river.artifacts.model.minfo; import gnu.trove.TDoubleArrayList; import java.util.ArrayList; import java.util.List; import org.apache.log4j.Logger; import org.dive4elements.river.artifacts.access.SedimentLoadAccess; import org.dive4elements.river.artifacts.model.Calculation; import org.dive4elements.river.artifacts.model.CalculationResult; import org.dive4elements.river.artifacts.model.Range; /** Calculate sediment load. */ public class SedimentLoadCalculation extends Calculation { /** Private logger. */ private static final Logger logger = Logger .getLogger(SedimentLoadCalculation.class); protected String river; protected String yearEpoch; protected double kmUp; protected double kmLow; protected int[] period; /** Years of chosen epochs. */ protected int[][] epoch; protected String unit; public SedimentLoadCalculation() { } /** Returns CalculationResult with array of SedimentLoadResults. */ public CalculationResult calculate(SedimentLoadAccess access) { logger.info("SedimentLoadCalculation.calculate"); String river = access.getRiver(); String yearEpoch = access.getYearEpoch(); String unit = access.getUnit(); int[] period = null; int[][] epoch = null; double kmUp = access.getUpperKM(); double kmLow = access.getLowerKM(); if (yearEpoch.equals("year")) { period = access.getPeriod(); epoch = null; } else if (yearEpoch.equals("epoch") || yearEpoch.equals("off_epoch")) { epoch = access.getEpochs(); period = null; } else { addProblem("minfo.missing.year_epoch"); } if (river == null) { // TODO: i18n addProblem("minfo.missing.river"); } if (period == null && epoch == null) { addProblem("minfo.missing.time"); } if (!hasProblems()) { this.river = river; this.yearEpoch = yearEpoch; this.unit = unit; this.period = period; this.epoch = epoch; this.kmUp = kmUp; this.kmLow = kmLow; return internalCalculate(); } return new CalculationResult(); } /** Returns CalculationResult with array of SedimentLoadResults. */ private CalculationResult internalCalculate() { logger.debug("internalCalulate; mode:" + yearEpoch); if (yearEpoch.equals("year")) { List<SedimentLoadResult> results = new ArrayList<SedimentLoadResult>(); for (int i = period[0]; i <= period[1]; i++) { SedimentLoadResult res = calculateYear(i); results.add(res); } return new CalculationResult( results.toArray(new SedimentLoadResult[results.size()]), this); } else if (yearEpoch.equals("epoch")) { List<SedimentLoadResult> results = new ArrayList<SedimentLoadResult>(); for (int i = 0; i < epoch.length; i++) { SedimentLoadResult res = calculateEpoch(i); results.add(res); } return new CalculationResult( results.toArray(new SedimentLoadResult[results.size()]), this); } else if (yearEpoch.equals("off_epoch")) { List<SedimentLoadResult> results = new ArrayList<SedimentLoadResult>(); for (int i = 0; i < epoch.length; i++) { SedimentLoadResult res = calculateOffEpoch(i); results.add(res); } return new CalculationResult( results.toArray(new SedimentLoadResult[results.size()]), this); } else { logger.error("Unknown mode " + yearEpoch); } return null; } /** * @param[out] resLoad resulting SedimentLoad */ private void calculateEpochKm( List<SedimentLoad> epochLoads, SedimentLoad resLoad, double km ) { int cSum = 0; int fmSum = 0; int sSum = 0; int ssSum = 0; int ssbSum = 0; int sseSum = 0; for (SedimentLoad load : epochLoads) { SedimentLoadFraction f = load.getFraction(km); if (f.getCoarse() > 0d) { double c = resLoad.getFraction(km).getCoarse(); resLoad.setCoarse(km, c + f.getCoarse(), f.getCoarseRange()); cSum++; } if (f.getFineMiddle() > 0d) { double fm = resLoad.getFraction(km).getFineMiddle(); resLoad.setFineMiddle(km, fm + f.getFineMiddle(), f.getFineMiddleRange()); fmSum++; } if (f.getSand() > 0d) { double s = resLoad.getFraction(km).getSand(); resLoad.setSand(km, s + f.getSand(), f.getSandRange()); sSum++; } if (f.getSuspSand() > 0d) { double s = resLoad.getFraction(km).getSuspSand(); resLoad.setSuspSand(km, s + f.getSuspSand(), f.getSuspSandRange()); ssSum++; } if (f.getSuspSandBed() > 0d) { double s = resLoad.getFraction(km).getSuspSandBed(); resLoad.setSuspSandBed(km, s + f.getSuspSandBed(), f.getSuspSandBedRange()); ssbSum++; } if (f.getSuspSediment() > 0d) { double s = resLoad.getFraction(km).getSuspSediment(); resLoad.setSuspSediment(km, s + f.getSuspSediment(), f.getSuspSedimentRange()); sseSum++; } } SedimentLoadFraction fr = resLoad.getFraction(km); // Prevent divisions by zero, the fraction defaults to 0d. if (cSum != 0) { resLoad.setCoarse(km, fr.getCoarse()/cSum, fr.getCoarseRange()); } if (fmSum != 0) { resLoad.setFineMiddle(km, fr.getFineMiddle()/fmSum, fr.getFineMiddleRange()); } if (sSum != 0) { resLoad.setSand(km, fr.getSand()/sSum, fr.getSandRange()); } if (ssSum != 0) { resLoad.setSuspSand(km, fr.getSuspSand()/ssSum, fr.getSuspSandRange()); } if (ssbSum != 0) { resLoad.setSuspSandBed(km, fr.getSuspSandBed()/ssbSum, fr.getSuspSandBedRange()); } if (sseSum != 0) { resLoad.setSuspSediment(km, fr.getSuspSediment()/sseSum, fr.getSuspSedimentRange()); } } private SedimentLoadResult calculateEpoch(int i) { List<SedimentLoad> epochLoads = new ArrayList<SedimentLoad>(); for (int j = epoch[i][0]; j < epoch[i][1]; j++) { epochLoads.add(SedimentLoadFactory.getLoadWithData( this.river, this.yearEpoch, this.kmLow, this.kmUp, j, j)); } SedimentLoad resLoad = new SedimentLoad(); TDoubleArrayList kms = new TDoubleArrayList(); for (SedimentLoad load : epochLoads) { for (double km : load.getKms()) { if (!kms.contains(km)) { kms.add(km); } } } for (int j = 0; j < kms.size(); j++) { calculateEpochKm(epochLoads, resLoad, kms.get(j)); } resLoad.setDescription(""); resLoad.setEpoch(true); SedimentLoadResult result; SedimentLoad sl = calculateTotalLoad(resLoad, this.epoch[i][0]); if (this.unit.equals("m3_per_a")) { SedimentLoad slu = calculateUnit(sl, this.epoch[i][0]); result = new SedimentLoadResult( this.epoch[i][0], this.epoch[i][1], slu); } else { result = new SedimentLoadResult( this.epoch[i][0], this.epoch[i][1], sl); } return result; } /** * Calculate/Fetch values at off. epochs. * @param i index in epochs. */ private SedimentLoadResult calculateOffEpoch(int i) { SedimentLoad load = SedimentLoadFactory.getLoadWithData( this.river, this.yearEpoch, kmLow, kmUp, this.epoch[i][0], this.epoch[i][1]); SedimentLoadResult result; SedimentLoad sl = calculateTotalLoad(load, this.epoch[i][0]); if (unit.equals("m3_per_a")) { SedimentLoad slu = calculateUnit(sl, epoch[i][0]); result = new SedimentLoadResult( this.epoch[i][0], this.epoch[i][1], slu); } else { result = new SedimentLoadResult( this.epoch[i][0], this.epoch[i][1], sl); } return result; } /** * Fetch loads for a single year, calculate total and * return the result containing both. * @param y year, e.g. 1980 */ private SedimentLoadResult calculateYear(int y) { SedimentLoad load = SedimentLoadFactory.getLoadWithData( this.river, this.yearEpoch, this.kmLow, this.kmUp, y, y); SedimentLoadResult result; SedimentLoad sl = calculateTotalLoad(load, y); if (unit.equals("m3_per_a")) { SedimentLoad slu = calculateUnit(sl, y); result = new SedimentLoadResult(y, 0, slu); } else { result = new SedimentLoadResult(y, 0, sl); } return result; } /** Add up the loads of a year. */ private SedimentLoad calculateTotalLoad(SedimentLoad load, int year) { logger.debug("calculateTotalLoad"); boolean problemThisYear = false; if (!load.hasCoarse()) { addProblem("missing.fraction.coarse", Integer.toString(year)); problemThisYear = true; } if (!load.hasFineMiddle()) { addProblem("missing.fraction.fine_middle", Integer.toString(year)); problemThisYear = true; } if (!load.hasSand()) { addProblem("missing.fraction.sand", Integer.toString(year)); problemThisYear = true; } if (!load.hasSuspSand()) { addProblem("missing.fraction.susp_sand", Integer.toString(year)); problemThisYear = true; } if (!load.hasSuspSediment()) { addProblem("missing.fraction.susp_sediment", Integer.toString(year)); problemThisYear = true; } if (problemThisYear) { logger.warn("Some problem, not calculating total load."); return load; } return partialTotal(load); } /** * Set total values in load. * * Therefore, run over the sorted kms and find ranges where either all * or all Geschiebe or just the Schwebstoff fractions are set. * Merge these ranges and add (maybe new) respective fractions to * load. In the process, remember any 'unfished' ends from measurements * where the km-ranges did not completely match. * * @param load SedimentLoad to add total values (and ranges) to. * @return input param load, with total values set. */ private SedimentLoad partialTotal(SedimentLoad load) { // The load with balanced ranges, will be returned. SedimentLoad fairLoad = load; Range lastOtherRange = null; double lastOtherValue = 0d; Range lastSuspRange = null; double lastSuspValue = 0d; for (double km: load.getKms()) { // kms are already sorted! logger.debug ("Trying to add at km " + km); SedimentLoadFraction fraction = load.getFraction(km); if (fraction.isComplete()) { double total = fraction.getCoarse() + fraction.getFineMiddle() + fraction.getSand() + fraction.getSuspSand() + fraction.getSuspSediment(); // Easiest case. Add values up and set'em. if (fraction.getCoarseRange().equals( fraction.getSuspSedimentRange())) { lastOtherRange = null; lastSuspRange = null; fairLoad.setTotal(km, total, fraction.getCoarseRange()); } else { // Need to split a range. if (fraction.getCoarseRange().getEnd() < fraction.getSuspSedimentRange().getEnd()) { // Schwebstoff is longer. // Adjust and remember schwebstoffs range and value. lastSuspRange = (Range) fraction.getSuspSedimentRange().clone(); lastSuspRange.setStart(fraction.getCoarseRange().getEnd()); lastSuspValue = fraction.getSuspSediment(); lastOtherRange = null; fairLoad.setTotal(km, total, fraction.getCoarseRange()); } else { // Geschiebe is longer. // Adjust and remember other values. lastOtherRange = (Range) fraction.getCoarseRange().clone(); lastOtherRange.setStart(fraction.getSuspSedimentRange().getEnd()); lastOtherValue = (total - fraction.getSuspSediment()); lastSuspRange = null; fairLoad.setTotal(km, total, fraction.getSuspSedimentRange()); } } } else if (fraction.hasOnlySuspValues() && lastOtherRange != null) { // Split stuff. Range suspSedimentRange = fraction.getSuspSedimentRange(); // if intersects with last other range, cool! merge and add! if (lastOtherRange.contains(km)) { double maxStart = 0d; double minEnd = 0d; maxStart = Math.max(suspSedimentRange.getStart(), lastOtherRange.getStart()); minEnd = Math.min(suspSedimentRange.getEnd(), lastOtherRange.getEnd()); double total = lastOtherValue + fraction.getSuspSediment(); Range totalRange = new Range(maxStart, minEnd); if (suspSedimentRange.getEnd() > lastOtherRange.getEnd()) { lastSuspRange = (Range) suspSedimentRange.clone(); lastSuspRange.setStart(lastOtherRange.getEnd()); lastSuspValue = fraction.getSuspSediment(); lastOtherRange = null; } else { // Other is "longer". lastOtherRange.setStart(suspSedimentRange.getEnd()); lastSuspRange = null; } if (lastOtherRange != null && Math.abs(suspSedimentRange.getEnd() - lastOtherRange.getEnd()) < 0.1d) { lastOtherRange = null; lastSuspRange = null; } fairLoad.setTotal(km, total, totalRange); } else { lastSuspRange = suspSedimentRange; lastSuspValue = fraction.getSuspSediment(); lastOtherRange = null; } } else if (fraction.hasButSuspValues() && lastSuspRange != null) { // If intersects with last suspsed range, merge and add double total = fraction.getCoarse() + fraction.getFineMiddle() + fraction.getSand() + fraction.getSuspSand() + lastSuspValue; double maxStart = Math.max(fraction.getCoarseRange().getStart(), lastSuspRange.getStart()); if (lastSuspRange.contains(km)) { double minEnd = Math.min(fraction.getCoarseRange().getEnd(), lastSuspRange.getEnd()); Range totalRange = new Range(maxStart, minEnd); if (lastSuspRange.getEnd() > fraction.getCoarseRange().getEnd()) { // SuspSed longer. lastSuspRange.setStart(fraction.getCoarseRange().getEnd()); lastOtherRange = null; } else { // Other longer lastOtherRange = (Range) fraction.getCoarseRange().clone(); lastOtherRange.setStart(lastSuspRange.getEnd()); lastSuspRange = null; lastOtherValue = total - lastSuspValue; } if (lastSuspRange != null && lastOtherRange != null && Math.abs(lastSuspRange.getEnd() - lastOtherRange.getEnd()) < 0.1d) { lastOtherRange = null; lastSuspRange = null; } fairLoad.setTotal(km, total, totalRange); } else { // Ranges are disjoint. lastOtherRange = fraction.getCoarseRange(); lastOtherValue = total - fraction.getSuspSediment(); lastSuspRange = null; } } else { // Some values are missing or no intersection with former values. // Stay as we are. if (fraction.hasButSuspValues()) { double total = fraction.getCoarse() + fraction.getFineMiddle() + fraction.getSand() + fraction.getSuspSand(); lastOtherRange = fraction.getCoarseRange(); lastOtherValue = total; lastSuspRange = null; } else if (fraction.hasOnlySuspValues()) { lastSuspRange = fraction.getSuspSedimentRange(); lastSuspValue = fraction.getSuspSediment(); lastOtherRange = null; } } } return fairLoad; } /** * Transform values in load. * Background is to transform values measured in * t/a to m^3/a using the specific measured densities. * * @param load The load of which values should be transformed. * @param year The year at which to look at density (e.g. 2003). * * @return parameter load with transformed values. */ private SedimentLoad calculateUnit(SedimentLoad load, int year) { SedimentDensity density = SedimentDensityFactory.getSedimentDensity(river, kmLow, kmUp, year); for (double km: load.getKms()) { double dens = density.getDensity(km, year); SedimentLoadFraction fraction = load.getFraction(km); double coarse = fraction.getCoarse(); double fineMiddle = fraction.getFineMiddle(); double sand = fraction.getSand(); double suspSand = fraction.getSuspSand(); double bedSand = fraction.getSuspSandBed(); double sediment = fraction.getSuspSediment(); double total = fraction.getTotal(); double loadTotal = fraction.getLoadTotal(); load.setCoarse(km, (coarse * dens), fraction.getCoarseRange()); load.setFineMiddle(km, (fineMiddle * dens), fraction.getFineMiddleRange()); load.setSand(km, (sand * dens), fraction.getSandRange()); load.setSuspSand(km, (suspSand * dens), fraction.getSuspSandRange()); load.setSuspSandBed(km, (bedSand * dens), fraction.getSuspSandBedRange()); load.setSuspSediment(km, (sediment * dens), fraction.getSuspSedimentRange()); load.setTotal(km, (total * dens), fraction.getTotalRange()); load.setLoadTotal(km, (loadTotal * dens), fraction.getLoadTotalRange()); } return load; } } // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :