From 3e0141cdac7b62435ab9c8948c2902bd82bc76d5 Mon Sep 17 00:00:00 2001 From: Matthew Oslan Date: Tue, 11 Jun 2019 19:10:32 -0400 Subject: [PATCH] Add ETo AdjustmentMethod --- .../EToAdjustmentMethod.spec.ts | 43 ++++ .../adjustmentMethods/EToAdjustmentMethod.ts | 224 ++++++++++++++++++ routes/weather.spec.ts | 21 +- routes/weatherProviders/WeatherProvider.ts | 12 + test/etoTest.json | 42 ++++ 5 files changed, 334 insertions(+), 8 deletions(-) create mode 100644 routes/adjustmentMethods/EToAdjustmentMethod.spec.ts create mode 100644 routes/adjustmentMethods/EToAdjustmentMethod.ts create mode 100644 test/etoTest.json diff --git a/routes/adjustmentMethods/EToAdjustmentMethod.spec.ts b/routes/adjustmentMethods/EToAdjustmentMethod.spec.ts new file mode 100644 index 0000000..4698cc6 --- /dev/null +++ b/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, + timestamp: 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 + }[]; +} diff --git a/routes/adjustmentMethods/EToAdjustmentMethod.ts b/routes/adjustmentMethods/EToAdjustmentMethod.ts new file mode 100644 index 0000000..47e0980 --- /dev/null +++ b/routes/adjustmentMethods/EToAdjustmentMethod.ts @@ -0,0 +1,224 @@ +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 + } + } +} + +/* 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; + + 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.timestamp ).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 solarRadiation = etoData.solarRadiation; + + 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 ); +} + +// The time at which the formula for clear sky isolation will start/stop yielding a non-negative result. +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 megajoules 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 megajoules. + const clearSkyIsolation = ( 990 * Math.sin( solarElevationAngle ) - 30 ) * 0.0036 * 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. */ + timestamp: 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 megajoules 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; diff --git a/routes/weather.spec.ts b/routes/weather.spec.ts index 36e92e5..de6bea1 100644 --- a/routes/weather.spec.ts +++ b/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 } diff --git a/routes/weatherProviders/WeatherProvider.ts b/routes/weatherProviders/WeatherProvider.ts index 4826df2..e8ea210 100644 --- a/routes/weatherProviders/WeatherProvider.ts +++ b/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"; + }; } diff --git a/test/etoTest.json b/test/etoTest.json new file mode 100644 index 0000000..7519b82 --- /dev/null +++ b/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":10.4}}, + {"eto":0.063,"data":{"maxTemp":77,"minTemp":56.84,"maxHumidity":100,"minHumidity":63,"windSpeed":1.707,"solarRadiation":8.66}}, + {"eto":0.035,"data":{"maxTemp":68.36,"minTemp":56.84,"maxHumidity":100,"minHumidity":91,"windSpeed":2.309,"solarRadiation":4.27}}, + {"eto":0.11,"data":{"maxTemp":72.86,"minTemp":58.46,"maxHumidity":100,"minHumidity":36,"windSpeed":5.254,"solarRadiation":12.15}}, + {"eto":0.098,"data":{"maxTemp":69.44,"minTemp":48.56,"maxHumidity":96,"minHumidity":46,"windSpeed":6.324,"solarRadiation":10.61}}, + {"eto":0.098,"data":{"maxTemp":70.16,"minTemp":47.84,"maxHumidity":97,"minHumidity":39,"windSpeed":4.551,"solarRadiation":13.68}}, + {"eto":0.075,"data":{"maxTemp":71.42,"minTemp":39.74,"maxHumidity":100,"minHumidity":37,"windSpeed":2.259,"solarRadiation":13.56}}, + {"eto":0.114,"data":{"maxTemp":68.36,"minTemp":41.36,"maxHumidity":99,"minHumidity":34,"windSpeed":6.676,"solarRadiation":12.96}}, + {"eto":0.063,"data":{"maxTemp":68.72,"minTemp":36.32,"maxHumidity":99,"minHumidity":36,"windSpeed":1.673,"solarRadiation":13.14}}, + {"eto":0.071,"data":{"maxTemp":65.66,"minTemp":41.18,"maxHumidity":100,"minHumidity":43,"windSpeed":3.999,"solarRadiation":6.76}}, + {"eto":0.13,"data":{"maxTemp":69.08,"minTemp":42.08,"maxHumidity":78,"minHumidity":38,"windSpeed":7.88,"solarRadiation":12.99}}, + {"eto":0.071,"data":{"maxTemp":71.6,"minTemp":38.48,"maxHumidity":99,"minHumidity":35,"windSpeed":2.158,"solarRadiation":12.98}}, + {"eto":0.067,"data":{"maxTemp":73.04,"minTemp":38.84,"maxHumidity":100,"minHumidity":51,"windSpeed":2.326,"solarRadiation":12.49}}, + {"eto":0.079,"data":{"maxTemp":75.74,"minTemp":43.52,"maxHumidity":100,"minHumidity":33,"windSpeed":2.242,"solarRadiation":12.75}}, + {"eto":0.067,"data":{"maxTemp":72.68,"minTemp":44.42,"maxHumidity":100,"minHumidity":45,"windSpeed":1.991,"solarRadiation":12.62}}, + {"eto":0.067,"data":{"maxTemp":71.6,"minTemp":44.06,"maxHumidity":100,"minHumidity":47,"windSpeed":2.326,"solarRadiation":12.47}}, + {"eto":0.071,"data":{"maxTemp":73.94,"minTemp":43.16,"maxHumidity":100,"minHumidity":45,"windSpeed":2.393,"solarRadiation":12.28}}, + {"eto":0.071,"data":{"maxTemp":73.4,"minTemp":45.5,"maxHumidity":100,"minHumidity":50,"windSpeed":2.56,"solarRadiation":12.3}}, + {"eto":0.063,"data":{"maxTemp":73.22,"minTemp":51.44,"maxHumidity":100,"minHumidity":51,"windSpeed":2.342,"solarRadiation":10.02}}, + {"eto":0.055,"data":{"maxTemp":74.12,"minTemp":46.58,"maxHumidity":100,"minHumidity":51,"windSpeed":1.69,"solarRadiation":9.74}}, + {"eto":0.067,"data":{"maxTemp":78.44,"minTemp":44.06,"maxHumidity":100,"minHumidity":43,"windSpeed":1.723,"solarRadiation":11.84}}, + {"eto":0.071,"data":{"maxTemp":77.36,"minTemp":47.3,"maxHumidity":100,"minHumidity":40,"windSpeed":2.125,"solarRadiation":11.76}}, + {"eto":0.063,"data":{"maxTemp":74.48,"minTemp":53.06,"maxHumidity":100,"minHumidity":53,"windSpeed":1.991,"solarRadiation":11.43}}, + {"eto":0.059,"data":{"maxTemp":73.58,"minTemp":44.42,"maxHumidity":100,"minHumidity":48,"windSpeed":2.008,"solarRadiation":11.19}}, + {"eto":0.087,"data":{"maxTemp":77.9,"minTemp":42.8,"maxHumidity":100,"minHumidity":26,"windSpeed":2.828,"solarRadiation":11.78}}, + {"eto":0.091,"data":{"maxTemp":72.68,"minTemp":44.24,"maxHumidity":92,"minHumidity":29,"windSpeed":3.865,"solarRadiation":9.89}}, + {"eto":0.13,"data":{"maxTemp":66.02,"minTemp":39.74,"maxHumidity":82,"minHumidity":35,"windSpeed":9.905,"solarRadiation":8.73}}, + {"eto":0.106,"data":{"maxTemp":65.66,"minTemp":37.58,"maxHumidity":69,"minHumidity":31,"windSpeed":5.739,"solarRadiation":11.56}}, + {"eto":0.161,"data":{"maxTemp":65.48,"minTemp":47.66,"maxHumidity":52,"minHumidity":31,"windSpeed":10.859,"solarRadiation":10.79}}, + {"eto":0.102,"data":{"maxTemp":60.08,"minTemp":36.68,"maxHumidity":70,"minHumidity":31,"windSpeed":6.743,"solarRadiation":11.42}}, + {"eto":0.087,"data":{"maxTemp":68,"minTemp":34.34,"maxHumidity":82,"minHumidity":34,"windSpeed":4.149,"solarRadiation":11.34}} + ] + } +]