Created
October 16, 2021 02:26
-
-
Save electricshaman/a108619bc28965730485011f9b0d404f to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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') |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment