eto_service.py

from flask import Flask, request, jsonify, current_app
from math import *
from numpy import *
from datetime import datetime
import string
import collections

SB_CONSTANT     = 4.903*1e-09
SOLAR_CONSTANT  = 0.0820 # MJ m^2/min
ALBEDO          = 0.23
kRS_CO_COASTAL  = 0.19
kRS_CO_INTERIOR = 0.16

app = Flask(__name__)
app.config.from_envvar('ETO_SETTINGS', silent=True)

Daily_ETo = collections.namedtuple('Daily_ETo', ['mm', 'inches'])

#@app.before_request
#def log_request():
   #current_app.logger.debug(request.query_string)

@app.route('/daily_eto', methods=['GET'])
def calculate_ETo():
    missing_inputs = []

    # Temperature
    # Supplied in Fahrenheit
    t_min_f = get_float_param(request, 't_min_f')
    t_max_f = get_float_param(request, 't_max_f')
    # Supplied in Celsius
    t_min_c = get_float_param(request, 't_min_c')
    t_max_c = get_float_param(request, 't_max_c')
    # Supplied in Kelvin
    t_min_k = get_float_param(request, 't_min_k')
    t_max_k = get_float_param(request, 't_max_k')

    if t_min_f and t_max_f and not t_min_c and not t_max_c:
        # Calculate C from F
        t_min_c = (t_min_f - 32.0) * (5.0 / 9.0)
        t_max_c = (t_max_f - 32.0) * (5.0 / 9.0)

    if t_min_k and t_max_k and not t_min_c and not t_max_c:
        # Calculate C from K
        t_min_c = t_min_k - 273.15
        t_max_c = t_max_k - 273.15

    if not t_min_c and not t_min_f and not t_min_k:
        missing_inputs.append("t_min_c")
    if not t_max_c and not t_max_f and not t_max_k:
        missing_inputs.append("t_max_c")

    # Relative Humidity
    rh_min = get_float_param(request, 'rh_min')
    rh_max = get_float_param(request, 'rh_max')

    if rh_min > 1:
       # Convert to percentage
       rh_min = rh_min / 100

    if rh_max > 1:
       rh_max = rh_max / 100

    if not rh_min:
        missing_inputs.append("rh_min")
    if not rh_max:
        missing_inputs.append("rh_max")

    # Region for calculating approximate solar radiation
    region = get_string_param(request, 'region', 'interior')

    # Wind speed
    u2_mph = get_float_param(request, 'u2_mph')
    u2_ms = get_float_param(request, 'u2_ms')
    u10_mph = get_float_param(request, 'u10_mph')
    u10_ms = get_float_param(request, 'u10_ms')
    u10_kts = get_float_param(request, 'u10_kts')

    if u10_kts and not u10_ms and not u2_ms:
        u10_ms = u10_kts * 0.51444444444

    if u2_mph and not u2_ms:
        # Calculate u2_ms from u2_mph
        u2_ms = u2_mph * 0.44704 # m/s
    elif u10_mph and not u10_ms:
        # Calculate u2_ms from u10_mph
        u10_ms = u10_mph * 0.44704
        u2_ms = convert_wind_speed_to_2m(u10_ms, 10)
    elif u10_ms:
        # Calculate u2_ms from u10_ms
        u2_ms = convert_wind_speed_to_2m(u10_ms, 10)

    if not u2_ms:
        missing_inputs.append("u2_ms")

    # Elevation (meters)
    elev_ft = get_float_param(request, 'elev_ft')
    elev_m = get_float_param(request, 'elev_m')

    if elev_ft and not elev_m:
        # Calculate elev_m from elev_ft
        elev_m = elev_ft * 0.3048

    if not elev_m:
        missing_inputs.append("elev_m")

    # Latitude in decimal degrees
    latitude = get_float_param(request, 'latitude')

    if not latitude:
        missing_inputs.append("latitude")

    # Julian day - defaults to current day if not provided
    julian_day = get_int_param(request, 'julian_day')
    year = get_int_param(request, 'year')
    month = get_int_param(request, 'month')
    day = get_int_param(request, 'day')

    if not julian_day and year and month and day:
       julian_day = datetime(year, month, day).timetuple().tm_yday
    else:
       julian_day = datetime.now().timetuple().tm_yday

    if len(missing_inputs) > 0:
        response = jsonify(error="Missing input parameters", missing=missing_inputs)
        response.status_code = 400
        return response
    else:
        result = calculate_daily_ETo(t_min_c, t_max_c, rh_min, rh_max, region, u2_ms, elev_m, latitude, julian_day)
        return jsonify(mm=result.mm, inches=result.inches)

def calculate_daily_ETo(t_min_c, t_max_c, rh_min, rh_max, region, u2_ms, elev_m, latitude, julian_day):
    t_min_k = t_min_c + 273.15
    t_max_k = t_max_c + 273.15
    kRS = get_kRS_coefficient_for_region(region)

    # Atmospheric pressure
    P = 101.3 * ((293 - 0.0065 * elev_m) / 293)**5.26

    # Mean temp
    t_mean = (t_min_c + t_max_c) / 2;

    # Slope of saturation vapor pressure curve
    delta = (4098.0 * (0.6108 * exp(((17.27 * t_mean) / (t_mean + 237.3))))) / (t_mean + 273.3)**2

    # Psychrometric constant
    lhv = 2.501 - 2.361e-3 * t_mean
    gamma = 1.013e-3 * P / (0.622 * lhv)

    # Mean saturation vapor pressure and actual vapor pressure
    e_tmax = 0.6108 * exp((17.27 * t_max_c) / (t_max_c + 237.3))
    e_tmin = 0.6108 * exp((17.27 * t_min_c) / (t_min_c + 237.3))
    e_sat = (e_tmax + e_tmin) / 2
    e_act = ((e_tmin * rh_max) + (e_tmax * rh_min)) / 2

    # Angles for radiation terms
    inv_dist_earth_sun_rad = 1 + 0.033 * cos(((2 * pi) / 365) * julian_day)
    solar_decl_rad = 0.409 * sin(((2 * pi / 365) * julian_day) - 1.39)

    # Latitude
    lat_rad = (pi / 180) * latitude

    # Sunset hour angle
    sunset_hr_ang_rad = arccos(-tan(lat_rad) * tan(solar_decl_rad))

    # Extraterrestrial radiation (Ra)
    Ra = ((24*60) / pi) * SOLAR_CONSTANT * inv_dist_earth_sun_rad * (sunset_hr_ang_rad * sin(lat_rad) * sin(solar_decl_rad) + cos(lat_rad) * cos(solar_decl_rad) * sin(sunset_hr_ang_rad))

    # Solar radiation (Longwave)
    Rs = kRS * sqrt(t_max_c - t_min_c) * Ra

    # Clear sky radiation
    Rso = (0.75 + 2e-05 * elev_m) * Ra

    # Net shortwave radiation
    Rns = (1 - ALBEDO) * Rs

    # Net outgoing long wave solar radiation
    Rnl = SB_CONSTANT * ((t_max_k**4 + t_min_k**4) / 2) * (0.34 - (0.14 * sqrt(e_act))) * (1.35 * (Rs / Rso) - 0.35)

    # Net radiation
    Rn = Rns - Rnl
    Rng = 0.408 * Rn # mm

    # Radiation term
    DT = delta / (delta + gamma * (1 + 0.34 * u2_ms))
    ETrad = DT * Rng

    # Wind term
    PT = gamma / (delta + (gamma * (1 + 0.34 * u2_ms)))
    TT = (900 / (t_mean + 273)) * u2_ms
    ETwind = PT * TT * (e_sat - e_act)

   # Final ref ET
    ETo_mm = round(ETwind + ETrad, 4) # mm per day
    ETo_in = round(ETo_mm * 0.0393701, 4)

    return Daily_ETo(ETo_mm, ETo_in)

def convert_wind_speed_to_2m(uz_ms, z):
   # log below = ln
   return uz_ms * (4.87 / log(67.8 * z - 5.42))

def get_kRS_coefficient_for_region(region):
    if string.lower(region) == 'coastal':
        return kRS_CO_COASTAL
    else:
        return kRS_CO_INTERIOR

def get_float_param(request, key):
    return request.args.get(key, type=float)

def get_string_param(request, key):
    return request.args.get(key, type=str)

def get_string_param(request, key, default):
    return request.args.get(key, default, type=str)

def get_int_param(request, key):
    return request.args.get(key, type=int)

if __name__ == '__main__':
    app.run(debug=True, host='0.0.0.0')