View Javadoc
1   /*
2    * This file is part of Indicators.
3    *
4    * Indicators is free software: you can redistribute it and/or modify
5    * it under the terms of the GNU General Public License as published by
6    * the Free Software Foundation, either version 3 of the License, or
7    * (at your option) any later version.
8    *
9    * Indicators is distributed in the hope that it will be useful,
10   * but WITHOUT ANY WARRANTY; without even the implied warranty of
11   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12   * GNU General Public License for more details.
13   *
14   * You should have received a copy of the GNU General Public License
15   * along with Indicators. If not, see <https://www.gnu.org/licenses/>.
16   */
17  package fr.inrae.agroclim.indicators.model.data.soil;
18  
19  import java.text.SimpleDateFormat;
20  import java.util.ArrayList;
21  import java.util.Arrays;
22  import java.util.Calendar;
23  import java.util.Collection;
24  import java.util.Collections;
25  import java.util.Date;
26  import java.util.HashMap;
27  import java.util.HashSet;
28  import java.util.List;
29  import java.util.Map;
30  import java.util.Objects;
31  import java.util.Set;
32  import java.util.StringJoiner;
33  import java.util.TimeZone;
34  
35  import fr.inrae.agroclim.indicators.model.TimeScale;
36  import fr.inrae.agroclim.indicators.model.data.DataLoadingListenerHandler;
37  import fr.inrae.agroclim.indicators.model.data.Variable;
38  import fr.inrae.agroclim.indicators.model.data.climate.ClimaticDailyData;
39  import fr.inrae.agroclim.indicators.model.data.phenology.Stage;
40  import fr.inrae.agroclim.indicators.resources.Messages;
41  import lombok.Getter;
42  import lombok.Setter;
43  import lombok.extern.log4j.Log4j2;
44  
45  /**
46   * Calculator for Soil water balance.
47   *
48   * Last changed : $Date$
49   *
50   * @author $Author$
51   * @version $Revision$
52   */
53  @Log4j2
54  public final class SoilCalculator extends DataLoadingListenerHandler implements HasSoilCalculatorParams, SoilLoader {
55      /**
56       * UID for serialization.
57       */
58      private static final long serialVersionUID = 9107886175528703325L;
59  
60      /**
61       * Four stages.
62       */
63      private static final int FOUR = 4;
64  
65      /**
66       * Converting meter to centimeter.
67       */
68      private static final int M_TO_CM = 10;
69  
70      /**
71       * Third stage.
72       */
73      private static final int THREE = 3;
74  
75      /**
76       * Climatic daily data used to compute soil water balance.
77       */
78      @Setter
79      private List<ClimaticDailyData> climaticDailyData;
80  
81      /**
82       * Coefficient cultural en situation non stressée durant la période initiale
83       * de la saison [0-1] (semis-levée selon notre hypothèse).
84       */
85      @Setter
86      @Getter
87      private Double kcIni;
88  
89      /**
90       * Coefficient cultural en situation non stressée durant la late season
91       * [0-1] (à partir de floraison selon nous).
92       */
93      @Setter
94      @Getter
95      private Double kcLate;
96  
97      /**
98       * Coefficient cultural en situation non stressée durant la middle season
99       * [0-1] (épi 1 cm - floraison selon notre hypothèse).
100      */
101     @Setter
102     @Getter
103     private Double kcMid;
104 
105     /**
106      * Proportion of easily available water in the soil.
107      *
108      * Average fraction of TAW, pour obtenir RAW.
109      *
110      * http://www.kimberly.uidaho.edu/water/fao56/fao56.pdf#189
111      */
112     @Setter
113     private Double p;
114 
115     /**
116      * Réserve utile du sol (mm).
117      */
118     @Setter
119     private Double ru;
120 
121     /**
122      * Soil depth (cm).
123      */
124     @Setter
125     private Double soilDepth;
126 
127     /**
128      * Date of the 4 stages by year needed to compute SWB and R.
129      */
130     private List<Date> stages;
131 
132     /**
133      * Date of the 4 stages by year needed to compute SWB and R.
134      *
135      * Computed from stages.
136      */
137     private final Map<Integer, List<Date>> stagesByYear = new HashMap<>();
138 
139     /**
140      * Teneur en eau du sol à la capacité au champ (% volumique).
141      */
142     @Setter
143     @Getter
144     private Double swcFc;
145 
146     /**
147      * Teneur en eau du sol maxi (% volumique).
148      */
149     @Setter
150     @Getter
151     private Double swcMax;
152 
153     /**
154      * Teneur en eau du sol au point de flétrissement permanent (% volumique).
155      */
156     @Setter
157     @Getter
158     private Double swcWp;
159 
160     /**
161      * Constructor.
162      */
163     public SoilCalculator() {
164         // Do nothing
165     }
166 
167     @Override
168     public SoilCalculator clone() {
169         final SoilCalculator clone = new SoilCalculator();
170         if (climaticDailyData != null) {
171             clone.climaticDailyData = new ArrayList<>();
172             Collections.copy(clone.climaticDailyData, climaticDailyData);
173         }
174         clone.kcIni = kcIni;
175         clone.kcLate = kcLate;
176         clone.kcMid = kcMid;
177         clone.p = p;
178         clone.ru = ru;
179         clone.soilDepth = soilDepth;
180         if (stages != null) {
181             clone.stages = new ArrayList<>();
182             Collections.copy(clone.stages, stages);
183         }
184         clone.swcFc = swcFc;
185         clone.swcMax = swcMax;
186         clone.swcWp = swcWp;
187         return clone;
188     }
189 
190     /**
191      * Compute Kc.
192      *
193      * @param date date of climatic data
194      * @param previous date of previous stage
195      * @param current date of current stage
196      * @return Kc value
197      */
198     private double computeKc(final Date date, final Date previous,
199             final Date current) {
200         final long timeInPhase = current.getTime() - previous.getTime();
201         final double daysInPhase = Math.ceil(timeInPhase
202                 / (1000 * 3600 * 24.));
203         final long timeFromStage3 = date.getTime() - previous.getTime();
204         final double daysFromStage3 = timeFromStage3
205                 / (1000 * 3600 * 24.);
206         return (kcLate - kcMid) / daysInPhase * daysFromStage3 + kcLate;
207     }
208 
209     /**
210      * Compute SoilDailyData swc and waterReserve.
211      *
212      * @param soilDailyData daily data to update
213      * @param dailyData climatic data for the day
214      * @param stage1 stage 1
215      * @param stage2 stage 2
216      * @param stage3 stage 3
217      * @param stage4 stage 4
218      * @param firstDay is it the first day of climatic series
219      * @param rU Réserve utile du sol (mm).
220      * @param raw Readily available water en cm (au-delà, la plante commence à
221      * être stressée).
222      * @param rSat Réserve à saturation (mm).
223      * @param comptVeille compteur pour le ressuyage du jour précédent
224      * @param rVeille ressuyage du jour précédent
225      * @return new comptVeille
226      */
227     @SuppressWarnings("checkstyle:ParameterNumber")
228     private int computeSoilDailyData(final SoilDailyData soilDailyData,
229             final ClimaticDailyData dailyData, final Date stage1,
230             final Date stage2, final Date stage3, final Date stage4,
231             final boolean firstDay, final double rU, final double raw,
232             final double rSat, final int comptVeille, final double rVeille) {
233         final Date date = dailyData.getDate();
234         // Crop evapotranspiration under non-standard conditions
235         double etCAdj;
236         final double et0 = dailyData.getEtp();
237         final double rain = dailyData.getRain();
238         final double kC;
239         // Calcul de Kc (X valeurs pour X phases phénologiques).
240         // SWBressuyage : lignes 151-159
241         if (date.before(stage1) || date.equals(stage1)) {
242             kC = kcIni;
243         } else if (date.before(stage2)) {
244             kC = computeKc(date, stage1, stage2);
245         } else if (date.before(stage3) || date.equals(stage3)) {
246             kC = kcMid;
247         } else if (date.before(stage4)) {
248             kC = computeKc(date, stage3, stage4);
249         } else {
250             kC = kcIni;
251         }
252         // Calcul de Ks en journalier à partir de la réserve de la veille
253         // (si le déficit est supérieure a la RAW, on estime que cela a un
254         // effet stress sur la transpiration de la plante).
255         int compt;
256         double kS;
257         double r;
258         double ressuyage;
259         if (firstDay) {
260             // cas particulier du premier jour de la première année
261             kS = rU / (rU - raw);
262             if (kS > 1) {
263                 kS = 1;
264             }
265             //
266             etCAdj = kS * kC * et0;
267             // réserve en mm, bornée à la réserve à saturation
268             r = rU + rain - etCAdj;
269             if (r > rSat) {
270                 r = rSat;
271             }
272             // compteur pour le ressuyage du jour suivant
273             if (r <= rU) {
274                 compt = 0;
275             } else {
276                 compt = 1;
277             }
278         } else {
279             // cas des autres jours
280             if (comptVeille == 0) {
281                 // cas pas de ressuyage
282                 if (rVeille < raw) {
283                     kS = rVeille / (rU - raw);
284                 } else {
285                     kS = 1;
286                 }
287                 etCAdj = kS * kC * et0;
288                 r = rVeille + rain - etCAdj;
289             } else {
290                 // cas ressuyage
291                 etCAdj = kC * et0;
292                 switch (comptVeille) {
293                 case 1:
294                     // cas 1er jour de ressuyage
295                     ressuyage = Math.min((rSat - rU) / 2, rVeille - rU);
296                     // calcul de la réserve
297                     r = rVeille - ressuyage - etCAdj + rain;
298                     break;
299                 case 2:
300                     // cas 2e jour de ressuyage
301                     r = rU + rain - etCAdj;
302                     break;
303                 default:
304                     throw new RuntimeException("compt_veille != 1 && != 2");
305                 }
306             }
307             // réserve en mm, bornée à la réserve à saturation
308             if (r > rSat) {
309                 r = rSat;
310             }
311             // on calcule le compteur ressuyage pour le jour d'après
312             if (r <= rU) {
313                 compt = 0;
314             } else {
315                 switch (comptVeille) {
316                 case 0:
317                     compt = 1;
318                     break;
319                 case 1:
320                     if (rain == 0) {
321                         compt = 2;
322                     } else {
323                         compt = 1;
324                     }
325                     break;
326                 case 2:
327                     compt = 1;
328                     break;
329                 default:
330                     throw new RuntimeException(
331                             "compt_veille != 0 && != 1 && != 2");
332                 }
333             }
334         }
335         // Calcul de SWC, dans le cas où SWC_WP est renseigné.
336         // On ne considère pas la profondeur maximale d'enracinement mais la
337         // profondeur de sol
338         if (swcWp != null) {
339             final double swc = r / (soilDepth / M_TO_CM) + swcWp;
340             soilDailyData.setSwc(swc);
341         }
342         soilDailyData.setWaterReserve(r);
343         return compt;
344     }
345 
346     @Override
347     public Map<String, String> getConfigurationErrors() {
348         final Map<String, String> errors = new HashMap<>();
349         if (climaticDailyData == null) {
350             errors.put("soil.climaticDailyData", "error.evaluation.soil.climaticDailydata.missing");
351         } else if (climaticDailyData.isEmpty()) {
352             errors.put("soil.climaticDailyData", "error.evaluation.soil.climaticDailydata.empty");
353         }
354         // assert that the climatic data set is complete (1 data per day, all
355         // days per year, following years).
356         if (stages == null) {
357             errors.put("soil.stages", "error.evaluation.soil.stages.missing");
358         } else if (climaticDailyData != null && !climaticDailyData.isEmpty()) {
359             // check there are 4 stages for each year
360             setStagesByYear();
361             for (final Map.Entry<Integer, List<Date>> entry : stagesByYear.entrySet()) {
362                 final Integer year = entry.getKey();
363                 final List<Date> dates = entry.getValue();
364                 if (dates.isEmpty()) {
365                     LOGGER.error("The year {} does not have any stages.", year);
366                     errors.put("soil.stages", "error.evaluation.soil.stages.4stagesperyear " + year + " 0");
367                 } else if (dates.size() != Stage.FOUR) {
368                     LOGGER.error("The year {} does not have {} stages {}, filling stages.", year, Stage.FOUR,
369                             stagesByYear.get(year));
370                     for (int nbOfStages = dates.size(); nbOfStages <= Stage.FOUR; nbOfStages++) {
371                         stagesByYear.get(year).add(stagesByYear.get(year).get(nbOfStages - 1));
372                     }
373 
374                 }
375             }
376         }
377         if (ru == null) {
378             if (swcFc == null) {
379                 errors.put("soil.swcFc", "error.evaluation.soil.swcFc.missing");
380             }
381             if (swcWp == null) {
382                 errors.put("soil.swcWp", "error.evaluation.soil.swcWp.missing");
383             }
384             if (soilDepth == null) {
385                 errors.put("soil.soilDepth", "error.evaluation.soil.soilDepth.missing");
386             }
387         }
388         if (kcIni == null) {
389             errors.put("soil.kcIni", "error.evaluation.soil.kCini.missing");
390         }
391         if (kcMid == null) {
392             errors.put("soil.kcMid", "error.evaluation.soil.kCmid.missing");
393         }
394         if (kcLate == null) {
395             errors.put("soil.kcLate", "error.evaluation.soil.kClate.missing");
396         }
397         if (p == null) {
398             errors.put("soil.p", "error.evaluation.soil.p.missing");
399         }
400         if (swcMax == null) {
401             errors.put("soil.swcMax", "error.evaluation.soil.swcMax.missing");
402         }
403 
404         if (errors.isEmpty()) {
405             return null;
406         }
407         return errors;
408     }
409 
410     @Override
411     public Collection<String> getMissingVariables() {
412         throw new RuntimeException("Not implemented for soil!");
413     }
414 
415     @Override
416     public Set<Variable> getProvidedVariables() {
417         final Set<Variable> variables = new HashSet<>();
418         variables.add(Variable.SOILWATERCONTENT);
419         variables.add(Variable.WATER_RESERVE);
420         return variables;
421     }
422 
423     /**
424      * @return Readily available water en cm (au-delà, la plante commence à être
425      *         stressée).
426      */
427     private double getRAW() {
428         return p * getRu();
429     }
430 
431     /**
432      * Calcul de la réserve à saturation.
433      *
434      * @return Réserve à saturation (mm).
435      */
436     private double getRsat() {
437         if (swcMax == null) {
438             throw new IllegalStateException(
439                     "swcMax must not be null to compute Rsat!");
440         }
441         if (swcWp == null) {
442             throw new IllegalStateException(
443                     "swcWp must not be null to compute Rsat!");
444         }
445         if (soilDepth == null) {
446             throw new IllegalStateException(
447                     "soilDepth must not be null to compute Rsat!");
448         }
449         return (swcMax - swcWp) * (soilDepth / M_TO_CM);
450     }
451 
452     /**
453      * On calcule la RU sur la profondeur du sol à partir de soilDepth, SWC_FC
454      * et SWC_WP si elle n'a pas été renseignée.
455      *
456      * @return Réserve utile du sol (mm).
457      */
458     private double getRu() {
459         if (ru == null) {
460             if (swcFc == null) {
461                 throw new IllegalStateException(
462                         "swcFc must not be null to compute RU!");
463             }
464             if (swcWp == null) {
465                 throw new IllegalStateException(
466                         "swcWp must not be null to compute RU!");
467             }
468             if (soilDepth == null) {
469                 throw new IllegalStateException(
470                         "soilDepth must not be null to compute RU!");
471             }
472             return (swcFc - swcWp) * (soilDepth / M_TO_CM);
473         }
474         return ru;
475     }
476 
477     @Override
478     public Set<Variable> getVariables() {
479         final Set<Variable> variables = new HashSet<>();
480         variables.addAll(Arrays.asList(Variable.ETP, Variable.RAIN));
481         return variables;
482     }
483 
484     /**
485      * Raise exception when there is no stage for a year.
486      *
487      * @param year the year for which there is no stage.
488      */
489     private void handleMissingStages(final Integer year) {
490         final SimpleDateFormat sdf = new SimpleDateFormat("yyyy-MM-dd");
491         sdf.setTimeZone(TimeZone.getTimeZone("UTC"));
492         final StringJoiner sj = new StringJoiner(" ");
493         sj.add(Messages.format("error.soilcalculator.nostage", "" + year));
494         stages.stream() //
495         .filter(Objects::nonNull) //
496         .map(sdf::format) //
497         .forEach(sj::add);
498         throw new IllegalStateException(sj.toString());
499     }
500 
501     @Override
502     public List<SoilDailyData> load() {
503         if (climaticDailyData == null) {
504             throw new IllegalStateException("No climatic daily data are set!");
505         }
506         if (climaticDailyData.isEmpty()) {
507             throw new IllegalStateException("Climatic daily data is empty!");
508         }
509         setStagesByYear();
510         final double rU = getRu();
511         final double raw = getRAW();
512         final double rSat = getRsat();
513         final List<SoilDailyData> data = new ArrayList<>();
514         double rVeille = 0.;
515         int comptVeille = 0;
516         boolean firstDay = true;
517         Integer year = climaticDailyData.get(0).getYear();
518         if (stagesByYear.get(year) == null) {
519             handleMissingStages(year);
520         }
521         Date stage1 = null;
522         Date stage2 = null;
523         Date stage3 = null;
524         Date stage4 = null;
525         if (stagesByYear.get(year).size() >= FOUR) {
526             stage1 = stagesByYear.get(year).get(0);
527             stage2 = stagesByYear.get(year).get(1);
528             stage3 = stagesByYear.get(year).get(2);
529             stage4 = stagesByYear.get(year).get(THREE);
530         }
531         // KC for the beginning of the year,
532         // this could be kcLate for the 2nd year
533         for (final ClimaticDailyData dailyData : climaticDailyData) {
534             if (year != dailyData.getYear()) {
535                 year = dailyData.getYear();
536                 if (stagesByYear.get(year) == null) {
537                     handleMissingStages(year);
538                 }
539                 if (stagesByYear.get(year).size() >= FOUR) {
540                     stage1 = stagesByYear.get(year).get(0);
541                     stage2 = stagesByYear.get(year).get(1);
542                     stage3 = stagesByYear.get(year).get(2);
543                     stage4 = stagesByYear.get(year).get(THREE);
544                 } else {
545                     stage1 = null;
546                     stage2 = null;
547                     stage3 = null;
548                     stage4 = null;
549                 }
550             }
551             final SoilDailyData soilDailyData = new SoilDailyData();
552             soilDailyData.setYear(dailyData.getYear());
553             soilDailyData.setMonth(dailyData.getMonth());
554             soilDailyData.setDay(dailyData.getDay());
555             if (stage4 != null && stage3 != null && stage2 != null
556                     && stage1 != null) {
557                 comptVeille = computeSoilDailyData(soilDailyData, dailyData,
558                         stage1, stage2, stage3, stage4, firstDay, rU, raw, rSat,
559                         comptVeille, rVeille);
560                 rVeille = soilDailyData.getWaterReserve();
561             }
562             data.add(soilDailyData);
563             if (firstDay) {
564                 firstDay = false;
565             }
566         }
567         return data;
568     }
569 
570     /**
571      * @param fourStages
572      *            the 4 stages (per year) needed to compute SWB and R
573      */
574     public void setStages(final List<Date> fourStages) {
575         this.stages = fourStages;
576         stagesByYear.clear();
577     }
578 
579     /**
580      * Divide all stages by year.
581      */
582     private void setStagesByYear() {
583         if (stagesByYear.isEmpty()) {
584             final Calendar cal = Calendar.getInstance(TimeZone.getTimeZone("UTC"));
585             Integer year;
586             for (int i = 0; i < stages.size(); i++) {
587                 final Date date = stages.get(i);
588                 cal.setTime(date);
589                 year = cal.get(Calendar.YEAR);
590                 stagesByYear.computeIfAbsent(year, y -> new ArrayList<>());
591                 stagesByYear.get(year).add(date);
592                 Collections.sort(stagesByYear.get(year));
593             }
594         }
595     }
596 
597     @Override
598     public void setTimeScale(final TimeScale timeScale) {
599         // do nothing
600     }
601 
602 }