SoilCalculator.java
/*
* This file is part of Indicators.
*
* Indicators is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Indicators is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Indicators. If not, see <https://www.gnu.org/licenses/>.
*/
package fr.inrae.agroclim.indicators.model.data.soil;
import java.text.SimpleDateFormat;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Calendar;
import java.util.Collection;
import java.util.Collections;
import java.util.Date;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Objects;
import java.util.Set;
import java.util.StringJoiner;
import java.util.TimeZone;
import fr.inrae.agroclim.indicators.model.TimeScale;
import fr.inrae.agroclim.indicators.model.data.DataLoadingListenerHandler;
import fr.inrae.agroclim.indicators.model.data.Variable;
import fr.inrae.agroclim.indicators.model.data.climate.ClimaticDailyData;
import fr.inrae.agroclim.indicators.model.data.phenology.Stage;
import fr.inrae.agroclim.indicators.resources.Messages;
import lombok.Getter;
import lombok.Setter;
import lombok.extern.log4j.Log4j2;
/**
* Calculator for Soil water balance.
*
* Last changed : $Date$
*
* @author $Author$
* @version $Revision$
*/
@Log4j2
public final class SoilCalculator extends DataLoadingListenerHandler implements HasSoilCalculatorParams, SoilLoader {
/**
* UID for serialization.
*/
private static final long serialVersionUID = 9107886175528703325L;
/**
* Four stages.
*/
private static final int FOUR = 4;
/**
* Converting meter to centimeter.
*/
private static final int M_TO_CM = 10;
/**
* Third stage.
*/
private static final int THREE = 3;
/**
* Climatic daily data used to compute soil water balance.
*/
@Setter
private List<ClimaticDailyData> climaticDailyData;
/**
* Coefficient cultural en situation non stressée durant la période initiale
* de la saison [0-1] (semis-levée selon notre hypothèse).
*/
@Setter
@Getter
private Double kcIni;
/**
* Coefficient cultural en situation non stressée durant la late season
* [0-1] (à partir de floraison selon nous).
*/
@Setter
@Getter
private Double kcLate;
/**
* Coefficient cultural en situation non stressée durant la middle season
* [0-1] (épi 1 cm - floraison selon notre hypothèse).
*/
@Setter
@Getter
private Double kcMid;
/**
* Proportion of easily available water in the soil.
*
* Average fraction of TAW, pour obtenir RAW.
*
* http://www.kimberly.uidaho.edu/water/fao56/fao56.pdf#189
*/
@Setter
private Double p;
/**
* Réserve utile du sol (mm).
*/
@Setter
private Double ru;
/**
* Soil depth (cm).
*/
@Setter
private Double soilDepth;
/**
* Date of the 4 stages by year needed to compute SWB and R.
*/
private List<Date> stages;
/**
* Date of the 4 stages by year needed to compute SWB and R.
*
* Computed from stages.
*/
private final Map<Integer, List<Date>> stagesByYear = new HashMap<>();
/**
* Teneur en eau du sol à la capacité au champ (% volumique).
*/
@Setter
@Getter
private Double swcFc;
/**
* Teneur en eau du sol maxi (% volumique).
*/
@Setter
@Getter
private Double swcMax;
/**
* Teneur en eau du sol au point de flétrissement permanent (% volumique).
*/
@Setter
@Getter
private Double swcWp;
/**
* Constructor.
*/
public SoilCalculator() {
// Do nothing
}
@Override
public SoilCalculator clone() {
final SoilCalculator clone = new SoilCalculator();
if (climaticDailyData != null) {
clone.climaticDailyData = new ArrayList<>();
Collections.copy(clone.climaticDailyData, climaticDailyData);
}
clone.kcIni = kcIni;
clone.kcLate = kcLate;
clone.kcMid = kcMid;
clone.p = p;
clone.ru = ru;
clone.soilDepth = soilDepth;
if (stages != null) {
clone.stages = new ArrayList<>();
Collections.copy(clone.stages, stages);
}
clone.swcFc = swcFc;
clone.swcMax = swcMax;
clone.swcWp = swcWp;
return clone;
}
/**
* Compute Kc.
*
* @param date date of climatic data
* @param previous date of previous stage
* @param current date of current stage
* @return Kc value
*/
private double computeKc(final Date date, final Date previous,
final Date current) {
final long timeInPhase = current.getTime() - previous.getTime();
final double daysInPhase = Math.ceil(timeInPhase
/ (1000 * 3600 * 24.));
final long timeFromStage3 = date.getTime() - previous.getTime();
final double daysFromStage3 = timeFromStage3
/ (1000 * 3600 * 24.);
return (kcLate - kcMid) / daysInPhase * daysFromStage3 + kcLate;
}
/**
* Compute SoilDailyData swc and waterReserve.
*
* @param soilDailyData daily data to update
* @param dailyData climatic data for the day
* @param stage1 stage 1
* @param stage2 stage 2
* @param stage3 stage 3
* @param stage4 stage 4
* @param firstDay is it the first day of climatic series
* @param rU Réserve utile du sol (mm).
* @param raw Readily available water en cm (au-delà, la plante commence à
* être stressée).
* @param rSat Réserve à saturation (mm).
* @param comptVeille compteur pour le ressuyage du jour précédent
* @param rVeille ressuyage du jour précédent
* @return new comptVeille
*/
@SuppressWarnings("checkstyle:ParameterNumber")
private int computeSoilDailyData(final SoilDailyData soilDailyData,
final ClimaticDailyData dailyData, final Date stage1,
final Date stage2, final Date stage3, final Date stage4,
final boolean firstDay, final double rU, final double raw,
final double rSat, final int comptVeille, final double rVeille) {
final Date date = dailyData.getDate();
// Crop evapotranspiration under non-standard conditions
double etCAdj;
final double et0 = dailyData.getEtp();
final double rain = dailyData.getRain();
final double kC;
// Calcul de Kc (X valeurs pour X phases phénologiques).
// SWBressuyage : lignes 151-159
if (date.before(stage1) || date.equals(stage1)) {
kC = kcIni;
} else if (date.before(stage2)) {
kC = computeKc(date, stage1, stage2);
} else if (date.before(stage3) || date.equals(stage3)) {
kC = kcMid;
} else if (date.before(stage4)) {
kC = computeKc(date, stage3, stage4);
} else {
kC = kcIni;
}
// Calcul de Ks en journalier à partir de la réserve de la veille
// (si le déficit est supérieure a la RAW, on estime que cela a un
// effet stress sur la transpiration de la plante).
int compt;
double kS;
double r;
double ressuyage;
if (firstDay) {
// cas particulier du premier jour de la première année
kS = rU / (rU - raw);
if (kS > 1) {
kS = 1;
}
//
etCAdj = kS * kC * et0;
// réserve en mm, bornée à la réserve à saturation
r = rU + rain - etCAdj;
if (r > rSat) {
r = rSat;
}
// compteur pour le ressuyage du jour suivant
if (r <= rU) {
compt = 0;
} else {
compt = 1;
}
} else {
// cas des autres jours
if (comptVeille == 0) {
// cas pas de ressuyage
if (rVeille < raw) {
kS = rVeille / (rU - raw);
} else {
kS = 1;
}
etCAdj = kS * kC * et0;
r = rVeille + rain - etCAdj;
} else {
// cas ressuyage
etCAdj = kC * et0;
switch (comptVeille) {
case 1:
// cas 1er jour de ressuyage
ressuyage = Math.min((rSat - rU) / 2, rVeille - rU);
// calcul de la réserve
r = rVeille - ressuyage - etCAdj + rain;
break;
case 2:
// cas 2e jour de ressuyage
r = rU + rain - etCAdj;
break;
default:
throw new RuntimeException("compt_veille != 1 && != 2");
}
}
// réserve en mm, bornée à la réserve à saturation
if (r > rSat) {
r = rSat;
}
// on calcule le compteur ressuyage pour le jour d'après
if (r <= rU) {
compt = 0;
} else {
switch (comptVeille) {
case 0:
compt = 1;
break;
case 1:
if (rain == 0) {
compt = 2;
} else {
compt = 1;
}
break;
case 2:
compt = 1;
break;
default:
throw new RuntimeException(
"compt_veille != 0 && != 1 && != 2");
}
}
}
// Calcul de SWC, dans le cas où SWC_WP est renseigné.
// On ne considère pas la profondeur maximale d'enracinement mais la
// profondeur de sol
if (swcWp != null) {
final double swc = r / (soilDepth / M_TO_CM) + swcWp;
soilDailyData.setSwc(swc);
}
soilDailyData.setWaterReserve(r);
return compt;
}
@Override
public Map<String, String> getConfigurationErrors() {
final Map<String, String> errors = new HashMap<>();
if (climaticDailyData == null) {
errors.put("soil.climaticDailyData", "error.evaluation.soil.climaticDailydata.missing");
} else if (climaticDailyData.isEmpty()) {
errors.put("soil.climaticDailyData", "error.evaluation.soil.climaticDailydata.empty");
}
// assert that the climatic data set is complete (1 data per day, all
// days per year, following years).
if (stages == null) {
errors.put("soil.stages", "error.evaluation.soil.stages.missing");
} else if (climaticDailyData != null && !climaticDailyData.isEmpty()) {
// check there are 4 stages for each year
setStagesByYear();
for (final Map.Entry<Integer, List<Date>> entry : stagesByYear.entrySet()) {
final Integer year = entry.getKey();
final List<Date> dates = entry.getValue();
if (dates.isEmpty()) {
LOGGER.error("The year {} does not have any stages.", year);
errors.put("soil.stages", "error.evaluation.soil.stages.4stagesperyear " + year + " 0");
} else if (dates.size() != Stage.FOUR) {
LOGGER.error("The year {} does not have {} stages {}, filling stages.", year, Stage.FOUR,
stagesByYear.get(year));
for (int nbOfStages = dates.size(); nbOfStages <= Stage.FOUR; nbOfStages++) {
stagesByYear.get(year).add(stagesByYear.get(year).get(nbOfStages - 1));
}
}
}
}
if (ru == null) {
if (swcFc == null) {
errors.put("soil.swcFc", "error.evaluation.soil.swcFc.missing");
}
if (swcWp == null) {
errors.put("soil.swcWp", "error.evaluation.soil.swcWp.missing");
}
if (soilDepth == null) {
errors.put("soil.soilDepth", "error.evaluation.soil.soilDepth.missing");
}
}
if (kcIni == null) {
errors.put("soil.kcIni", "error.evaluation.soil.kCini.missing");
}
if (kcMid == null) {
errors.put("soil.kcMid", "error.evaluation.soil.kCmid.missing");
}
if (kcLate == null) {
errors.put("soil.kcLate", "error.evaluation.soil.kClate.missing");
}
if (p == null) {
errors.put("soil.p", "error.evaluation.soil.p.missing");
}
if (swcMax == null) {
errors.put("soil.swcMax", "error.evaluation.soil.swcMax.missing");
}
if (errors.isEmpty()) {
return null;
}
return errors;
}
@Override
public Collection<String> getMissingVariables() {
throw new RuntimeException("Not implemented for soil!");
}
@Override
public Set<Variable> getProvidedVariables() {
final Set<Variable> variables = new HashSet<>();
variables.add(Variable.SOILWATERCONTENT);
variables.add(Variable.WATER_RESERVE);
return variables;
}
/**
* @return Readily available water en cm (au-delà, la plante commence à être
* stressée).
*/
private double getRAW() {
return p * getRu();
}
/**
* Calcul de la réserve à saturation.
*
* @return Réserve à saturation (mm).
*/
private double getRsat() {
if (swcMax == null) {
throw new IllegalStateException(
"swcMax must not be null to compute Rsat!");
}
if (swcWp == null) {
throw new IllegalStateException(
"swcWp must not be null to compute Rsat!");
}
if (soilDepth == null) {
throw new IllegalStateException(
"soilDepth must not be null to compute Rsat!");
}
return (swcMax - swcWp) * (soilDepth / M_TO_CM);
}
/**
* On calcule la RU sur la profondeur du sol à partir de soilDepth, SWC_FC
* et SWC_WP si elle n'a pas été renseignée.
*
* @return Réserve utile du sol (mm).
*/
private double getRu() {
if (ru == null) {
if (swcFc == null) {
throw new IllegalStateException(
"swcFc must not be null to compute RU!");
}
if (swcWp == null) {
throw new IllegalStateException(
"swcWp must not be null to compute RU!");
}
if (soilDepth == null) {
throw new IllegalStateException(
"soilDepth must not be null to compute RU!");
}
return (swcFc - swcWp) * (soilDepth / M_TO_CM);
}
return ru;
}
@Override
public Set<Variable> getVariables() {
final Set<Variable> variables = new HashSet<>();
variables.addAll(Arrays.asList(Variable.ETP, Variable.RAIN));
return variables;
}
/**
* Raise exception when there is no stage for a year.
*
* @param year the year for which there is no stage.
*/
private void handleMissingStages(final Integer year) {
final SimpleDateFormat sdf = new SimpleDateFormat("yyyy-MM-dd");
sdf.setTimeZone(TimeZone.getTimeZone("UTC"));
final StringJoiner sj = new StringJoiner(" ");
sj.add(Messages.format("error.soilcalculator.nostage", "" + year));
stages.stream() //
.filter(Objects::nonNull) //
.map(sdf::format) //
.forEach(sj::add);
throw new IllegalStateException(sj.toString());
}
@Override
public List<SoilDailyData> load() {
if (climaticDailyData == null) {
throw new IllegalStateException("No climatic daily data are set!");
}
if (climaticDailyData.isEmpty()) {
throw new IllegalStateException("Climatic daily data is empty!");
}
setStagesByYear();
final double rU = getRu();
final double raw = getRAW();
final double rSat = getRsat();
final List<SoilDailyData> data = new ArrayList<>();
double rVeille = 0.;
int comptVeille = 0;
boolean firstDay = true;
Integer year = climaticDailyData.get(0).getYear();
if (stagesByYear.get(year) == null) {
handleMissingStages(year);
}
Date stage1 = null;
Date stage2 = null;
Date stage3 = null;
Date stage4 = null;
if (stagesByYear.get(year).size() >= FOUR) {
stage1 = stagesByYear.get(year).get(0);
stage2 = stagesByYear.get(year).get(1);
stage3 = stagesByYear.get(year).get(2);
stage4 = stagesByYear.get(year).get(THREE);
}
// KC for the beginning of the year,
// this could be kcLate for the 2nd year
for (final ClimaticDailyData dailyData : climaticDailyData) {
if (year != dailyData.getYear()) {
year = dailyData.getYear();
if (stagesByYear.get(year) == null) {
handleMissingStages(year);
}
if (stagesByYear.get(year).size() >= FOUR) {
stage1 = stagesByYear.get(year).get(0);
stage2 = stagesByYear.get(year).get(1);
stage3 = stagesByYear.get(year).get(2);
stage4 = stagesByYear.get(year).get(THREE);
} else {
stage1 = null;
stage2 = null;
stage3 = null;
stage4 = null;
}
}
final SoilDailyData soilDailyData = new SoilDailyData();
soilDailyData.setYear(dailyData.getYear());
soilDailyData.setMonth(dailyData.getMonth());
soilDailyData.setDay(dailyData.getDay());
if (stage4 != null && stage3 != null && stage2 != null
&& stage1 != null) {
comptVeille = computeSoilDailyData(soilDailyData, dailyData,
stage1, stage2, stage3, stage4, firstDay, rU, raw, rSat,
comptVeille, rVeille);
rVeille = soilDailyData.getWaterReserve();
}
data.add(soilDailyData);
if (firstDay) {
firstDay = false;
}
}
return data;
}
/**
* @param fourStages
* the 4 stages (per year) needed to compute SWB and R
*/
public void setStages(final List<Date> fourStages) {
this.stages = fourStages;
stagesByYear.clear();
}
/**
* Divide all stages by year.
*/
private void setStagesByYear() {
if (stagesByYear.isEmpty()) {
final Calendar cal = Calendar.getInstance(TimeZone.getTimeZone("UTC"));
Integer year;
for (int i = 0; i < stages.size(); i++) {
final Date date = stages.get(i);
cal.setTime(date);
year = cal.get(Calendar.YEAR);
stagesByYear.computeIfAbsent(year, y -> new ArrayList<>());
stagesByYear.get(year).add(date);
Collections.sort(stagesByYear.get(year));
}
}
}
@Override
public void setTimeScale(final TimeScale timeScale) {
// do nothing
}
}