PhenologyCalculator.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.phenology;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
import java.util.Collections;
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 fr.inrae.agroclim.indicators.model.TimeScale;
import fr.inrae.agroclim.indicators.model.data.DataLoadingListenerHandler;
import fr.inrae.agroclim.indicators.model.data.ResourcesLoader;
import fr.inrae.agroclim.indicators.model.data.Variable;
import fr.inrae.agroclim.indicators.model.data.climate.ClimaticDailyData;
import jakarta.xml.bind.annotation.XmlAccessType;
import jakarta.xml.bind.annotation.XmlAccessorType;
import jakarta.xml.bind.annotation.XmlElement;
import jakarta.xml.bind.annotation.XmlTransient;
import jakarta.xml.bind.annotation.XmlType;
import lombok.EqualsAndHashCode;
import lombok.Getter;
import lombok.NonNull;
import lombok.Setter;
import lombok.extern.log4j.Log4j2;
/**
* Calculator to provide phenological stages.
*
* Definitions : Vernalization is the induction of a plant's flowering process
* by exposure to the prolonged cold of winter, or by an artificial equivalent.
* After vernalization, plants have acquired the ability to flower, but they may
* require additional seasonal cues or weeks of growth before they will actually
* flower. Vernalization is sometimes used to refer to herbal (non-woody) plants
* requiring a cold dormancy to produce new shoots and leaves[1] but this usage
* is discouraged
*
* Phenology is related to the harvesting year.
*
* Soil water balance needs 4-stages phenological model.
*
* Curvilinear model uses Tmin, Topt, Tmax. Linear model uses Tbase.
*
* @author Olivier Maury
*/
@XmlAccessorType(XmlAccessType.FIELD)
@XmlType(propOrder = {"crop", "variety", "allYears", "chuineA", "chuineB",
"chuineC", "baseTemperature", "isPhotoperiodToCompute",
"isVernalizationToCompute", "maxTemperature", "maxTemperatureDeb",
"maxTemperatureVer", "minNbOfVernalizationDays", "minTemperature",
"minTemperatureDeb", "minTemperatureVer", "model", "nbOfStages",
"nbOfVernalizationDays", "nbOfYears", "optimalTemperature",
"optimalTemperatureDeb", "optimalTemperatureVer",
"optimalVernalizationAmplitude", "optimalVernalizationTemperature",
"phase1Sum", "phase2Sum", "phase3Sum", "phase4Sum", "phase5Sum",
"phase6Sum", "photoperiodBase", "photoperiodSaturating",
"photoperiodSensitivity", "s3ToS4", "siteLatitude", "sowingDate"})
@EqualsAndHashCode(
callSuper = false,
of = {"crop", "variety", "allYears", "chuineA", "chuineB",
"chuineC", "baseTemperature", "isPhotoperiodToCompute",
"isVernalizationToCompute", "maxTemperature", "maxTemperatureDeb",
"maxTemperatureVer", "minNbOfVernalizationDays", "minTemperature",
"minTemperatureDeb", "minTemperatureVer", "model",
"nbOfVernalizationDays", "nbOfYears", "optimalTemperature",
"optimalTemperatureDeb", "optimalTemperatureVer",
"optimalVernalizationAmplitude", "optimalVernalizationTemperature",
"photoperiodBase", "photoperiodSaturating",
"photoperiodSensitivity", "s3ToS4", "siteLatitude", "sowingDate"}
)
@Log4j2
public final class PhenologyCalculator extends DataLoadingListenerHandler
implements ResourcesLoader<List<AnnualStageData>> {
/**
* UID for serialization.
*/
private static final long serialVersionUID = -5403460610691492829L;
/**
* All years have the same number of days.
*/
private static final int DAYS_IN_YEAR = 365;
/**
* Default number of days between S3 and S4 for curve_grapevine.
*/
private static final int DEFAULT_S3_TO_S4_CURVE_GRAPEVINE = 35;
/**
* Compute the cumulative sum of JVI.
*
* @param optiVernAmp Thermal half-amplitude around the optimal temperature
* for vernalization (°C) (ver_ampfroid).
* @param optiVernTemp Optimal temperature for vernalization (°C)
* (Ver_tfroid).
* @param tmeans average temperatures
* @param startDate start of phase (day of year)
* @param initialCumjvi value of vernalization effect before phase
* @return list for all days of cumulative sum of JVI
*/
public static double[] cumjvi(final double optiVernAmp,
final double optiVernTemp, @NonNull final List<Double> tmeans,
final int startDate, final double initialCumjvi) {
double cumjvi = initialCumjvi;
final double[] results = new double[tmeans.size() - startDate + 1];
for (int i = startDate; i < tmeans.size(); i++) {
final Double tmean = tmeans.get(i - 1);
final double d = 1
- Math.pow((optiVernTemp - tmean) / optiVernAmp, 2);
if (d > 0) {
cumjvi += d;
}
results[i - startDate] = cumjvi;
}
return results;
}
/**
* Compute the photoperiod number of hours according to a given latitude on
* a specific julian day.
*
* @param day
* julian day
* @param latitude
* given latitude
* @return number of hours
*/
public static double photoperiod(final int day, final double latitude) {
final int hoursInDay = 24;
final int daysInYear = 365;
final double alat = latitude / 57.296;
int doy;
if (day > daysInYear) {
doy = day - daysInYear;
} else {
doy = day;
}
final double theta1 = 2 * Math.PI * (doy - 80) / daysInYear;
final double theta2 = 0.034 * (Math.sin(2 * Math.PI * doy / daysInYear)
- Math.sin(2 * Math.PI * 80 / daysInYear));
final double theta = theta1 + theta2;
final double dec = Math.asin(0.3978 * Math.sin(theta));
final double d = -1.0 * 0.10453 / (Math.cos(alat) * Math.cos(dec));
double p = d - Math.tan(alat) * Math.tan(dec);
p = Math.acos(p);
final double photoperiod = hoursInDay * p / Math.PI;
if (photoperiod > hoursInDay) {
return hoursInDay;
} else {
return photoperiod;
}
}
/**
* Compute the RFPI (Effect of the photoperiod).
*
* RFPI = "effet de la photopériode, coefficient compris entre 0 et 1 qui
* tient compte des paramètres plante et de la latitude." ⮕ tiré du code R
*
* @param day
* julian day
* @param latitude
* given latitude
* @param sensitivity
* photoperiod sensitivity (1=insensitive)
* @param saturating
* photoperiod saturating (hours)
* @param base
* photoperiod base (hours)
* @return RFPI
*/
public static double rfpi(final int day, final double latitude,
final double sensitivity, final double saturating,
final double base) {
final double photoperiod = photoperiod(day, latitude);
double cRFPI = (1.0 - sensitivity) / (saturating - base)
* (photoperiod - saturating) + 1;
cRFPI = Math.max(Math.min(cRFPI, 1), sensitivity);
return cRFPI;
}
/**
* Effect of vernalization for a day.
*
* @param nbOfVernalizationDays number of vernalization days
* @param minNbOfVernalizationDays minimal number of vernalization days
* @param cumjvi cumulative sum of JVI for the day
* @return effect of vernalization
*/
public static double rfvi(final int nbOfVernalizationDays,
final int minNbOfVernalizationDays, final double cumjvi) {
if (cumjvi < nbOfVernalizationDays) {
double rfvi = (cumjvi - minNbOfVernalizationDays)
/ (nbOfVernalizationDays - minNbOfVernalizationDays);
if (rfvi < 0) {
rfvi = 0;
}
if (rfvi > 1) {
rfvi = 1;
}
return rfvi;
}
return 1.;
}
/**
* Compute stages for all years, whatever is nbYears.
*
* True for SoilCalculator.
*/
@Getter
@Setter
private boolean allYears = false;
/**
* Minimal temperature for development (°C) used in the curvilinear model.
*/
@Getter
@Setter
private Double minTemperature;
/**
* Minimal temperature for development (°C) for DEB stage, used in the
* curvilinear grapevine model.
*/
@Getter
@Setter
private Double minTemperatureDeb;
/**
* Minimal temperature for development (°C) for VER stage, used in the
* curvilinear grapevine model.
*/
@Getter
@Setter
private Double minTemperatureVer;
/**
* Maximal temperature for development (°C) used in the curvilinear model.
*/
@Getter
@Setter
private Double maxTemperature;
/**
* Maximal temperatures for development (°C) for DEB stage, used in the
* curvilinear grapevine model.
*/
@Getter
@Setter
private Double maxTemperatureDeb;
/**
* Maximal temperatures for development (°C) for VER stage, used in the
* curvilinear grapevine model.
*/
@Getter
@Setter
private Double maxTemperatureVer;
/**
* Optimal temperature for development (°C) used in the curvilinear model.
*/
@Getter
@Setter
private Double optimalTemperature;
/**
* Optimal temperatures for development (°C) for DEB stage, used in the
* curvilinear grapevine model.
*/
@Getter
@Setter
private Double optimalTemperatureDeb;
/**
* Optimal temperatures for development (°C) for VER stage, used in the
* curvilinear grapevine model.
*/
@Getter
@Setter
private Double optimalTemperatureVer;
/**
* Base temperature for development (°C) used in linear model.
*/
@Getter
@Setter
private Double baseTemperature;
/**
* Parameter for Chuine function, used in curve_grapevine* models.
*/
@Getter
@Setter
private Double chuineA;
/**
* Parameter for Chuine function, used in curve_grapevine* models.
*/
@Getter
@Setter
private Double chuineB;
/**
* Parameter for Chuine function, used in curve_grapevine* models.
*/
@Getter
@Setter
private Double chuineC;
/**
* Climatic daily data used to compute phenology.
*/
@XmlTransient
@Getter
@Setter
private List<ClimaticDailyData> climaticDailyData;
/**
* Given crop.
*/
@Getter
@Setter
private String crop;
/**
* Compute photoperiod in phenology model ?
*
* Photoperiod = day duration.
*/
@Getter
@Setter
private boolean isPhotoperiodToCompute;
/**
* Compute vernalization in phenology model ?
*
* Vernalization = need cold.
*/
@Getter
@Setter
private boolean isVernalizationToCompute;
/**
* Minimum number of vernalization days.
*/
@Getter
@Setter
private int minNbOfVernalizationDays;
/**
* Given model.
*/
@Getter
@Setter
private PhenologicalModelType model;
/**
* Number of vernalization days.
*/
@Getter
@Setter
private int nbOfVernalizationDays;
/**
* Type of crop.
*
* nbOfYears = 1 (sowing and harvest the same year)
*
* nbOfYears = 2 (harvest the year following the year of sowing)
*/
@Getter
@Setter
private int nbOfYears;
/**
* Thermal half-amplitude around the optimal temperature for vernalization
* (°C).
*/
@Getter
@Setter
private Double optimalVernalizationAmplitude;
/**
* Optimal temperature for vernalization (°C).
*/
@Getter
@Setter
private Double optimalVernalizationTemperature;
/**
* Sum of effective temperature from
*
* - sowing to emergence (°C), aka ST_Phase1_4 OR ST_Phase1_6.
*
* - emergence to maximum acceleration of leaf growth (°C), aka ST_Phase2_4,
*
* - emergence to début tallage (°C), aka ST_Phase2_6.
*
* - maximum acceleration of leaf growth to flowering (°C), aka ST_Phase3_4,
*
* - début tallage to epi 1 cm (°C), aka ST_Phase3_6
*
* - flowering to physiological maturity (°C), aka ST_Phase4_4,
*
* - epi 1 cm to meiose (°C), aka ST_Phase4_6.
*
* - meiose to floraison (°C), aka ST_Phase5_6,
*
* - floraison to physiological maturity (°C), aka ST_Phase6_6.
*/
@XmlTransient
private Double[] phaseSums;
/**
* Base photoperiod (number of hours).
*/
@Getter
@Setter
private double photoperiodBase;
/**
* Saturation photoperiod (number of hours).
*/
@Getter
@Setter
private double photoperiodSaturating;
/**
* Sensibility to the photoperiod (from 0 to 1 with 1 = insensitive).
*/
@Getter
@Setter
private double photoperiodSensitivity;
/**
* Number of days between S3 and S4 for curve_grapevine.
*
* By default: 35 days.
*/
@Getter
@Setter
private int s3ToS4 = DEFAULT_S3_TO_S4_CURVE_GRAPEVINE;
/**
* Site latitude (°N) to compute photoperiod.
*/
@Getter
@Setter
private double siteLatitude;
/**
* Sowing date (julian day).
*/
@Getter
@Setter
private int sowingDate;
/**
* Given crop variety.
*/
@Getter
@Setter
private String variety;
/**
* Constructor.
*/
public PhenologyCalculator() {
// Do nothing
}
@Override
public PhenologyCalculator clone() {
final PhenologyCalculator clone = new PhenologyCalculator();
clone.baseTemperature = baseTemperature;
if (climaticDailyData != null) {
clone.climaticDailyData = new ArrayList<>();
clone.climaticDailyData.addAll(climaticDailyData);
}
clone.crop = crop;
clone.isPhotoperiodToCompute = isPhotoperiodToCompute;
clone.isVernalizationToCompute = isVernalizationToCompute;
clone.maxTemperature = maxTemperature;
clone.minNbOfVernalizationDays = minNbOfVernalizationDays;
clone.minTemperature = minTemperature;
clone.model = model;
clone.nbOfVernalizationDays = nbOfVernalizationDays;
clone.nbOfYears = nbOfYears;
clone.optimalTemperature = optimalTemperature;
clone.optimalVernalizationAmplitude = optimalVernalizationAmplitude;
clone.optimalVernalizationTemperature = optimalVernalizationTemperature;
clone.phaseSums = Arrays.copyOf(phaseSums, phaseSums.length);
clone.photoperiodBase = photoperiodBase;
clone.photoperiodSaturating = photoperiodSaturating;
clone.photoperiodSensitivity = photoperiodSensitivity;
clone.sowingDate = sowingDate;
clone.variety = variety;
return clone;
}
/**
* @return year ⮕ Tmean
*/
private Map<Integer, List<Double>> getClimaticDataByYear() {
if (climaticDailyData == null) {
throw new RuntimeException("climaticDailyData is not defined!");
}
if (climaticDailyData.isEmpty()) {
throw new RuntimeException("climaticDailyData is empty!");
}
final Map<Integer, List<Double>> dataByYear = new HashMap<>();
for (final ClimaticDailyData climaticData : climaticDailyData) {
final int year = climaticData.getYear();
dataByYear.computeIfAbsent(year, y -> new ArrayList<>());
if (climaticData.getTmean() == null) {
throw new RuntimeException("Null value for TMEAN at "
+ climaticData.getYear() + "/"
+ climaticData.getMonth() + "/"
+ climaticData.getDay());
}
dataByYear.get(year).add(climaticData.getTmean());
}
return dataByYear;
}
@Override
public Map<String, String> getConfigurationErrors() {
final Map<String, String> errors = new HashMap<>();
if (climaticDailyData == null) {
errors.put("phenology.climaticDailyData", "error.evaluation.phenology.climaticDailydata.missing");
} else if (climaticDailyData.isEmpty()) {
errors.put("phenology.climaticDailyData", "error.evaluation.phenology.climaticDailydata.empty");
}
if (errors.isEmpty()) {
return null;
}
return errors;
}
@Override
public Collection<String> getMissingVariables() {
throw new UnsupportedOperationException("Not implemented for phenology!");
}
/**
* @return number of stages to compute.
*/
public int getNbOfStages() {
return phaseSums.length;
}
/**
* @return ST_Phase1_4, ST_Phase1_6
*/
public Double getPhase1Sum() {
final int phase = 1;
return getPhaseSumOrNull(phase);
}
/**
* @return ST_Phase2_4, ST_Phase2_6
*/
public Double getPhase2Sum() {
final int phase = 2;
return getPhaseSumOrNull(phase);
}
/**
* @return ST_Phase3_4, ST_Phase3_6
*/
public Double getPhase3Sum() {
final int phase = 3;
return getPhaseSumOrNull(phase);
}
/**
* @return ST_Phase4_4, ST_Phase4_6
*/
public Double getPhase4Sum() {
final int phase = 4;
return getPhaseSumOrNull(phase);
}
/**
* @return ST_Phase5_6
*/
public Double getPhase5Sum() {
final int phase = 5;
return getPhaseSumOrNull(phase);
}
/**
* @return ST_Phase6_6
*/
public Double getPhase6Sum() {
final int phase = 6;
return getPhaseSumOrNull(phase);
}
/**
* @param phase
* phase number
* @return Sum of effective temperature for the phase
*/
public double getPhaseSum(final int phase) {
if (phase < 1 || phase > phaseSums.length) {
throw new IllegalArgumentException("Phase must be between 1 and "
+ phaseSums.length);
}
return phaseSums[phase - 1];
}
/**
* @param phase phase number
* @return Sum of effective temperature for the phase
*/
private Double getPhaseSumOrNull(final int phase) {
if (phase < 1 || phase > phaseSums.length) {
return null;
}
return phaseSums[phase - 1];
}
@Override
public Set<Variable> getVariables() {
final Set<Variable> variables = new HashSet<>();
variables.add(Variable.TMEAN);
return variables;
}
/**
* This is a translation of the R code from the R "indicateurs" project.
*
* General idea: 1/ get average temperatures for 1 or 2 years according to
* crop type, 2/ Compute Phase 1 to N: sum of difference between base
* temperature (* photoperiod effect * vernalization effect) and average
* until threshold to get end date of phase
*
* Phases for linear 4 stages: 1 : sowing - emergence, 2 : emergence - LAI
* inflexion, 3 : inflexion - flowering, 4 : flowering - maturity
*
* The original algorithms are : pheno_curve_model_4stades.r,
* pheno_linear_model_4stades.r, pheno_linear_model_5stades.r,
* pheno_linear_model_6stades.r.
*/
@Override
public List<AnnualStageData> load() {
LOGGER.traceEntry();
if (model == null) {
throw new RuntimeException("model is not defined!");
}
if (phaseSums == null) {
throw new RuntimeException("phaseSums is not defined!");
}
for (int i = 0; i < phaseSums.length; i++) {
if (phaseSums[i] == null) {
final StringJoiner sj = new StringJoiner("\n");
sj.add("Missing PhaseSums[" + i + "] for variety " + variety);
for (int j = 0; j < phaseSums.length; j++) {
sj.add(String.format(
"PhaseSums[%d]=%s", j, phaseSums[j]));
}
LOGGER.error(sj.toString());
throw new RuntimeException("phaseSum is not defined for stage "
+ (i + 1));
}
}
final List<AnnualStageData> phenologicalData = new ArrayList<>();
Map<Integer, List<Double>> climaticDataByYear;
climaticDataByYear = getClimaticDataByYear();
final List<Integer> years;
years = new ArrayList<>(climaticDataByYear.keySet());
Collections.sort(years);
final Integer lastYear = years.get(years.size() - 1);
final String stageNamePrefix = "s";
for (final Integer year : years) {
final boolean isLastYear = Objects.equals(year, lastYear);
if (!allYears && nbOfYears != 1 && isLastYear) {
LOGGER.trace("Ignoring last year {}", year);
continue;
}
// The calculation is made for each climatic year
final List<Double> dataForCurrentYear;
dataForCurrentYear = climaticDataByYear.get(year);
final List<Double> tmeans = new ArrayList<>(dataForCurrentYear);
if (nbOfYears == 2 && !(allYears && isLastYear)) {
final List<Double> dataForNextYear;
dataForNextYear = climaticDataByYear.get(year + 1);
if (dataForNextYear == null) {
throw new RuntimeException("No climatic data for year "
+ (year + 1));
}
// If the crop grows between 2 years
tmeans.addAll(dataForNextYear);
}
final List<Stage> stages = new ArrayList<>();
// Stage 0
final Stage s0 = new Stage(stageNamePrefix + "0", sowingDate);
stages.add(s0);
// Stages 1 - nbOfStages
int phaseStart = sowingDate;
double cumjvi = 0.;
for (int stage = 1; stage <= phaseSums.length; stage++) {
final String stageName = stageNamePrefix + stage;
PhaseEnd phaseEnd;
phaseEnd = model.getEndPhase(this, tmeans, phaseStart, cumjvi, stage);
cumjvi = phaseEnd.getCumjvi();
if (phaseEnd.getDay() == null) {
LOGGER.warn("No end date for phase {}/{} in {}!", stage, phaseSums.length, year);
final Stage s = new Stage();
s.setName(stageName);
stages.add(s);
break;
}
phaseStart = phaseEnd.getDay() + 1;
Stage s;
if (model == PhenologicalModelType.curve_grapevine_sw) {
// for this model, dates are computed for year N and N+1,
// but used on 1 year only
s = new Stage(stageName, phaseEnd.getDay() - DAYS_IN_YEAR);
} else {
s = new Stage(stageName, phaseEnd.getDay());
}
stages.add(s);
}
//
final AnnualStageData data = new AnnualStageData();
if (nbOfYears == 1 || allYears) {
data.setYear(year);
} else {
data.setYear(year + 1);
}
data.setStages(stages);
data.setComplete(stages.size() == phaseSums.length + 1);
phenologicalData.add(data);
}
LOGGER.traceExit();
return phenologicalData;
}
/**
* @param nb
* number of stages to compute.
*/
@XmlElement
public void setNbOfStages(final int nb) {
phaseSums = new Double[nb];
}
/**
* @param value
* Sum of effective temperature from sowing to emergence (°C),
* aka ST_Phase1_4 OR ST_Phase1_6.
*/
@XmlElement
public void setPhase1Sum(final Double value) {
setPhaseSum(1, value);
}
/**
* @param value
* Sum of effective temperature from emergence to maximum
* acceleration of leaf growth (°C), aka ST_Phase2_4, OR Sum of
* effective temperature from emergence to début tallage (°C),
* aka ST_Phase2_6.
*/
@XmlElement
public void setPhase2Sum(final Double value) {
setPhaseSum(2, value);
}
/**
* @param value
* Sum of effective temperature from maximum acceleration of leaf
* growth to flowering (°C), aka ST_Phase3_4, OR Sum of effective
* temperature from début tallage to epi 1 cm (°C), aka
* ST_Phase3_6.
*/
@XmlElement
public void setPhase3Sum(final Double value) {
final int phase = 3;
setPhaseSum(phase, value);
}
/**
* @param value
* Sum of effective temperature from flowering to physiological
* maturity (°C), aka ST_Phase4_4, OR Sum of effective
* temperature from epi 1 cm to meiose (°C), aka ST_Phase4_6.
*/
@XmlElement
public void setPhase4Sum(final Double value) {
final int phase = 4;
setPhaseSum(phase, value);
}
/**
* @param value
* Sum of effective temperature from meiose to floraison (°C),
* aka ST_Phase5_6.
*/
@XmlElement
public void setPhase5Sum(final Double value) {
final int phase = 5;
setPhaseSum(phase, value);
}
/**
* @param value
* Sum of effective temperature from floraison to physiological
* maturity (°C), aka ST_Phase6_6.
*/
@XmlElement
public void setPhase6Sum(final Double value) {
final int phase = 6;
setPhaseSum(phase, value);
}
/**
* @param phase
* phase number
* @param value
* Sum of effective temperature (°C)
*/
private void setPhaseSum(final int phase, final Double value) {
if (value == null) {
throw new IllegalArgumentException("phase" + phase
+ "Sum must not be null!");
}
if (phase < 1) {
throw new IllegalArgumentException("phase number must be >= 1");
}
if (phase > this.phaseSums.length) {
throw new IllegalArgumentException("phase number must be <= "
+ this.phaseSums.length);
}
this.phaseSums[phase - 1] = value;
}
@Override
public void setTimeScale(final TimeScale timeScale) {
// do nothing
}
}