Improved leaft wetness calculation algorithm

src/main/java/no/bioforsk/vips/entity/WeatherObservation.java

@@ -41,6 +41,20 @@ public class WeatherObservation implements Comparable{
     /** Creates a new instance of WeatherObservation */
     public WeatherObservation() {
     }
+    
+    /**
+     * 
+     * @param timeMeasured
+     * @param elementMeasurementTypeId
+     * @param logIntervalId
+     * @param value 
+     */
+    public WeatherObservation(Date timeMeasured, String elementMeasurementTypeId, Integer logIntervalId, Double value){
+        this.timeMeasured = timeMeasured;
+        this.elementMeasurementTypeId = elementMeasurementTypeId;
+        this.logIntervalId = logIntervalId;
+        this.value = value;
+    }
         
     public void setTimeMeasured(Date timeMeasured) { this.timeMeasured = timeMeasured; }
     public Date getTimeMeasured() { return this.timeMeasured; }
@@ -58,6 +72,7 @@ public class WeatherObservation implements Comparable{
     /**
      * Compares by time measured in ascending order
      */
+    @Override
     public int compareTo(Object o) {
         WeatherObservation other = (WeatherObservation) o;
         return this.compareTo(other);

src/main/java/no/bioforsk/vips/util/WeatherElements.java

@@ -54,7 +54,10 @@ public class WeatherElements {
      * Leaf wetness (Minutes per hour)
      */
     public static final String LEAF_WETNESS = "BT";
-
+    /**
+     * 
+     */
+    public static final String WIND_SPEED_2M = "FM2";
     
     
 }

src/main/java/no/bioforsk/vips/util/WeatherUtil.java

@@ -330,4 +330,169 @@ public class WeatherUtil {
         cal.set(Calendar.MILLISECOND,0);
         return cal.getTime();
     }
+
+    /**
+     * Given input data, attempts to calculate leaf wetness. Does not overwrite provided leaf wetness data.
+     * Priority (depending on what data are provided):
+     * <ol>
+     * <li>Provided leaf wetness data</li>
+     * <li>Nærstad's leaf wetness calculation algorithm</li>
+     * <li>Trapman's simple leaf wetness calculation algorithm</li>
+     * </ol>
+     * <p>Requires that at least TM, RR, and UM are complete. And that all lists are sorted!</p>
+     * 
+     * @param BT Leaf wetness (minutes/hour)
+     * @param TM Mean temperature (hourly mean, degrees Celcius)
+     * @param RR Rainfall (mm, hourly aggregate)
+     * @param FM2 Wind speed at 2m height (average, m/s)
+     * @param Q0 Global radiation (unit??)
+     * @param UM Relative humidity (average, %)
+     * @return 
+     */
+    public List<WeatherObservation> calculateLeafWetnessHourSeriesBestEffort(List<WeatherObservation> BT, List<WeatherObservation> TM, List<WeatherObservation> RR, List<WeatherObservation> Q0, List<WeatherObservation> FM2, List<WeatherObservation> UM) throws ConfigValidationException {
+        // If BT is same length as TM, return BT unmodified
+        if(BT.size() == TM.size())
+        {
+            return BT;
+        }
+        
+        // Possible TODO: Validate that TM==RR==UM
+        // First: Find time and index for when to start calculating
+        Integer listIndex = BT.size();
+        Date startTime = TM.get(listIndex).getTimeMeasured();
+        
+        // If all parameters are complete, we can use Nærstad's algorithm. Otherwise, if we only have UM,
+        // we must resort to Trapman's algorithm
+        List<WeatherObservation> calculatedBT;
+        if(TM.size() == RR.size() && TM.size() == Q0.size() && TM.size() == FM2.size() && TM.size() == UM.size())
+        {
+            // The algorithm uses the first hour of input values as "startup"
+            // So if we have some leaf wetness, we can provide the first hour.
+            if(listIndex > 0)
+            {
+                calculatedBT = this.calculateLeafWetnessHourSeriesNaerstad(
+                                TM.subList(listIndex-1, TM.size()),
+                                RR.subList(listIndex-1, TM.size()),
+                                Q0.subList(listIndex-1, TM.size()),
+                                FM2.subList(listIndex-1, TM.size()),
+                                UM.subList(listIndex-1, TM.size()),
+                                BT.get(listIndex-1)
+                                );
+            }
+            // Otherwise, we have to create a dry hour as the first hour
+            else
+            {
+                WeatherObservation emptyObs = new WeatherObservation(
+                                                startTime,
+                                                WeatherElements.LEAF_WETNESS,
+                                                WeatherObservation.LOG_INTERVAL_ID_1H,
+                                                0d
+                                                );
+                calculatedBT = new ArrayList<>();
+                calculatedBT.add(emptyObs);
+                calculatedBT.addAll(this.calculateLeafWetnessHourSeriesNaerstad(TM, RR, Q0, FM2, UM, null));
+            }
+        }
+        else if(TM.size() == RR.size() && TM.size() == UM.size())
+        {
+            calculatedBT = this.calculateLeafWetnessHourSeriesSimple(
+                    TM.subList(listIndex, TM.size()), 
+                    RR.subList(listIndex, TM.size()),
+                    UM.subList(listIndex, TM.size())
+            );
+        }
+        else
+        {
+            throw new ConfigValidationException("Missing weather data");
+        }
+        BT.addAll(calculatedBT);
+        return BT;
+    }
+
+    /**
+     * 
+     * @param temperature
+     * @param precipitation
+     * @param globalRadiation
+     * @param wind2m
+     * @param relativeHumidity
+     * @param precedingLeafWetness
+     * @return Leaf wetness for the same time period as provided by input data, 
+     * EXCEPT for the first hour, as that is used for "startup" of the algorithm.
+     * So the size of the returned list is 1 less than the input lists.
+     */
+    private List<WeatherObservation> calculateLeafWetnessHourSeriesNaerstad(
+            List<WeatherObservation> temperature, 
+            List<WeatherObservation> precipitation, 
+            List<WeatherObservation> globalRadiation, 
+            List<WeatherObservation> wind2m, 
+            List<WeatherObservation> relativeHumidity,
+            WeatherObservation precedingLeafWetness) 
+    {
+        // TODO Validation
+        
+        List<WeatherObservation> calculatedLeafWetnessSeries = new ArrayList<>();
+        Double lastLatentHeatFlux = this.calculateLatentHeatFlux(
+                temperature.get(0).getValue(), 
+                relativeHumidity.get(0).getValue(), 
+                wind2m.get(0).getValue(), 
+                globalRadiation.get(0).getValue()
+        );
+        Integer precedingLeafWetnessValue = precedingLeafWetness != null ?  (int) precedingLeafWetness.getValue() : 0;
+        // Iterating list with index (instead of using Iterator pattern) has 
+        // performance gains on very long lists
+        for(int i=1;i<temperature.size();i++)
+        {
+            WeatherObservation calculatedLeafWetness = new WeatherObservation(
+                    temperature.get(i).getTimeMeasured(),
+                    WeatherElements.LEAF_WETNESS,
+                    WeatherObservation.LOG_INTERVAL_ID_1H,
+                    (double) this.calculateLeafWetnessNaerstad(lastLatentHeatFlux, precedingLeafWetnessValue, precipitation.get(i).getValue())
+            );
+            calculatedLeafWetnessSeries.add(calculatedLeafWetness);
+        }
+        return calculatedLeafWetnessSeries;
+    }
+    
+    /**
+     * Nærstad's model of leaf wetness calculation
+     * @param lastLatentHeatFlux
+     * @param lastLeafWetness
+     * @param precipitation
+     * @return
+     */
+    public Long calculateLeafWetnessNaerstad(Double lastLatentHeatFlux, Integer lastLeafWetness, Double precipitation) {
+
+        Long DLeafWetness = Math.round(lastLeafWetness - lastLatentHeatFlux);
+        if (DLeafWetness < 0) {
+            DLeafWetness = 0l;
+        } else if (DLeafWetness > 60) {
+            DLeafWetness = 60l;
+        }
+
+
+        Long leafWetness = Math.round(DLeafWetness + precipitation * 150);
+        if (leafWetness < 0) {
+            leafWetness = 0l;
+        } else if (leafWetness > 60) {
+            leafWetness = 60l;
+        }
+
+        return leafWetness;
+    }
+    
+    /**
+     * Part of the Nærstad leaf wetness calculation
+     * @return
+     */
+    public Double calculateLatentHeatFlux(Double temperature, Double relativeHumidity, Double wind2m, Double globalRadiation) {
+        // s(t) = 0.61078*EKSP(17.269* T(t)/( T(t)+237.3))*4.097*10000/(T(t)+237.3)2
+        Double s = 0.61078 * Math.exp(17.269 * temperature / (temperature + 237.3)) * 4.097 * 10000 / Math.pow(temperature + 237.3, 2);
+        // resistance(t) = 307*(0.07/( (windFF(t)   + 0.1) / 2 ))0.5
+        Double resistance = 307 * Math.pow(0.07 / ((wind2m + 0.1) / 2), 0.5);
+        // LE(t)  = -( (s(t) * Rn(t)  + (12(WVD (t) )/( resistance(t)))/(s(t) +0.64)
+        Double LE = (s * globalRadiation + (12 * this.getWVD(temperature, relativeHumidity) / resistance)) / (s + 0.64);
+
+        return LE;
+    }
 }

src/test/java/no/bioforsk/vips/util/WeatherUtilTest.java

@@ -26,10 +26,13 @@ import java.util.Collections;
 import java.util.Date;
 import java.util.List;
 import java.util.TimeZone;
+import java.util.logging.Level;
+import java.util.logging.Logger;
 import static junit.framework.Assert.assertEquals;
 import static junit.framework.Assert.fail;
 import junit.framework.TestCase;
 import no.bioforsk.vips.entity.WeatherObservation;
+import no.bioforsk.vips.model.ConfigValidationException;
 import org.codehaus.jackson.JsonFactory;
 import org.codehaus.jackson.JsonNode;
 import org.codehaus.jackson.JsonParser;
@@ -155,6 +158,80 @@ public class WeatherUtilTest extends TestCase {
         assertEquals(expMin, minResults.get(0).getValue());
         assertEquals(expMax, maxResults.get(0).getValue());
     }
+    
+    public void testCalculateLeafWetnessHourSeriesNaerstad()
+    {
+        System.out.println("testCalculateLeafWetnessHourSeriesNaerstad");
+        List<WeatherObservation> allObservations = this.getObservations("/weatherDataLeafWetnessCalc.json");
+        List<WeatherObservation> TM = new ArrayList<>();
+        List<WeatherObservation> UM = new ArrayList<>();
+        List<WeatherObservation> FM2 = new ArrayList<>();
+        List<WeatherObservation> Q0 = new ArrayList<>();
+        List<WeatherObservation> RR = new ArrayList<>();
+        List<WeatherObservation> BT = new ArrayList<>();
+        // Add collections here
+        for(WeatherObservation obs:allObservations)
+        {
+            if(obs.getElementMeasurementTypeId().equals(WeatherElements.TEMPERATURE_MEAN))
+                TM.add(obs);
+            else if(obs.getElementMeasurementTypeId().equals(WeatherElements.RELATIVE_HUMIDITY))
+                UM.add(obs);
+            else if(obs.getElementMeasurementTypeId().equals(WeatherElements.WIND_SPEED_2M))
+                FM2.add(obs);
+            else if(obs.getElementMeasurementTypeId().equals(WeatherElements.GLOBAL_RADIATION))
+                Q0.add(obs);
+            else if(obs.getElementMeasurementTypeId().equals(WeatherElements.PRECIPITATION))
+                RR.add(obs);
+            else if(obs.getElementMeasurementTypeId().equals(WeatherElements.LEAF_WETNESS))
+                BT.add(obs);
+        }
+        
+        // Trying "Best effort"
+        List<WeatherObservation> BTtest = BT.size() > 1 ? new ArrayList(BT.subList(0, BT.size()/2)) : new ArrayList(BT);
+        System.out.println("BTtest.size=" + BTtest.size());
+        System.out.println("BT.size=" + BT.size());
+        
+        WeatherUtil instance = new WeatherUtil();
+        try {
+            // Naerstad should be called now
+            List<WeatherObservation> result = instance.calculateLeafWetnessHourSeriesBestEffort(BTtest, TM, RR, Q0, FM2, UM);
+            System.out.println("result.size=" + result.size());
+            // Test 1 // Result have same length as input data
+            assertEquals(TM.size(), result.size());
+            // Starts with same timestamp
+            assertEquals(TM.get(0).getTimeMeasured(), result.get(0).getTimeMeasured());
+            // Ends with same timestamp
+            assertEquals(TM.get(TM.size()-1).getTimeMeasured(), result.get(TM.size()-1).getTimeMeasured());
+            /*
+            for(int i=0; i< result.size();i++)
+            {
+                WeatherObservation calc = result.get(i);
+                WeatherObservation orig = i < BT.size() ? BT.get(i) : null;
+                System.out.println(calc.getElementMeasurementTypeId() + " [" +calc.getTimeMeasured() + "]: " + calc.getValue() + " | " + (orig != null ? orig.getValue():"NULL"));
+            }*/
+            
+            // Trapman should be called now
+            Q0 = new ArrayList<>();
+            BTtest = BT.size() > 1 ? new ArrayList(BT.subList(0, BT.size()/2)) : new ArrayList(BT);
+            result = instance.calculateLeafWetnessHourSeriesBestEffort(BTtest, TM, RR, Q0, FM2, UM);
+            System.out.println("result.size=" + result.size());
+            // Test 1 // Result have same length as input data
+            assertEquals(TM.size(), result.size());
+            // Starts with same timestamp
+            assertEquals(TM.get(0).getTimeMeasured(), result.get(0).getTimeMeasured());
+            // Ends with same timestamp
+            assertEquals(TM.get(TM.size()-1).getTimeMeasured(), result.get(TM.size()-1).getTimeMeasured());
+            /*
+            for(int i=0; i< result.size();i++)
+            {
+                WeatherObservation calc = result.get(i);
+                WeatherObservation orig = i < BT.size() ? BT.get(i) : null;
+                System.out.println(calc.getElementMeasurementTypeId() + " [" +calc.getTimeMeasured() + "]: " + calc.getValue() + " | " + (orig != null ? orig.getValue():"NULL"));
+            }*/
+        } catch (ConfigValidationException ex) {
+            Logger.getLogger(WeatherUtilTest.class.getName()).log(Level.SEVERE, null, ex);
+        }
+    }
 
     /**
      * Test of normalizeToExactDate method, of class WeatherUtil.