Merge pull request #62 from Derpthemeus/add-eto-adjustment-method

Add ETo AdjustmentMethod

routes/adjustmentMethods/EToAdjustmentMethod.spec.ts

@@ -0,0 +1,43 @@
+import * as moment from "moment";
+import { expect } from "chai";
+import { GeoCoordinates } from "../../types";
+import { calculateETo, EToData } from "./EToAdjustmentMethod";
+
+
+const testData: TestData[] = require( "../../test/etoTest.json" );
+
+describe( "ETo AdjustmentMethod", () => {
+	describe( "Should correctly calculate ETo", async () => {
+		for ( const locationData of testData ) {
+			it( "Using data from " + locationData.description, async () => {
+				let date = moment.unix( locationData.startTimestamp );
+				for ( const entry of locationData.entries ) {
+					const etoData: EToData = {
+						...entry.data,
+						precip: 0,
+						periodStartTime: date.unix(),
+						weatherProvider: "mock"
+					};
+					const calculatedETo = calculateETo( etoData, locationData.elevation, locationData.coordinates );
+					// Allow a small margin of error for rounding, unit conversions, and approximations.
+					expect( calculatedETo ).approximately( entry.eto, 0.003 );
+
+					date = date.add( 1, "days" );
+				}
+			} );
+		}
+	} );
+} );
+
+interface TestData {
+	description: string;
+	source: string;
+	startTimestamp: number;
+	elevation: number;
+	coordinates: GeoCoordinates;
+	entries: {
+		eto: number,
+		/** This is not actually full EToData - it is missing `timestamp`, `weatherProvider`, and `precip`. */
+		data: EToData
+	}[];
+}

routes/adjustmentMethods/EToAdjustmentMethod.ts

@@ -0,0 +1,234 @@
+import * as SunCalc from "suncalc";
+import * as moment from "moment";
+import { AdjustmentMethod, AdjustmentMethodResponse, AdjustmentOptions } from "./AdjustmentMethod";
+import { GeoCoordinates, WateringData, WeatherProviderId } from "../../types";
+import { WeatherProvider } from "../weatherProviders/WeatherProvider";
+
+
+/**
+ * Calculates how much watering should be scaled based on weather and adjustment options by comparing the recent
+ * potential ETo to the baseline potential ETo that the watering program was designed for.
+ */
+async function calculateEToWateringScale(
+	adjustmentOptions: EToScalingAdjustmentOptions,
+	wateringData: WateringData | undefined,
+	coordinates: GeoCoordinates,
+	weatherProvider: WeatherProvider
+): Promise< AdjustmentMethodResponse > {
+
+	// Temporarily disabled since OWM forecast data is checking if rain is forecasted for 3 hours in the future.
+	/*
+	if ( wateringData && wateringData.raining ) {
+		return {
+			scale: 0,
+			rawData: { raining: 1 }
+		}
+	}
+	 */
+
+	// This will throw an error message if ETo data cannot be retrieved.
+	const etoData: EToData = await weatherProvider.getEToData( coordinates );
+
+	let baseETo: number;
+	// Default elevation is based on data from https://www.pnas.org/content/95/24/14009.
+	let elevation = 600;
+
+	if ( adjustmentOptions && "baseETo" in adjustmentOptions ) {
+		baseETo = adjustmentOptions.baseETo
+	} else {
+		throw "A baseline potential ETo must be provided.";
+	}
+
+	if ( adjustmentOptions && "elevation" in adjustmentOptions ) {
+		elevation = adjustmentOptions.elevation;
+	}
+
+	const eto: number = calculateETo( etoData, elevation, coordinates );
+
+	const scale =  Math.floor( Math.min( Math.max( 0, ( eto - etoData.precip ) / baseETo * 100 ), 200 ) );
+	return {
+		scale: scale,
+		rawData: {
+			eto: Math.round( eto * 1000) / 1000,
+			radiation: Math.round( etoData.solarRadiation * 100) / 100,
+			minT: Math.round( etoData.minTemp ),
+			maxT: Math.round( etoData.maxTemp ),
+			minH: Math.round( etoData.minHumidity ),
+			maxH: Math.round( etoData.maxHumidity ),
+			wind: Math.round( etoData.windSpeed * 10 ) / 10
+		}
+	}
+}
+
+/* The implementation of this algorithm was guided by a step-by-step breakdown
+	(http://edis.ifas.ufl.edu/pdffiles/ae/ae45900.pdf) */
+/**
+ * Calculates the reference potential evapotranspiration using the Penman-Monteith (FAO-56) method
+ * (http://www.fao.org/3/X0490E/x0490e07.htm).
+ *
+ * @param etoData The data to calculate the ETo with.
+ * @param elevation The elevation above sea level of the watering site (in feet).
+ * @param coordinates The coordinates of the watering site.
+ * @return The reference potential evapotranspiration (in inches per day).
+ */
+export function calculateETo( etoData: EToData, elevation: number, coordinates: GeoCoordinates ): number {
+	// Convert to Celsius.
+	const minTemp = ( etoData.minTemp - 32 ) * 5 / 9;
+	const maxTemp = ( etoData.maxTemp - 32 ) * 5 / 9;
+	// Convert to meters.
+	elevation = elevation / 3.281;
+	// Convert to meters per second.
+	const windSpeed = etoData.windSpeed / 2.237;
+	// Convert to megajoules.
+	const solarRadiation = etoData.solarRadiation * 3.6;
+
+	const avgTemp = ( maxTemp + minTemp ) / 2;
+
+	const saturationVaporPressureCurveSlope = 4098 * 0.6108 * Math.exp( 17.27 * avgTemp / ( avgTemp + 237.3 ) ) / Math.pow( avgTemp + 237.3, 2 );
+
+	const pressure = 101.3 * Math.pow( ( 293 - 0.0065 * elevation ) / 293, 5.26 );
+
+	const psychrometricConstant = 0.000665 * pressure;
+
+	const deltaTerm = saturationVaporPressureCurveSlope / ( saturationVaporPressureCurveSlope + psychrometricConstant * ( 1 + 0.34 * windSpeed ) );
+
+	const psiTerm = psychrometricConstant / ( saturationVaporPressureCurveSlope + psychrometricConstant * ( 1 + 0.34 * windSpeed ) );
+
+	const tempTerm = ( 900 / ( avgTemp + 273 ) ) * windSpeed;
+
+	const minSaturationVaporPressure = 0.6108 * Math.exp( 17.27 * minTemp / ( minTemp + 237.3 ) );
+
+	const maxSaturationVaporPressure = 0.6108 * Math.exp( 17.27 * maxTemp / ( maxTemp + 237.3 ) );
+
+	const avgSaturationVaporPressure = ( minSaturationVaporPressure + maxSaturationVaporPressure ) / 2;
+
+	const actualVaporPressure = ( minSaturationVaporPressure * etoData.maxHumidity / 100 + maxSaturationVaporPressure * etoData.minHumidity / 100 ) / 2;
+
+	const dayOfYear = moment.unix( etoData.periodStartTime ).dayOfYear();
+
+	const inverseRelativeEarthSunDistance = 1 + 0.033 * Math.cos( 2 * Math.PI / 365 * dayOfYear );
+
+	const solarDeclination = 0.409 * Math.sin( 2 * Math.PI / 365 * dayOfYear - 1.39 );
+
+	const latitudeRads = Math.PI / 180 * coordinates[ 0 ];
+
+	const sunsetHourAngle = Math.acos( -Math.tan( latitudeRads ) * Math.tan( solarDeclination ) );
+
+	const extraterrestrialRadiation = 24 * 60 / Math.PI * 0.082 * inverseRelativeEarthSunDistance * ( sunsetHourAngle * Math.sin( latitudeRads ) * Math.sin( solarDeclination ) + Math.cos( latitudeRads ) * Math.cos( solarDeclination ) * Math.sin( sunsetHourAngle ) );
+
+	const clearSkyRadiation = ( 0.75 + 2e-5 * elevation ) * extraterrestrialRadiation;
+
+	const netShortWaveRadiation = ( 1 - 0.23 ) * solarRadiation;
+
+	const netOutgoingLongWaveRadiation = 4.903e-9 * ( Math.pow( maxTemp + 273.16, 4 ) + Math.pow( minTemp + 273.16, 4 ) ) / 2 * ( 0.34 - 0.14 * Math.sqrt( actualVaporPressure ) ) * ( 1.35 * solarRadiation / clearSkyRadiation - 0.35);
+
+	const netRadiation = netShortWaveRadiation - netOutgoingLongWaveRadiation;
+
+	const radiationTerm = deltaTerm * 0.408 * netRadiation;
+
+	const windTerm = psiTerm * tempTerm * ( avgSaturationVaporPressure - actualVaporPressure );
+
+	return ( windTerm + radiationTerm ) / 25.4;
+}
+
+/**
+ * Approximates the wind speed at 2 meters using the wind speed measured at another height.
+ * @param speed The wind speed measured at the specified height (in miles per hour).
+ * @param height The height of the measurement (in feet).
+ * @returns The approximate wind speed at 2 meters (in miles per hour).
+ */
+export function standardizeWindSpeed( speed: number, height: number ) {
+	return speed * 4.87 / Math.log( 67.8 * height / 3.281 - 5.42 );
+}
+
+/* For hours where the Sun is too low to emit significant radiation, the formula for clear sky isolation will yield a
+ * negative value. "radiationStart" marks the times of day when the Sun will rise high for solar isolation formula to
+ * become positive, and "radiationEnd" marks the time of day when the Sun sets low enough that the equation will yield
+ * a negative result. For any times outside of these ranges, the formula will yield incorrect results (they should be
+ * clamped at 0 instead of being negative).
+ */
+SunCalc.addTime( Math.asin( 30 / 990 ) * 180 / Math.PI, "radiationStart", "radiationEnd" );
+
+/**
+ * Approximates total solar radiation for a day given cloud coverage information using a formula from
+ * http://www.shodor.org/os411/courses/_master/tools/calculators/solarrad/
+ * @param cloudCoverInfo Information about the cloud coverage for several periods that span the entire day.
+ * @param coordinates The coordinates of the location the data is from.
+ * @return The total solar radiation for the day (in kilowatt hours per square meter per day).
+ */
+export function approximateSolarRadiation(cloudCoverInfo: CloudCoverInfo[], coordinates: GeoCoordinates ): number {
+	return cloudCoverInfo.reduce( ( total, window: CloudCoverInfo ) => {
+		const radiationStart: moment.Moment = moment( SunCalc.getTimes( window.endTime.toDate(), coordinates[ 0 ], coordinates[ 1 ])[ "radiationStart" ] );
+		const radiationEnd: moment.Moment = moment( SunCalc.getTimes( window.startTime.toDate(), coordinates[ 0 ], coordinates[ 1 ])[ "radiationEnd" ] );
+
+		// Clamp the start and end times of the window within time when the sun was emitting significant radiation.
+		const startTime: moment.Moment = radiationStart.isAfter( window.startTime ) ? radiationStart : window.startTime;
+		const endTime: moment.Moment = radiationEnd.isBefore( window.endTime ) ? radiationEnd: window.endTime;
+
+		// The length of the window that will actually be used (in hours).
+		const windowLength = ( endTime.unix() - startTime.unix() ) / 60 / 60;
+
+		// Skip the window if there is no significant radiation during the time period.
+		if ( windowLength <= 0 ) {
+			return total;
+		}
+
+		const startPosition = SunCalc.getPosition( startTime.toDate(), coordinates[ 0 ], coordinates[ 1 ] );
+		const endPosition = SunCalc.getPosition( endTime.toDate(), coordinates[ 0 ], coordinates[ 1 ] );
+		const solarElevationAngle = ( startPosition.altitude + endPosition.altitude ) / 2;
+
+		// Calculate radiation and convert from watts to kilowatts.
+		const clearSkyIsolation = ( 990 * Math.sin( solarElevationAngle ) - 30 ) / 1000 * windowLength;
+
+		return total + clearSkyIsolation * ( 1 - 0.75 * Math.pow( window.cloudCover, 3.4 ) );
+	}, 0 );
+}
+
+export interface EToScalingAdjustmentOptions extends AdjustmentOptions {
+	/** The watering site's height above sea level (in feet). */
+	elevation?: number;
+	/** Baseline potential ETo (in inches per day). */
+	baseETo?: number;
+}
+
+/** Data about the cloud coverage for a period of time. */
+export interface CloudCoverInfo {
+	/** The start of this period of time. */
+	startTime: moment.Moment;
+	/** The end of this period of time. */
+	endTime: moment.Moment;
+	/** The average fraction of the sky covered by clouds during this time period. */
+	cloudCover: number;
+}
+
+/**
+ * Data used to calculate ETo. This data should be taken from a 24 hour time window.
+ */
+export interface EToData {
+	/** The WeatherProvider that generated this data. */
+	weatherProvider: WeatherProviderId;
+	/** The Unix epoch seconds timestamp of the start of this 24 hour time window. */
+	periodStartTime: number;
+	/** The minimum temperature over the time period (in Fahrenheit). */
+	minTemp: number;
+	/** The maximum temperature over the time period (in Fahrenheit). */
+	maxTemp: number;
+	/** The minimum relative humidity over the time period (as a percentage). */
+	minHumidity: number;
+	/** The maximum relative humidity over the time period (as a percentage). */
+	maxHumidity: number;
+	/** The solar radiation, accounting for cloud coverage (in kilowatt hours per square meter per day). */
+	solarRadiation: number;
+	/**
+	 * The average wind speed measured at 2 meters over the time period (in miles per hour). A measurement taken at a
+	 * different height can be standardized to 2m using the `standardizeWindSpeed` function in EToAdjustmentMethod.
+	 */
+	windSpeed: number;
+	/** The total precipitation over the time period (in inches). */
+	precip: number;
+}
+
+const EToAdjustmentMethod: AdjustmentMethod = {
+	calculateWateringScale: calculateEToWateringScale
+};
+export default EToAdjustmentMethod;

routes/weather.spec.ts

@@ -7,6 +7,7 @@ import * as MockDate from 'mockdate';
 import { getWateringData } from './weather';
 import { GeoCoordinates, WateringData, WeatherData } from "../types";
 import { WeatherProvider } from "./weatherProviders/WeatherProvider";
+import { EToData } from "./adjustmentMethods/EToAdjustmentMethod";
 
 const expected = require( '../test/expected.json' );
 const replies = require( '../test/replies.json' );
@@ -78,16 +79,19 @@ export class MockWeatherProvider extends WeatherProvider {
     }
 
     public async getWateringData( coordinates: GeoCoordinates ): Promise< WateringData > {
-        const data = this.mockData.wateringData;
-        if ( !data.weatherProvider ) {
-            data.weatherProvider = "mock";
-        }
-
-        return data;
+        return await this.getData( "wateringData" ) as WateringData;
     }
 
     public async getWeatherData( coordinates: GeoCoordinates ): Promise< WeatherData > {
-        const data = this.mockData.weatherData;
+        return await this.getData( "weatherData" ) as WeatherData;
+    }
+
+    public async getEToData( coordinates: GeoCoordinates ): Promise< EToData > {
+        return await this.getData( "etoData" ) as EToData;
+    }
+
+    private async getData( type: "wateringData" | "weatherData" | "etoData" ) {
+        const data = this.mockData[ type ];
         if ( !data.weatherProvider ) {
             data.weatherProvider = "mock";
         }
@@ -98,5 +102,6 @@ export class MockWeatherProvider extends WeatherProvider {
 
 interface MockWeatherData {
     wateringData?: WateringData,
-    weatherData?: WeatherData
+    weatherData?: WeatherData,
+    etoData?: EToData
 }

routes/weather.ts

@@ -11,6 +11,7 @@ import { AdjustmentMethod, AdjustmentMethodResponse, AdjustmentOptions } from ".
 import ManualAdjustmentMethod from "./adjustmentMethods/ManualAdjustmentMethod";
 import ZimmermanAdjustmentMethod from "./adjustmentMethods/ZimmermanAdjustmentMethod";
 import RainDelayAdjustmentMethod from "./adjustmentMethods/RainDelayAdjustmentMethod";
+import EToAdjustmentMethod from "./adjustmentMethods/EToAdjustmentMethod";
 const weatherProvider: WeatherProvider = new ( require("./weatherProviders/" + ( process.env.WEATHER_PROVIDER || "OWM" ) ).default )();
 
 // Define regex filters to match against location
@@ -26,7 +27,8 @@ const filters = {
 const ADJUSTMENT_METHOD: { [ key: number ] : AdjustmentMethod } = {
 	0: ManualAdjustmentMethod,
 	1: ZimmermanAdjustmentMethod,
-	2: RainDelayAdjustmentMethod
+	2: RainDelayAdjustmentMethod,
+	3: EToAdjustmentMethod
 };
 
 /**
@@ -187,6 +189,11 @@ export const getWateringData = async function( req: express.Request, res: expres
 	// the string is split against a comma and the first value is selected
 	remoteAddress = remoteAddress.split( "," )[ 0 ];
 
+	if ( !adjustmentMethod ) {
+		res.send( "Error: Unknown AdjustmentMethod ID" );
+		return;
+	}
+
 	// Parse weather adjustment options
 	try {
 

routes/weatherProviders/WeatherProvider.ts

@@ -1,4 +1,5 @@
 import { GeoCoordinates, WateringData, WeatherData } from "../../types";
+import { EToData } from "../adjustmentMethods/EToAdjustmentMethod";
 
 export class WeatherProvider {
 	/**
@@ -22,4 +23,15 @@ export class WeatherProvider {
 	getWeatherData( coordinates : GeoCoordinates ): Promise< WeatherData > {
 		throw "Selected WeatherProvider does not support getWeatherData";
 	}
+
+	/**
+	 * Retrieves the data necessary for calculating potential ETo.
+	 * @param coordinates The coordinates to retrieve the data for.
+	 * @return A Promise that will be resolved with the EToData if it is successfully retrieved,
+	 * or rejected with an error message if an error occurs while retrieving the EToData or the WeatherProvider does
+	 * not support this method.
+	 */
+	getEToData( coordinates: GeoCoordinates ): Promise< EToData > {
+		throw "Selected WeatherProvider does not support getEToData";
+	};
 }

test/etoTest.json

@@ -0,0 +1,42 @@
+[
+  {
+    "description": "Badgerys Creek, AU for May 2019",
+    "source": "http://www.bom.gov.au/watl/eto/tables/nsw/badgerys_creek/badgerys_creek-201905.csv",
+    "elevation": 266,
+    "coordinates": [ -33.90, 150.73 ],
+    "startTimestamp": 1556668800,
+    "entries": [
+      {"eto":0.075,"data":{"maxTemp":76.46,"minTemp":55.04,"maxHumidity":100,"minHumidity":58,"windSpeed":2.309,"solarRadiation":2.889}},
+      {"eto":0.063,"data":{"maxTemp":77,"minTemp":56.84,"maxHumidity":100,"minHumidity":63,"windSpeed":1.707,"solarRadiation":2.406}},
+      {"eto":0.035,"data":{"maxTemp":68.36,"minTemp":56.84,"maxHumidity":100,"minHumidity":91,"windSpeed":2.309,"solarRadiation":1.186}},
+      {"eto":0.11,"data":{"maxTemp":72.86,"minTemp":58.46,"maxHumidity":100,"minHumidity":36,"windSpeed":5.254,"solarRadiation":3.375}},
+      {"eto":0.098,"data":{"maxTemp":69.44,"minTemp":48.56,"maxHumidity":96,"minHumidity":46,"windSpeed":6.324,"solarRadiation":2.947}},
+      {"eto":0.098,"data":{"maxTemp":70.16,"minTemp":47.84,"maxHumidity":97,"minHumidity":39,"windSpeed":4.551,"solarRadiation":3.8}},
+      {"eto":0.075,"data":{"maxTemp":71.42,"minTemp":39.74,"maxHumidity":100,"minHumidity":37,"windSpeed":2.259,"solarRadiation":3.767}},
+      {"eto":0.114,"data":{"maxTemp":68.36,"minTemp":41.36,"maxHumidity":99,"minHumidity":34,"windSpeed":6.676,"solarRadiation":3.6}},
+      {"eto":0.063,"data":{"maxTemp":68.72,"minTemp":36.32,"maxHumidity":99,"minHumidity":36,"windSpeed":1.673,"solarRadiation":3.65}},
+      {"eto":0.071,"data":{"maxTemp":65.66,"minTemp":41.18,"maxHumidity":100,"minHumidity":43,"windSpeed":3.999,"solarRadiation":1.878}},
+      {"eto":0.13,"data":{"maxTemp":69.08,"minTemp":42.08,"maxHumidity":78,"minHumidity":38,"windSpeed":7.88,"solarRadiation":3.608}},
+      {"eto":0.071,"data":{"maxTemp":71.6,"minTemp":38.48,"maxHumidity":99,"minHumidity":35,"windSpeed":2.158,"solarRadiation":3.606}},
+      {"eto":0.067,"data":{"maxTemp":73.04,"minTemp":38.84,"maxHumidity":100,"minHumidity":51,"windSpeed":2.326,"solarRadiation":3.469}},
+      {"eto":0.079,"data":{"maxTemp":75.74,"minTemp":43.52,"maxHumidity":100,"minHumidity":33,"windSpeed":2.242,"solarRadiation":3.542}},
+      {"eto":0.067,"data":{"maxTemp":72.68,"minTemp":44.42,"maxHumidity":100,"minHumidity":45,"windSpeed":1.991,"solarRadiation":3.506}},
+      {"eto":0.067,"data":{"maxTemp":71.6,"minTemp":44.06,"maxHumidity":100,"minHumidity":47,"windSpeed":2.326,"solarRadiation":3.464}},
+      {"eto":0.071,"data":{"maxTemp":73.94,"minTemp":43.16,"maxHumidity":100,"minHumidity":45,"windSpeed":2.393,"solarRadiation":3.411}},
+      {"eto":0.071,"data":{"maxTemp":73.4,"minTemp":45.5,"maxHumidity":100,"minHumidity":50,"windSpeed":2.56,"solarRadiation":3.417}},
+      {"eto":0.063,"data":{"maxTemp":73.22,"minTemp":51.44,"maxHumidity":100,"minHumidity":51,"windSpeed":2.342,"solarRadiation":2.783}},
+      {"eto":0.055,"data":{"maxTemp":74.12,"minTemp":46.58,"maxHumidity":100,"minHumidity":51,"windSpeed":1.69,"solarRadiation":2.706}},
+      {"eto":0.067,"data":{"maxTemp":78.44,"minTemp":44.06,"maxHumidity":100,"minHumidity":43,"windSpeed":1.723,"solarRadiation":3.289}},
+      {"eto":0.071,"data":{"maxTemp":77.36,"minTemp":47.3,"maxHumidity":100,"minHumidity":40,"windSpeed":2.125,"solarRadiation":3.267}},
+      {"eto":0.063,"data":{"maxTemp":74.48,"minTemp":53.06,"maxHumidity":100,"minHumidity":53,"windSpeed":1.991,"solarRadiation":3.175}},
+      {"eto":0.059,"data":{"maxTemp":73.58,"minTemp":44.42,"maxHumidity":100,"minHumidity":48,"windSpeed":2.008,"solarRadiation":3.108}},
+      {"eto":0.087,"data":{"maxTemp":77.9,"minTemp":42.8,"maxHumidity":100,"minHumidity":26,"windSpeed":2.828,"solarRadiation":3.272}},
+      {"eto":0.091,"data":{"maxTemp":72.68,"minTemp":44.24,"maxHumidity":92,"minHumidity":29,"windSpeed":3.865,"solarRadiation":2.747}},
+      {"eto":0.13,"data":{"maxTemp":66.02,"minTemp":39.74,"maxHumidity":82,"minHumidity":35,"windSpeed":9.905,"solarRadiation":2.425}},
+      {"eto":0.106,"data":{"maxTemp":65.66,"minTemp":37.58,"maxHumidity":69,"minHumidity":31,"windSpeed":5.739,"solarRadiation":3.211}},
+      {"eto":0.161,"data":{"maxTemp":65.48,"minTemp":47.66,"maxHumidity":52,"minHumidity":31,"windSpeed":10.859,"solarRadiation":2.997}},
+      {"eto":0.102,"data":{"maxTemp":60.08,"minTemp":36.68,"maxHumidity":70,"minHumidity":31,"windSpeed":6.743,"solarRadiation":3.172}},
+      {"eto":0.087,"data":{"maxTemp":68,"minTemp":34.34,"maxHumidity":82,"minHumidity":34,"windSpeed":4.149,"solarRadiation":3.15}}
+    ]
+  }
+]