diff --git a/code/common/calc_eto.pl b/code/common/calc_eto.pl new file mode 100644 index 000000000..b27301b10 --- /dev/null +++ b/code/common/calc_eto.pl @@ -0,0 +1,1089 @@ +# Category = Irrigation + +# June 2016 +# v1.2 +# - Added in minimum time calculation + +#@ This module allows MisterHouse to calculate daily EvapoTranspiration based on a +#@ Data feed from Weatherunderground (WU). To use it you need to sign up for a weatherundeground key +#@ The free developer account provides a low-volume hobbyist service that will work for these calcs +#@ A location is also required. Best is a lat/long pair. +#@ By default wuData is written to $Data_Dir/wuData and the eto logs are written to $Data_Dir/eto +#@ +#@ The ET programs can be automatically uploaded to an OpenSprinkler. (need v1.1 of the lib) + +########################################################################################################### +## Credits ## +########################################################################################################### +## portions of code provided by Zimmerman method used by OpenSprinkler ## +## portions of code provided/edited by Ray and Samer of OpenSprinkler ## +## eto library provided by Mark Richards ## +## Original python code provided by Shawn Harte 2014 no copyright reserved ## +## python cleanup and first pass my Neil Cherry ## +## Code was used with utmost respect to the original authors, your efforts have prevented the ## +## re-invention of the wheel ## +########################################################################################################### + +#The script accounts for wind, freezing conditions, and current/recent +#rainfall when considering start and run times. It will avoid watering +#during midday, unless early morning winds prevent earlier start times. +#The starts are serialized so no odd overlaps should occur. Mornings +#are preferred to evenings to allow for the best use of water and +#absorption without causing mold and fungus growth by leaving grass wet +#overnight. The script is commented quite heavily, so that anyone can +#edit or use it to their liking. Please be mindful that other authors +#work was used or modified when the code seemed generalized enough that +#I shouldn’t be stepping on toes. Please do not pester the original +#author if something doesn’t work for you, as they will probably have +#enough on their own plate with their own original works. +# +#Everything is done based off your latitude and longitude, however, the +#script can find the info when provided with a city/state or country, a +#US Zip Code, or a PWS ID. + +#Usage: +#Create a config_parm{eto_zone_1mm} with the number of seconds to distribute 1mm of water in the zone. +#Create a config_parm{eto_zone_crop} with 1's and 0's (1=grass, 0=shrubs/garden) +#Find your closest weatherunderground location and store it in config_parms{eto_location} +#Get your wu api key and store it in config_parms{wu_key} + +#TODO +# - the safefloat and safeint subs are from python. don't know if they're needed +# - if run after sunrise, then use the sunset times and max 2 cycles (line 677) +# - populate yesterday's rain value in the RRD + +#VERIFY +# - line 430 sub getConditionsData chkcond array isn't checked yet +# - line 63 Use use JSON qw(decode_json) instead of JSON::XS +# - line 360 test timezone subroutine, confirm that it actually works +# - line 610 read in multiple water times for overall aggregate +# - line 711 when multiple times are scheduled, only one entry was written to the logs. + +use eto; +use LWP::UserAgent; +use HTTP::Request::Common; +#use JSON::XS; +use JSON qw(decode_json); +use List::Util qw(min max sum); +#use Data::Dumper; +use Time::Local; +use Date::Calc qw(Day_of_Year); +my $debug = 0; +my $msg_string; + +$p_wu_forecast = new Process_Item qq[get_url --quiet "http://api.wunderground.com/api/$config_parms{wu_key}/astronomy/yesterday/conditions/forecast/q/$config_parms{eto_location}.json" "$config_parms{data_dir}/wuData/wu_data.json"]; +$v_get_eto = new Voice_Cmd("Update ETO Programs"); +$t_wu_forecast_timer = new Timer; + +my $eto_data_dir = $config_parms{data_dir} . "/eto"; +$eto_data_dir = $config_parms{eto_data_dir} if (defined $config_parms{eto_data_dir}); + +my $eto_calc_time = "3:00 AM"; +$eto_calc_time = $config_parms{eto_calc_time} if (defined $config_parms{eto_calc_time}); +$eto_calc_time = " " . $eto_calc_time if ($eto_calc_time =~ m/^\d:\d\d\s/); #time_now has a space in front if only a single digit hour + +my $eto_retries = 3; +$eto_retries = $config_parms{eto_retries} if (defined $config_parms{eto_retries}); +my $eto_retries_today; + +$config_parms{eto_rainfallsatpoint} = 25 unless (defined $config_parms{eto_rainfallsatpoint}); +$config_parms{eto_minmax} = "5,15" unless (defined $config_parms{eto_minmax}); + +my $eto_ready; + +if ($Startup or $Reload) { + $eto_ready = 1; + print_log "[calc_eto] v1.2 Startup. Checking Configuration..."; + mkdir "$eto_data_dir" unless ( -d "$eto_data_dir" ); + mkdir "$eto_data_dir/ET" unless ( -d "$eto_data_dir/ET" ); + mkdir "$eto_data_dir/logs" unless ( -d "$eto_data_dir/logs" ); + mkdir "$config_parms{data_dir}/wuData" unless ( -d "$config_parms{data_dir}/wuData" ); + mkdir "$eto_data_dir/weatherprograms" unless ( -d "$eto_data_dir/weatherprograms" ); + + if (defined $config_parms{eto_location}) { + print_log "[calc_eto] Location : $config_parms{eto_location}"; + } else { + print_log "[calc_eto] ERROR! eto_location undefined!!"; + $eto_ready = 0; + } + if (defined $config_parms{eto_zone_1mm}) { + print_log "[calc_eto] 1mm zone runtimes (in seconds) : $config_parms{eto_zone_1mm}"; + } else { + print_log "[calc_eto] ERROR! eto_zone_1mm undefined!!"; + $eto_ready = 0; + } + if (defined $config_parms{eto_zone_crop}) { + print_log "[calc_eto] 1mm crop definitions : $config_parms{eto_zone_crop}"; + } else { + print_log "[calc_eto] ERROR! eto_zone_crop undefined!!"; + $eto_ready = 0; + } + unless (defined $config_parms{wu_key}) { + print_log "[calc_eto] ERROR! wu key undefined!!"; + $eto_ready = 0; + } + if (defined $config_parms{eto_irrigation}) { + print_log "[calc_eto] $config_parms{eto_irrigation} set as programmable irrigation system"; + } else { + print_log "[calc_eto] WARNING! no sprinkler system defined!"; + } + if ($eto_ready) { + print_log "[calc_eto] Configuration good. ETo Calcuations Ready"; + print_log "[calc_eto] Will email results to $config_parms{eto_email}" if (defined $config_parms{eto_email}); + } else { + print_log "[calc_eto] ERROR! ETo configuration problem. ETo will not calcuate"; + } +} + +if ((said $v_get_eto) or ($New_Minute and ($Time_Now eq $eto_calc_time))) { + if ($eto_ready) { + print_log "[calc_eto] Starting Daily ETO Calculation Process..."; + $eto_retries_today = 0; + start $p_wu_forecast; + } else { + print_log "[calc_eto] ERROR! ETo configuration problem. ETo will not calcuate"; + } +} + +if (done_now $p_wu_forecast) { + my $write_secs = time() - (stat("$config_parms{data_dir}/wuData/wu_data.json"))[9]; + if ($write_secs > 300) { + print_log "[calc_eto] Stale Data, not calculating ETo. WU Data written $write_secs ago..."; + } else { + my $program_data = &calc_eto_runtimes($eto_data_dir,"file",$config_parms{eto_location}, "$config_parms{data_dir}/wuData/wu_data.json"); + if ($program_data) { + if (defined $config_parms{eto_irrigation}) { + if (lc $config_parms{eto_irrigation} eq "opensprinkler") { + my $os_program = &get_object_by_name($config_parms{eto_opensprinkler_program}); + my ($run_times,$run_seconds) = $program_data =~ /\[\[(.*)\],\[(.*)\]\]/; + print_log "[calc_eto] Loading values $run_times,$run_seconds into program $config_parms{eto_opensprinkler_program}"; + $os_program->set_program($Day,$run_times,$run_seconds); + } #elsif (other sprinkler system...) + } + } else { + if ($eto_retries_today < $eto_retries) { + $eto_retries_today++; + print_log "[calc_eto] WARNING! bad program data, retry attempt $eto_retries_today"; + set $t_wu_forecast_timer 600; + #start $p_wu_forecast; + } else { + print_log "[calc_eto] ERROR! retry max $eto_retries reaches. Aborting calculation attempt"; + } + } + } +} + + +if (expired $t_wu_forecast_timer) { + start $p_wu_forecast; +} + +#-------------------------------------------------------------------------------------------------------------------------------# +# Mapping of conditions to a level of shading. +# Since these are for sprinklers any hint of snow will be considered total cover (10) +# Don't worry about wet conditions like fog these are accounted for below we are only concerned with how much sunlight is blocked at ground level + + our $conditions = { + 'Clear' => 0, + 'Partial Fog' => 2, + 'Patches of Fog' => 2, + 'Haze' => 2, + 'Shallow Fog' => 3, + 'Scattered Clouds' => 4, + 'Unknown' => 5, + 'Fog' => 5, + 'Partly Cloudy' => 5, + 'Mostly Cloudy' => 8, + 'Mist' => 10, + 'Light Drizzle' => 10, + 'Light Freezing Drizzle' => 10, + 'Light Freezing Rain' => 10, + 'Light Freezing Fog' => 5, + 'Light Ice Pellets' => 10, + 'Light Rain' => 10, + 'Light Rain Showers' => 10, + 'Light Snow' => 10, + 'Light Snow Grains' => 10, + 'Light Snow Showers' => 10, + 'Light Thunderstorms and Rain' => 10, + 'Low Drifting Snow' => 10, + 'Rain' => 10, + 'Rain Showers' => 10, + 'Snow' => 10, + 'Snow Showers' => 10, + 'Thunderstorm' => 10, + 'Thunderstorms and Rain' => 10, + 'Blowing Snow' => 10, + 'Chance of Snow' => 10, + 'Freezing Rain' => 10, + 'Unknown Precipitation' => 10, + 'Overcast' => 10, + }; + +# List of precipitation conditions we don't want to water in, the conditions will be checked to see if they contain these phrases. + + our $chkcond = { + 'flurries' => 1, + 'rain' => 1, + 'sleet' => 1, + 'snow' => 1, + 'storm' => 1, + 'hail' => 1, + 'ice' => 1, + 'squall' => 1, + 'precip' => 1, + 'funnel' => 1, + 'drizzle' => 1, + 'mist' => 1, + 'freezing' => 1, + }; + + + +# +################################################################################ +# -[ Functions ]---------------------------------------------------------------- + +# define safe functions for variable conversion, preventing errors with NaN and Null as string values +# 's'=value to convert 'dv'=value to default to on error make sure this is a legal float or integer value + +#just stub for testing, have to fix up the floats and int +sub safe_float { + my ($arg,$val) = @_; +# $val = "0.0" unless ($val); +# $arg = $val unless ($arg); + return ($arg); +} + +sub safe_int { + my ($arg,$val) = @_; +# $val = "0" unless ($val); +# $arg = $val unless ($arg); + return ($arg); +} + +sub isInt { + my ($arg) = @_; + return ($arg - int($arg))? 0 : 1; +} + +sub isFloat { + my ($arg) = @_; + return 1; + #return ($arg - int($arg))? 1 : 0; +} + +sub round { + my ($number, $places) = @_; + my $sign = ($number < 0) ? '-' : ''; + my $abs = abs($number); + + if($places < 0) { + print_log "[calc_eto] ERROR! rounding to $places"; + return $number; + } else { + my $p10 = 10**$places; + return $sign . int($abs*$p10 + 0.5)/$p10; + } +} + +sub findwuLocation { + my ($loc) = @_; + my ($whttyp, $ploc, $noData, $tzone, $lat, $lon); + my $ua = new LWP::UserAgent( keep_alive => 1 ); + + my $request = HTTP::Request->new( GET => "http://autocomplete.wunderground.com/aq?format=json&query=$loc" ); + my $responseObj = $ua->request($request); + my $data; +# eval { $data = JSON::XS->new->decode( $responseObj->content ); }; + eval { $data = decode_json( $responseObj->content ); }; + my $responseCode = $responseObj->code; + my $isSuccessResponse = $responseCode < 400; + if ( $isSuccessResponse and defined $data->{RESULTS}) { + my $chk = $data->{RESULTS}->[0]->{ll}; # # ll has lat and lon in one spot no matter how we search + if ($chk) { + my @ll = split(' ', $chk); + if (scalar (@ll) == 2 and isFloat($ll[0]) and isFloat($ll[1])) { + $lat = $ll[0]; + $lon = $ll[1]; + } + } + + $chk = $data->{RESULTS}->[0]->{tz}; + if ($chk) { + $tzone = $chk; + } else { + my $chk2 = $data->{RESULTS}->[0]->{tz_long}; + if ($chk2) { + $tzone = $chk2; + } else { + $tzone = "None"; + } + } + + $chk = $data->{RESULTS}->[0]->{name}; # # this is great for showing a pretty name for the location + if ($chk) { + $ploc = $chk; + } + + $chk = $data->{RESULTS}->[0]->{type}; + if ($chk) { + $whttyp = $chk; + } + + } else { + $noData = 1; + $lat = "None"; + $lon = "None"; + $tzone = "None"; + $ploc = "None"; + $whttyp = "None"; + } + return ($whttyp, $ploc, $noData, $tzone, $lat, $lon); +} + +sub getwuData { + my ($loc, $key) = @_; + my $tloc = split(',', $loc); + #return if ($key == '' or (scalar ($tloc) < 2)); + my $ua = new LWP::UserAgent( keep_alive => 1 ); + + my $request = HTTP::Request->new( GET => "http://api.wunderground.com/api/$key/astronomy/yesterday/conditions/forecast/q/$loc.json" ); + + my $responseObj = $ua->request($request); + my $data; +# eval { $data = JSON::XS->new->decode( $responseObj->content ); }; + eval { $data = decode_json( $responseObj->content ); }; + if ($@) { + print_log "[calc_eto] ERROR problem parsing json from web call"; + } + my $responseCode = $responseObj->code; + my $isSuccessResponse = $responseCode < 400; + print "code=$responseCode\n" if ($debug); + + return ($data); + +} + +sub getwuDataTZOffset { + my ($data,$tzone) = @_; +#HP TODO - I'm not sure if this works as expected + if ($tzone eq "None" or $tzone eq "") { + $tzone = $data->{current_observation}->{local_tz_long}; + } + my $tdelta; + + if ($tzone) { + my @tnow = localtime(time); + $tdelta = timegm(@tnow) - timelocal(@tnow); +# tdelta = tnow.astimezone(pytz.timezone(tz)).utcoffset() + } + if ($tdelta) { + return ({'t' => ($tdelta / 900 + 48), 'gmt' => ($tdelta / 3600)}); + } else { + return ({'t' => "None", 'gmt' => "None"}); + } +} + +# Calculate an adjustment based on predicted rainfall +# Rain forecast should lessen current watering and reduce rain water runoff, making the best use of rain. +# returns tadjust (???) +sub getForecastData { + my ($data) = @_; +#HP TODO - I don't know why the python wanted to create a bunch of arrays (mm, cor, wfc). It seems like +#HP TODO - just the end result is needed + if (@{$data}) { + + my $fadjust = 0; + for (my $day = 1; $day < scalar (@{$data}); $day++) { + my $mm = $data->[$day]->{qpf_allday}->{mm}; + my $cor = $data->[$day]->{pop}; + my $wfc = 1 / $day ** 2; #HP I assume this is to modify that further days out are more volatile? + $fadjust += safe_float($mm, -1) * (safe_float($cor, -1) / 100) * safe_float($wfc, -1); + } + return $fadjust; + } + return -1; +} + +# Grab the sunrise and sunset times in minutes from midnight +sub getAstronomyData { + my ($data) = @_; + + if (not $data) { + return ({"rise" => -1, "set" => -1}); + } + + my $rHour = safe_int($data->{'sunrise'}->{'hour'}, 6); + my $rMin = safe_int($data->{'sunrise'}->{'minute'}); + my $sHour = safe_int($data->{'sunset'}->{'hour'}, 18); + my $sMin = safe_int($data->{'sunset'}->{'minute'}); + if ($rHour, $rMin, $sHour, $sMin) { + return ({"rise" => $rHour * 60 + $rMin, "set" => $sHour * 60 + $sMin}); + } else { + return ({"rise" => -1, "set" => -1}); + } +} + +# Let's check the current weather and make sure the wind is calm enough, it's not raining, and the temp is above freezing +# We will also look at what the rest of the day is supposed to look like, we want to stop watering if it is going to rain, +# or if the temperature will drop below freezing, as it would be bad for the pipes to contain water in these conditions. +# Windspeed for the rest of the day is used to determine best low wind watering time. + +sub getConditionsData { + my ($current, $predicted,$conditions) = @_; + + my $nowater = 1; + my $whynot = 'Unknown'; + unless ($current and $predicted) { + return (0, 1, 'No conditions data'); + } + + my $cWeather = ""; + $cWeather = safe_float($conditions->{$current->{weather}}, 5); + + unless (defined $conditions->{$current->{weather}}) { + # check if any of the chkcond words exist in the $current-{weather} + + my $badcond = 0; + foreach my $chkword (split(' ',lc $current->{weather})) { + $badcond = 1 if (defined $chkcond->{$chkword}); + } + +# if (defined $conditions->{$current->{weather}} ) { + if ($badcond) { + $cWeather = 10; + } else { + print_log '[calc_eto] INFO Cound not find current conditions ' . $current->{weather}; + $cWeather = 5; + } + } + + my $cWind = &eto::wind_speed_2m(safe_float($current->{wind_kph}), 10); + my $cTemp = safe_float($current->{temp_c}, 20); + + # current rain will only be used to adjust watering right before the start time + + my $cmm = safe_float($current->{precip_today_metric}); # Today's predicted rain (mm) + my $pWind = &eto::wind_speed_2m(safe_float($predicted->{avewind}->{kph}), 10); # Today's predicted wind (kph) + my $pLowTemp = safe_float($predicted->{low}->{celsius}); # Today's predicted low (C) + my $pCoR = safe_float($predicted->{pop}) / 100; # Today's predicted POP (%) (Probability of Precipitation) + my $pmm = safe_float($predicted->{qpf_allday}->{mm}); # Today's predicted QFP (mm) (Quantitative Precipitation Forecast) + # + + # Let's check to see if it's raining, windy, or freezing. Since watering is based on yesterday's data + # we will see how much it rained today and how much it might rain later today. This should + # help reduce excess watering, without stopping water when little rain is forecast. + + $nowater = 0; + $whynot = ''; + + # Its precipitating +#HP TODO - this triggered on 'Clear'? + if ($cWeather == 10 and lc $current->{weather} ne 'overcast') { + $nowater = 1; + $whynot .= 'precipitation (' . $current->{weather} . ') '; + } + + # Too windy + if ($cWind > $pWind and $pWind > 6 or $cWind > 8 ) { + $nowater = 1; + $whynot .= 'wind (' . round($cWind, 2) . ' kph) '; + } + + # Too cold + if ($cTemp < 4.5 or $pLowTemp < 1) { + $nowater = 1; + $whynot .= 'cold (current ' . round($cTemp, 2) . ' C / predicted ' . round($pLowTemp,2) . ' C) '; + } + + $cmm += $pmm * $pCoR if ($pCoR); + +#HP TODO - Don't know where this except comes from +#HP except: +#HP print 'we had a problem and just decided to water anyway' +#HP nowater = 0 + # +#print "[$cmm,$nowater,$whynot]\n"; + return ($cmm, $nowater, $whynot); +} + +sub sun_block { +# Difference from Python script. If there are multiple forecasts for a given hour (ie overcast and scattered clouds), then it will +# take the last entry for calculating cover. Could average it, but really the difference isn't that huge I don't think. + my ($wuData, $sunrise, $sunset, $conditions) = @_; + my $sh = 0; + my $previousCloudCover = 0; + + for (my $hour = int($sunrise / 60); $hour < int($sunset / 60 + 1); $hour++) { + # Set a default value so we know we found missing data and can handle the gaps + my $cloudCover = -1; + + # Now let's find the data for each hour there are more periods than hours so only grab the first + #in range(len(wuData['history']['observations'])): + for (my $period = 0; $period < scalar (@{$wuData->{history}->{observations}}); $period++) { + if (safe_int($wuData->{history}->{observations}->[$period]->{date}->{hour}, -1) == $hour ) { + if ($wuData->{history}->{observations}->[$period]->{conds}) { + print "[$hour," . $wuData->{history}->{observations}->[$period]->{conds} . "," . $conditions->{$wuData->{history}->{observations}->[$period]->{conds}} . "]\n" if ($debug); + $cloudCover = safe_float($conditions->{$wuData->{history}->{observations}->[$period]->{conds}}, 5) / 10; + unless (defined $cloudCover) { + $cloudCover = 10; + print_log '[calc_eto] INFO Sun Block Condition not found ' . $wuData->{history}->{observations}->[$period]->{conds}; + } + } + } + } + + # Found nothing, let's assume it was the same as last hour + $cloudCover = $previousCloudCover if ($cloudCover == -1); + print "[$hour,$cloudCover]\n" if ($debug); + # + + $previousCloudCover = $cloudCover; + + # Got something now? let's check + $sh += 1 - $cloudCover if ($cloudCover != -1); + print "total $sh $cloudCover\n" if ($debug); + + } + return ($sh); +} + +# We need to know how much it rained yesterday and how much we watered versus how much we required +sub mmFromLogs { + my ($_1mmProg,$logsPath,$ETPath) = @_; + + + my $prevLogFname = int((time - (time % 86400) -1) / 86400); + #look that the $prevLogFname exists for both logs and ET. if not look for up to 7 days back for a + #day that they both exist + my $tmpLogFname = $prevLogFname; + my $fnamefound = 0; + for (my $fx = $tmpLogFname; $fx > $tmpLogFname - 7; $fx--) { + print "***** Looking for $fx\n" if ($debug); + if ((-e "$logsPath/$fx") and (-e "$ETPath/$fx")) { + print "log file found $fx\n" if ($debug); + $prevLogFname = $fx; + $fnamefound = 1; + last; + } + } + print_log "[calc.eto] WARNING! Couldn't find log/ET files less than 7 days old" unless ($fnamefound); + + my $nStations = scalar( @{$_1mmProg->{mmTime}} ); + + my @ydur = (-1) x $nStations; + my @ymm = (-1) x $nStations; + + my @yET = (0,0); # Yesterday's Evap (evapotranspiration, moisture losses mm/day) + my @tET = (0,0); + + + # -[ Logs ]----------------------------------------------------------------- + my @logs = (); + + if (open (FILE, "$logsPath/$prevLogFname")) { + my $d_logs = ; + my $t_logs; +# eval { $t_logs = JSON::XS->new->decode($d_logs) }; + eval { $t_logs = decode_json( $d_logs ); }; + @logs = @$t_logs; + close (FILE); + } else { + print_log "[calc_eto] WARNING Can't open file $logsPath/$prevLogFname!"; + close (FILE); + } + + ### the original code first looked for yesterday's log file and used that + ### filename to get the json from ETPath and LogsPath. + ### I simply check for yesterday's files and if I get an exception I create + ### default vaules of 0 (in the appropriate array format) + ### We now look 7 days back to find the last file. If nothing exists for 7 days, then use 0's + + # -[ ET ]------------------------------------------------------------------- + if (open (FILE, "$ETPath/$prevLogFname")) { + my $d_yET = ; + my $t_yET; +# eval { $t_yET = JSON::XS->new->decode($d_yET) }; + eval { $t_yET = decode_json( $d_yET ) }; + @yET = @$t_yET; + close (FILE); + } else { + print_log "[calc_eto] WARNING Can't open file $ETPath/$prevLogFname!"; + close (FILE); + } + + # add all the run times together (a zone can have up to 4 daily runtimes) to get the overall amount of water + for (my $x = 0; $x < scalar (@logs); $x++) { + $ydur[$logs[$x][1]] += $logs[$x][2]; + print "[logs[$x][2] = " . $logs[$x][2] . " ydur[$logs[$x][1]] = " . $ydur[$logs[$x][1]] . "]\n"; + } + + for (my $x =0; $x < $nStations; $x++) { + if ($_1mmProg->{mmTime}[$x]) { + # 'mmTime': [15, 16, 20, 10, 30, 30] sec/mm + # 'crop': [1, 1, 1, 1, 0, 0] 1 = grasses, 0 = shrubs + #ymm[x] = round( (safe_float(yET[safe_int(_1mmProg['crop'][x])])) - (ydur[x]/safe_float(_1mmProg['mmTime'][x])), 4) * (-1) + # Rewritten to make it readable (nothing more) + my $yesterdaysET = safe_float($yET[safe_int($_1mmProg->{crop}[$x])]); # in seconds + my $yesterdaysDuration = $ydur[$x]; # in mm + my $mmProg = safe_float($_1mmProg->{mmTime}[$x]); # in seconds/mm + # ymm = yET - (ydur / mmTime) // mm - (sec / sec/mm) Units look correct! + $ymm[$x] = round(($yesterdaysET - ($yesterdaysDuration/$mmProg) ),4) * (-1); + $tET[int($_1mmProg->{crop}[$x])] = $ymm[$x]; + print "[$x yesterdaysET=$yesterdaysET yesterdaysDuration=$yesterdaysDuration mmProg=$mmProg ymm[$x] = ". $ymm[$x] . " tET[" . int($_1mmProg->{crop}[$x]) . "] = " . $ymm[$x] ."]\n" if ($debug); + print "E: $x $ymm[$x] = ( " . $yET[$_1mmProg->{crop}[$x]] . " ) - ( $ydur[$x] / " . $_1mmProg->{mmTime}[$x] . " ) * -1\n" if ($debug); + print "E: $x _1mmProg['crop'][$x] = ". $_1mmProg->{crop}[$x] ."\n" if ($debug); + print "E: $x tET[" . int($_1mmProg->{crop}[$x]) . "] = " . $tET[int($_1mmProg->{crop}[$x])] . "\n" if ($debug); + } else { + $ymm[$x] = 0; + } + } + + print "E: Done - mmFromLogs\n" if ($debug); + return (\@ymm, \@tET); +} + + +sub writeResults { + my ($ETG,$ETS, $sun,$todayRain,$tadjust,$noWater,$logsPath,$ETPath,$WPPath) = @_; + + my @ET = ($ETG,$ETS); + my $msg; + my $pid = 2; #for legacy purposes, can probably remove it, but for now want the log files + #to be similar to keep the file structure the same to validate against python + my $data_1mm; + +# Get 1mm & crop data from config_parms + @{$data_1mm->{mmTime}} = split (/,/,$config_parms{eto_zone_1mm}); + @{$data_1mm->{crop}} = split (/,/,$config_parms{eto_zone_crop}); + + + my @minmax; + if (defined $config_parms{eto_minmax}) { + @minmax = split (/,/,$config_parms{eto_minmax}); + } else { + @minmax = (5,15); + } + + my $fname = int((time) / 86400); + + my $minRunmm = 5; + $minRunmm = min @minmax if (scalar (@minmax) > 0) and ((min @minmax) >= 0); + my $maxRunmm = 15; + $maxRunmm = max @minmax if (scalar (@minmax) > 1) and ((max @minmax) >= $minRunmm); + my $times = 0; + + + my ($ymm, $yET) = mmFromLogs($data_1mm,$logsPath,$ETPath); + + print "ymm = " . join(',',@$ymm) . "\n" if ($debug); + print "yET = " . join(',',@$yET) . "\n" if ($debug); + + my @tET = [0] x scalar (@ET); + + for (my $x = 0; $x < scalar (@ET); $x++) { + print "[ET[$x] = $ET[$x] yET[$x] = @$yET[$x]]\n" if ($debug); + $ET[$x] -= @$yET[$x]; + } + + + my @runTime = (); + for (my $x = 0; $x < scalar (@{$data_1mm->{mmTime}}); $x++) { + my $aET = safe_float($ET[$data_1mm->{crop}[$x]] - $todayRain - @$ymm[$x] - $tadjust); # tadjust is global ? + my $pretimes = $times; +#HP TODO This will determine if a 2nd, 3rd or 4th time is required. + $times = 1 if (($aET / $minRunmm) > 1); #if the minium threshold is met, then run at least once. + $times = int(max ( min( $aET / $maxRunmm, 4), $times)); # int(.999999) = 0 + print "[calc_eto] DB: times=$times aET=$aET minRunm=$minRunmm maxRunm=$maxRunmm\n"; #if ($debug); + print "E: aET[$x] = $aET (" . $aET / $maxRunmm . ") // mm/Day\n" if ($debug); + print "E: times = $times (max " .max (min ( $aET / $maxRunmm,4),$times) . "/min " . min($aET/$maxRunmm,4) ." max(min(". $aET / $maxRunmm . ", 4), $pretimes))\n" if ($debug); + # + # + # @FIXME: this is way too hard to read + +# runTime.append(min(max(safe_int(data['mmTime'][x] * ((aET if aET >= minRunmm else 0)) * (not noWater)), 0), \ +# safe_int(data['mmTime'][x]) * maxRunmm)) + my $tminrun = safe_int($data_1mm->{mmTime}[$x]); + $tminrun = 0 unless $aET >= $minRunmm; + $tminrun = int($tminrun * $aET); + $tminrun = 0 if $noWater; + my $tmaxrun = safe_int($data_1mm->{mmTime}[$x]) * $maxRunmm; + print "E: HP mmTime = " . $data_1mm->{mmTime}[$x] . " tminrun=$tminrun tmaxrum=$tmaxrun\n" if ($debug); + push (@runTime, min ($tminrun, $tmaxrun)); + } + + # ######################################### + # # Real logs will be written already ## + # ######################################### + + print "runTime count [" . scalar (@runTime) . "]\n" if ($debug); + if (open (FILE, ">$logsPath/$fname")) { + my $logData = "["; + for (my $x = 0; $x < scalar (@runTime); $x++) { + for (my $y = 0; $y < $times; $y++) { + my $delim = ""; + $delim = ", " unless (($x == 0) and ($y == 0)); + $logData .= $delim . "[$pid, $x, " . $runTime[$x] . "]"; + } + } + $logData .= "]"; + print FILE $logData; + close (FILE); + } else { + print_log "[calc_eto] ERROR Can't open log file $logsPath/$fname!"; + close (FILE); + } + + + # ######################################### + # # Write final daily water balance ## + # ######################################### + if (open (FILE, ">$ETPath/$fname")) { + my $Data = "["; + for (my $x = 0; $x < scalar (@ET); $x++) { + my $delim = ""; + $delim = ", " unless $x == 0; + $Data .= $delim . $ET[$x]; + } + $Data .= "]"; + print FILE $Data; + close (FILE); + } else { + print_log "[calc_eto] ERROR Can't open ET file $ETPath/$fname!"; + close (FILE); + } + + + # ########################################## + +#HP - ok, this is explained by the opensprinker setup, a program can have up to 4 runtimes +#HP - useful to avoid grass saturation. So if really dry and needs lots of moisture, then run +#HP - multiple programs +#HP - also set the number of times to water to 0 if $noWater is set. + $times = 0 if ($noWater); + + my @startTime = (-1) x 4; + my @availTimes = ( $sun->{rise} - sum( @runTime) / 60, + $sun->{rise} + 60, + $sun->{set} - sum( @runTime) / 60, + $sun->{set} + 60 ); + + + print "[times=$times, sun->{rise}=" . $sun->{rise} . " sum=" . sum(@runTime) / 60 . "]\n";# if ($debug); + + for (my $i = 0; $i < $times; $i++) { + $startTime[$i] = int($availTimes[$i]); + } + my $runTime_str = "[[" . join(',',@startTime) . "],[" . join(',', @runTime) . "]]"; + $msg = "[calc_eto] Current logged ET [" . join(',', @ET) . "]"; + print_log $msg; + $msg_string .= $msg . "\n"; + $msg = "[calc_eto] Current 1mm times [" . join (',', @{$data_1mm->{mmTime}}) . "]"; + print_log $msg; + $msg_string .= $msg . "\n"; + $msg = "[calc_eto] Current Calc time $runTime_str"; + print_log $msg; + $msg_string .= $msg . "\n"; + + if (open (FILE, ">$WPPath/run")) {; + print FILE $runTime_str; + close (FILE); + } else { + print_log "[calc_eto] ERROR Can't open run file $WPPath/run!"; + close (FILE); + } + return $runTime_str; +} +# -[ Data ]--------------------------------------------------------------------- + +sub writewuData { + my ($wuData, $noWater, $wuDataPath) = @_; + my $fname = int((time - (time % 86400) -1) / 86400); + if (open (FILE, ">$wuDataPath/$fname")) { + print FILE "observation_epoch, " . $wuData->{current_observation}->{observation_epoch} . "\n"; + print FILE "weather, " . $wuData->{current_observation}->{weather} . "\n"; + print FILE "temp_c, " . $wuData->{current_observation}->{temp_c} . "\n"; + print FILE "temp_f, " . $wuData->{current_observation}->{temp_f} . "\n"; + print FILE "relative_humidity, " . $wuData->{current_observation}->{relative_humidity} . "\n"; + print FILE "wind_degrees, " . $wuData->{current_observation}->{wind_degrees} . "\n"; + print FILE "wind_mph, " . $wuData->{current_observation}->{wind_mph} . "\n"; + print FILE "wind_kph, " . $wuData->{current_observation}->{wind_kph} . "\n"; + print FILE "precip_1hr_in, " . $wuData->{current_observation}->{precip_1hr_in} . "\n"; + print FILE "precip_1hr_metric, " . $wuData->{current_observation}->{precip_1hr_metric} . "\n"; + print FILE "precip_today_string, " . $wuData->{current_observation}->{precip_today_string} . "\n"; + print FILE "precip_today_in, " . $wuData->{current_observation}->{precip_today_in} . "\n"; + print FILE "precip_today_metric, " . $wuData->{current_observation}->{precip_today_metric} . "\n"; + print FILE "noWater, " . $noWater . "\n"; + close (FILE); + } else { + print_log "[calc_eto] WARNING Can't open wuData file $wuDataPath/$fname for writing!\n"; + } +} + + +#calc_eto_runtimes(".","wu",$config_parms{eto_location},$config_parms{wu_key}); +sub calc_eto_runtimes { + my ($datadir,$method,$loc,$arg1,$arg2) = @_; + my ($rt,$data); + my $success = 0; + if (lc $method eq "file") { + if (open (my $fh, "$arg1")) { + local $/; #otherwise raw_data is empty? + my $raw_data = <$fh>; +# eval { $data = JSON::XS->new->decode($raw_data) }; + eval { $data = decode_json($raw_data) }; + if ($@) { + print_log "[calc_eto] ERROR Problem parsing data $arg1! $@\n"; + } else { + $success = 1; + } + close ($fh); + } else { + print_log "[calc_eto] ERROR Problem opening data $arg1\n"; + close ($fh); + } + } elsif (lc $method eq "wu") { + $data = getwuData($loc,$arg1); + $success = 1 if ($data); + } + if ($success) { + $rt = main_calc_eto($datadir,$loc,$data); + } else { + print_log "[calc_eto] ERROR Data not available.\n"; + } + return $rt; +} + + +sub main_calc_eto { + my ($datadir,$loc,$wuData) = @_; +# -[ Init ]--------------------------------------------------------------------- + $msg_string = ""; + my $msg; + $datadir .= '/' unless ((substr($datadir,-1) eq "/") or (substr($datadir,-1) eq "\\")); + my $logsPath = $datadir . 'logs'; + my $ETPath = $datadir . 'ET'; + my $wuDataPath = "$config_parms{data_dir}/wuData"; + my $WPPath = $datadir . 'weatherprograms'; + + my $tzone; + + my $rainfallsatpoint = 25; + $rainfallsatpoint = $config_parms{eto_rainfallsatpoint} if (defined $config_parms{eto_rainfallsatpoint}); + +######################################### +## We need your latitude and longitude ## +## Let's try to get it with no api call## +######################################### + +# Hey we were given what we needed let's work with it + our ($lat,$t1,$lon) = $loc =~ /^([-+]?\d{1,2}([.]\d+)?),\s*([-+]?\d{1,3}([.]\d+)?)$/; + $lat = "None" unless ($lat); + $lon = "None" unless ($lon); + +# We got a 5+4 zip code, we only need the 5 + $loc =~ s/\-\d\d\d\d//; +# + +# We got a pws id, we don't need to tell wunderground, +# they know how to deal with the id numbers + $loc =~ s/'pws:'//; +# + +# Okay we finally have our loc ready to look up + my $noData = 0; + my ($whttyp,$ploc); + if ($lat eq "None" and $lon eq "None") { + ($whttyp, $ploc, $noData, $tzone, $lat, $lon) = findwuLocation($loc); + } +# Okay if all went well we got what we needed and snuck in a few more items we'll store those somewhere + + if ($lat and $lon) { + if ($lat and $lon and $whttyp and $ploc) { + print_log "[calc_eto] INFO For the $whttyp named: $ploc the lat, lon is: $lat, $lon, and the timezone is $tzone"; + } else { + print_log "[calc_eto] INFO Resolved your lat:$lat, lon:$lon"; + } + $loc = $lat . ',' . $lon; + } else { + if ($noData) { + print_log "[calc_eto] ERROR couldn't reach Weather Underground check connection"; + } else { + print_log "[calc_eto] ERROR $loc can't resolved try another location"; + } + } + +# -[ Main ]--------------------------------------------------------------------- + + our ($offsets) = getwuDataTZOffset($wuData,$tzone); + + unless ($wuData) { + print_log "[calc_eto] ERROR WU data appears to be empty, exiting"; + return "[[-1,-1,-1,-1],[0]]"; + } + +# Calculate an adjustment based on predicted rainfall + my $tadjust = getForecastData($wuData->{forecast}->{simpleforecast}->{forecastday}); + my $sun = getAstronomyData($wuData->{sun_phase}); + my ($todayRain, $noWater, $whyNot) = getConditionsData($wuData->{current_observation}, $wuData->{forecast}->{simpleforecast}->{forecastday}[0],$conditions); + +######################## Quick Ref Names For wuData ######################################## + my $hist = $wuData->{history}->{dailysummary}[0]; + +########################### Required Data ################################################## + $lat = safe_float($lat); + my $tmin = safe_float($hist->{mintempm}); + my $tmax = safe_float($hist->{maxtempm}); + my $tmean = ($tmin + $tmax) / 2; + my $alt = safe_float($wuData->{current_observation}->{display_location}->{elevation}); + my $tdew = safe_float($hist->{meandewptm}); + my $doy = Day_of_Year($hist->{date}->{year}, $hist->{date}->{mon}, $hist->{date}->{mday}); + my $sun_hours = sun_block($wuData,$sun->{rise}, $sun->{set}, $conditions); + my $rh_min = safe_float($hist->{minhumidity}); + my $rh_max = safe_float($hist->{maxhumidity}); + my $rh_mean = ($rh_min + $rh_max) / 2; + my $meanwindspeed = safe_float($hist->{meanwindspdm}); + my $rainfall = min(safe_float($hist->{precipm}), safe_float($rainfallsatpoint)); + +############################################################################################ +## Calculations ## +############################################################################################ +# Calc Rn +print "pl1 [lat=$lat,tmin=$tmin,tmax=$tmax,tmean=$tmean,alt=$alt,tdew=$tdew,doy=$doy,shour=$sun_hours,rmin=$rh_min,rmax=$rh_max,$meanwindspeed,$rainfall,$rainfallsatpoint]\n" if ($debug); + my $e_tmin = &eto::delta_sat_vap_pres($tmin); + my $e_tmax = &eto::delta_sat_vap_pres($tmax); + my $sd = &eto::sol_dec($doy); + my $sha = &eto::sunset_hour_angle($lat, $sd); + my $dl_hours = &eto::daylight_hours($sha); + my $irl = &eto::inv_rel_dist_earth_sun($doy); + my $etrad = &eto::et_rad($lat, $sd, $sha, $irl); + my $cs_rad = &eto::clear_sky_rad($alt, $etrad); + my $Ra = ""; + +print "pl2 [e_tmin=$e_tmin e_tmax=$e_tmax sd=$sd sha=$sha dl_hours=$dl_hours irl=$irl etrad=$etrad cs_rad=$cs_rad]\n" if ($debug); + + my $sol_rad = &eto::sol_rad_from_sun_hours($dl_hours, $sun_hours, $etrad); + $sol_rad = &eto::sol_rad_from_t($etrad, $cs_rad, $tmin, $tmax) unless ($sol_rad); + unless ($sol_rad) { + print_log "[calc_eto] INFO Data for Penman-Monteith ETo not available reverting to Hargreaves ETo\n"; + # Calc Ra + $Ra = $etrad; + print_log "[calc_eto] WARNING Not enough data to complete calculations" unless ($Ra); + } + + $msg = "[calc_eto] RESULTS Sun hours today: $sun_hours"; # tomorrow+2 days forecast rain + print_log $msg; + $msg_string .= $msg . "\n"; + + + my $ea = &eto::ea_from_tdew($tdew); + $ea = &eto::ea_from_tmin($tmin) unless ($ea); + $ea = &eto::ea_from_rhmin_rhmax($e_tmin, $e_tmax, $rh_min, $rh_max) unless ($ea); + $ea = &eto::ea_from_rhmax($e_tmin, $rh_max) unless ($ea); + $ea = &eto::ea_from_rhmean($e_tmin, $e_tmax, $rh_mean) unless ($ea); + print_log "[calc_eto] INFO Failed to set actual vapor pressure" unless ($ea); + + + my $ni_sw_rad = &eto::net_in_sol_rad($sol_rad); + my $no_lw_rad = &eto::net_out_lw_rad($tmin, $tmax, $sol_rad, $cs_rad, $ea); + my $Rn = &eto::net_rad($ni_sw_rad, $no_lw_rad); + +# Calc t + + my $t = ($tmin + $tmax) / 2; + +# Calc ws (wind speed) + + my $ws = &eto::wind_speed_2m($meanwindspeed, 10); + +# Calc es + + my $es = &eto::mean_es($tmin, $tmax); + +print "pl3 [sol_rad=$sol_rad ra=$Ra ea=$ea ni_sw_rad=$ni_sw_rad no_lw_rad=$no_lw_rad rn=$Rn t=$t ws=$ws es=$es]\n" if ($debug); + +# ea done in Rn calcs +# Calc delta_es + + my $delta_es = &eto::delta_sat_vap_pres($t); + +# Calc psy + + my $atmospres = &eto::atmos_pres($alt); + my $psy = &eto::psy_const($atmospres); +print "pl4 [delta_es=$delta_es atmospres=$atmospres psy=$psy]\n" if ($debug); +############################## Print Results ################################### + + $msg = "[calc_eto] RESULTS " . round($tadjust,4) . " mm precipitation forecast for next 3 days"; # tomorrow+2 days forecast rain + print_log $msg; + $msg_string .= $msg . "\n"; + $msg = "[calc_eto] RESULTS " . round($todayRain,4) . " mm precipitation fallen and forecast for today"; # rain fallen today + forecast rain for today + print_log $msg; + $msg_string .= $msg . "\n"; + + # Binary watering determination based on 3 criteria: 1)Currently raining 2)Wind>8kph~5mph 3)Temp<4.5C ~ 40F + if ($noWater) { + $msg = "[calc_eto] RESULTS We will not water because: $whyNot"; + print_log $msg; + $msg_string .= $msg . "\n"; + } + + my ($ETdailyG, $ETdailyS); + if (not $Ra) { + $ETdailyG = round(&eto::ETo($Rn, $t, $ws, $es, $ea, $delta_es, $psy, 0) - $rainfall,4); #ETo for most lawn grasses + $ETdailyS = round(&eto::ETo($Rn, $t, $ws, $es, $ea, $delta_es, $psy, 1) - $rainfall,4); #ETo for decorative grasses, most shrubs and flowers + $msg = "[calc_eto] RESULTS P-M ETo"; + print_log $msg; + $msg_string .= $msg . "\n"; + $msg = "[calc_eto] RESULTS $ETdailyG mm lost by grass"; + print_log $msg; + $msg_string .= $msg . "\n"; + $msg = "[calc_eto] RESULTS $ETdailyS mm lost by shrubs"; + print_log $msg; + $msg_string .= $msg . "\n"; + + } else { + $ETdailyG = round(&eto::hargreaves_ETo($tmin, $tmax, $tmean, $Ra) - $rainfall,4); + $ETdailyS = $ETdailyG; + $msg = "[calc_eto] RESULTS H ETo"; + print_log $msg; + $msg_string .= $msg . "\n"; + + $msg = "[calc_eto] RESULTS $ETdailyG mm lost today"; + print_log $msg; + $msg_string .= $msg . "\n"; + + } + + $msg = "[calc_eto] RESULTS sunrise & sunset in minutes from midnight local time: " . $sun->{rise} . ' ' . $sun->{set}; + print_log $msg; + $msg_string .= $msg . "\n"; + + my $stationID = $wuData->{current_observation}->{station_id}; + $msg = '[calc_eto] RESULTS Weather Station ID: ' . $stationID; + print_log $msg; + $msg_string .= $msg . "\n"; + + my $updateTime = $wuData->{current_observation}->{observation_time}; + $msg = '[calc_eto] RESULTS Weather data ' . $updateTime; + print_log $msg; + $msg_string .= $msg . "\n"; + + my ($rtime) = writeResults($ETdailyG, $ETdailyS, $sun, $todayRain, $tadjust, $noWater, $logsPath, $ETPath, $WPPath); +#Write the WU data to a file. This can be used for the MH weather data and save an api call + writewuData($wuData, $noWater, $wuDataPath); + $msg = "[calc_eto] RESULTS Calculated Schedule: $rtime"; + print_log $msg; + $msg_string .= $msg . "\n"; + if (defined $config_parms{eto_email}) { + print_log "[calc_eto] Emailing results"; + net_mail_send(to => $config_parms{eto_email}, subject => "EvapoTranspiration Results for $Time_Now", text => $msg_string); + } + return ($rtime); +} + + + + + diff --git a/lib/eto.pm b/lib/eto.pm new file mode 100644 index 000000000..006a1c64d --- /dev/null +++ b/lib/eto.pm @@ -0,0 +1,1030 @@ +package eto; + +#Original Name: fao_eto.py +#Purpose: Library for calculating reference evapotransporation (ETo) for +# grass using the FAO Penman-Monteith equation +#Author: Mark Richards +#Copyright: (c) Mark Richards 2010 + + +#Licence +#======= +#This program is free software: you can redistribute it and/or modify +#it under the terms of the GNU General Public License as published by +#the Free Software Foundation, either version 3 of the License, or +#(at your option) any later version. +# +#This program is distributed in the hope that it will be useful, +#but WITHOUT ANY WARRANTY; without even the implied warranty of +#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +#GNU General Public License for more details. +# +#You should have received a copy of the GNU General Public License +#along with this program. If not, see . +# +#Description +#=========== +#A library of functions to allow calculation of reference evapotranspiration +#(ETo) for a grass crop using minimum meteorological data. The methods are based +#on guidelines published by the Food and Argiculture Organisation (FAO) of the +#United Nations in: +# +#Allen, R.G., Pereira, L.S., Raes, D. and Smith, M. (1998) Crop +# evapotranspiration. Guidelines for computing crop water requirements, +# FAO irrigation and drainage paper 56) +# +#Almost all of the functions have been tested against examples given in the FAO +#paper. +# +#Instructions +#============ +#These instructions are a brief summary of those given in Allen et al (1998). +#The data required to calculate the daily, ten-day or monthly evapotranspiration +#over grass using the FAO Penman-Monteith equation are specified below. If +#measured data are not available, many of the variables can be estimated using +#functions in this module. +# +#If insufficient data are available, the alternative, +#data light Hargreaves ETo equation can be used. However, in general, +#estimating solar radiation, vapor pressure and wind speed using the functions +#describedbelow and then calculating evapotranspiration using the Penman-Monteith +#method will provide somewhat more accurate estimates compared to the Hargreaves +#equation. This is dueto the ability of the estimation equations to incorporate +#general climatic characteristics such as high or low wind speed or high or low +#relative humidity into the ETo estimate made using Penman-Monteith. +# +#The Hargreaves equation has a tendency to underpredict under high wind +#conditions(u2 > 3m/s) and to overpredict under conditions of high relative +#humidity. +# +#Monthly (or ten-day) time step +#------------------------------ +#The value of reference evapotranspiration calculated with mean monthly weather +#data is very similar to the average of the daily ETo values calculated with +#average weather data for that month. The follwoing data are required (if using +#a ten-day period substitude the words 'ten-day' in place of 'monthly'): +# +#- monthly average daily maximum and minimum temperature +#- monthly avereage of actual vapour pressure derived from psychrometric, +# dewpoint or relative humidty data. +#- monthly average of daily wind speed data measured at 2 m height (can be +# estimated from measurements made at different heights) +#- monthly avereage of daily net radiation computed from monthly measured short- +# wave radiation or from actual duration of daily sunshine hours. The +# extraterrestrial radiation and daylight hours for a specific day of the +# month can be computed using functions in this module. +#- soil heat flux for monthly periods can be significant when soil is warming in +# spring or cooling in autumn so its value should be determined from the +# mean monthly air tmperatures of the previous and next month (see +# monthly_soil_heat_flux(). + +#Daily time step +#--------------- +#The required meteorological data are: +# +#- minimum and maximum daily air temperature +#- mean daily actual vapour pressure derived from psychrometric, dewpoint +# temperature or relative humidty data (or even just minimum temperature) +#- daily average wind speed measured at 2 m height (can be estimated from +# measurements made at different heights) +#- net radiation measured or computed from solar (shortwave) and longwave +# radiation or from the actual duration of sunshine. The extraterrestrial +# radiation for a specific day of the month should be computed using +# the et_rad() and daylight_hours() functions. +#- as the magnitude of daily soil heat flux beneath a reference grass crop +# is relatively small it may ignored (soil heat flux = 0) for daily time +# steps though if you wish you can calculate it using the +# daily soil_heat_flux() function. +# +#To calculate ETo using the penman_monteith_ETo() function gather the data +#necessary for the function's arguments. It is best to provide measured +#values for the inputs where possible but if some of the data is not +#available then use one of the other functions to estimate the input. +# +#For some input variables there is an order of preference for which function +#to use to estimate the values due to variation +#in the robustness/generality of the different methods. +# +#e.g. If you wish to calculate daily net radiation +#you can estimate it from measured sunshine hours (intermediate option) or +#from the minimum temperature (worst option). +# +#Below is a list of variables for which multiple functions exist along with the +#order of preference for their use: +# +#Actual vapour pressure +#---------------------- +#If measured values are not available then use the following functions +#to estimate AVP (in order of preference): +#1. If dewpoint temperature data are available use ea_from_tdew() +#2. If dry and wet bulb temperatures are available from a psychrometer +# use ea_from_twet_tdry() +#3. If reliable min and max relatiuve humidity data available use +# aea_from_rhmin_rh_max() +#4. If measurement errors of RH are large then use only RH max using +# ea_from_rhmax() +#5. If RH min and RH max are not available but RH mean is then use +# ea_from_rhmean() (but this is less reliable than options 3 or 4) +#6. If no data for the above are available then use ea_from_tmin(). +# This function is less reliable in arid areas where it is recommended that +# 2 deg C is subtracted from Tmin before it is passed to the function +# following Annex 6 of the FAO paper. +# +#Soil heat flux +#-------------- +#For a daily time step soil heat flux is small compared to net radiation +#when the soil is covered by vegetation so it can be assumed to be zero. +#However, it daily soil heat flux can be estimated using daily_soil_heat_flux(). +# +#For a monthy time step soil heat flux is significant and should be estimated +#using: +#1. monthly_soil_heat_flux if temperature data for the previous and next month +# is available or +#2. monthly_soil_heat_flux2 if temeprature for the next month is not available. +# +#Solar (shortwave) radiation +#--------------------------- +#The amount of incoming solar radiation (AKA shortwave radiation) reaching a +#horizontal plane after scattering by the atmosphere. +#If measured values of gross solar radiation are not available the following 2 +#methods are available (in order of preference) to estimate it: +#1. If sunshine duration data are available use sol_rad_from_sun_hours() +#2. Otherwise use sol_rad_from_t() which requires T min and T max data. +# Suitable for coastal or inland areas but not islands. +#3. For island locations (island <= 20 km wide) where no measured values +# are available from elsewhere on the island and the altitude is 0-100m use +# sol_rad_island(). Only suitable for monthly calculations. +# +#Net solar (shortwave) radiation +#------------------------------- +#The amount of solar radiation (sometimes referred to as shortwave radiation) +#that is not reflected by the surface. The methods listed below assume an +#albedo of 0.23 for a grass reference crop. +#Use function net_rad() to estimate net solar radiation for a grass crop. +# +#References +#---------- +#Allen, R.G., Pereira, L.S., Raes, D. and Smith, M. (1998) Crop +# evapotranspiration. Guidelines for computing crop water requirements. +# FAO irrigation and drainage paper 56,FAO, Rome. +#Hargreaves, G.H. and Z.A. Samani (1982) Estimating potential +# evapotranspiration. J. Irrig. and Drain Engr., ASCE, 108(IR3):223-230. +#Hargreaves, G.H. and Z.A. Samani (1985) Reference crop evapotranspiration from +# temperature. Transaction of ASAE 1(2):96-99. +# +#Version history +#--------------- +#1.0 based on python 1.2.01 (29/11/10) - Fixed minor error when converting +# deg C to Kelvin (was adding 273.16 instead of 273.15.) + + +# Global constants +my $PI = 3.14159265; + +sub atmos_pres { + +# Calculates atmospheric pressure (kPa) using equation (7) in +# the FAO paper, page 62. Calculated using a simplification +# of the ideal gas law, assuming 20 deg C for a standard atmosphere. +# +# Arguments: +# alt - elevation/altitude above sea level (m) + + my ($alt) = @_; + + # Raise exceptions + if ($alt < -20 or $alt > 11000) { + print "[eto.pl]: alt=$alt is not in range -20 to 11000 m"; + } else { + my $tmp1 = (293.0 - (0.0065 * $alt)) / 293.0; + my $tmp2 = $tmp1 ** 5.26; + my $atmos_pres = 101.3 * $tmp2; + return $atmos_pres; + } +} + +sub clear_sky_rad { + +# Calculates clear sky radiation [MJ m-2 day-1] based on FAO equation 37 +# which is recommended when calibrated Angstrom values are not available. +# +# Arguments: +# alt - elevation above sea level [m] +# et_rad - extraterrestrial radiation [MJ m-2 day-1] + + my ($alt, $et_rad) = @_; + + # Raise exceptions + if ($alt < -20 or $alt > 8850) { + print "[eto.pl]: altitude=$alt is not in range -20 to 8850 m"; + } elsif ($et_rad < 0.0 or $et_rad > 50.0) { + # TODO: put in a more realistic upper bound for et_rad than this + print "[eto.pl]: et_rad=$et_rad is not in range 0-50"; + } else { + my $clear_sky_rad = (0.00002 * $alt + 0.75) * $et_rad; + return $clear_sky_rad; + } +} + +sub daily_mean_t { + +# Calculates mean daily temperature [deg C] from the daily minimum and +# maximum temperatures. +# +# Arguments: +# tmin - minimum daily temperature [deg C] +# tmax - maximum daily temperature [deg C] + + my ($tmin, $tmax) = @_; + + # Raise exceptions + if ($tmin < -95.0 or $tmin > 60.0) { + print "[eto.pl]: tmin=$tmin is not in range -95 to +60"; + } elsif ($tmax < -95.0 or $tmax > 60.0) { + print "[eto.pl]: tmax=$tmax is not in range -95 to +60"; + + } else { + my $tmean = ($tmax + $tmin) / 2.0; + return $tmean; + } +} + +sub daily_soil_heat_flux { + +# Estimates the daily soil heat flux (Gday) [MJ m-2 day-1] +# assuming a grass crop from the curent air temperature +# and the previous air temperature. The length over time over which the +# current and previous air temperatures are measured are specified by t_len +# which should be greater than 1 day. The calculations are based on FAO +# equation 41. The soil heat capacity is related to its mineral composition +# and water content. The effective soil depth (z) is only 0.10-0.20 m for one +# day. The resluting heat flux can be converted to +# equivalent evaporation [mm day-1] using the equiv_evap() function. +# +# Arguments: +# t_cur - air temperature at tim i (current) [deg C] +# t_prev - air temperature at time i-1 [deg C] +# delta_t - length of time interval between t_cur and t_prev [day] +# soil_heat_cap - soil heat capacity [MJ m-3 degC-1] (default value is 2.1) +# delta_z - effective soil depth [m] (default - 0.1 m following FAO +# recommendation for daily calculations + + my ($t_cur, $t_prev, $delta_t, $soil_heat_cap, $delta_z) = @_; + $soil_heat_cap = 2.1 unless ($soil_heat); + $delta_z = 0.10 unless ($delta_z); + + # Raise exceptions + if ($t_prev < -95.0 or $t_prev > 60.0) { + print "[eto.pl]: t_prev=$t_prev is not in range -95 to +60"; + } elsif ($t_cur < -95.0 or $t_cur > 60.0) { + print "[eto.pl]: t_cur=$t_cur is not in range -95 to +60"; + # for dilay calc delta_t should be greater than 1 day + } elsif ($delta_t < 1.0) { + print "[eto.pl]: delta_t=$delta_t is less than 1 day"; + } elsif ($soil_heat_cap < 1.0 or $soil_heat_cap > 4.5) { + print "[eto.pl]: soil_heat_cap=$soil_heat_cap is not in range 1-4.5"; + } elsif ($delta_z < 0.0 or $delta_z > 200.0) { + print "[eto.pl]: delta_z=$delta_z is not in range 0-200 m"; + } else { + # Assume an effective soil depth of 0.10 m for a daily calculation as per + # FAO recommendation + my $soil_heat_flux = $soil_heat_cap * (($t_cur - $t_prev) / $delta_t) * $delta_z; + return soil_heat_flux + } +} + +sub daylight_hours { + +# Calculates the number of daylight hours from sunset hour angle +# based on FAO equation 34. +# +# Arguments: +# sha - sunset hour angle [rad] + + my ($sha) = @_; + # Raise exceptions + # TODO: Put in check for sunset hour angle + my $daylight_hours = (24.0 / $PI) * $sha; + return $daylight_hours; +} + +sub delta_sat_vap_pres { + +# Calculates the slope of the saturation vapour pressure curve at a given +# temperature (t) [kPa degC-1] based on equation 13 from the FAO paper. For +# use in the Penman-Monteith equation the slope should be calculated using +# mean air temperature. +# +# Arguments: +# t - air temperature (deg C) (use mean air temp for use in Penman-Monteith) + + my ($t) = @_; + + # Raise exceptions + if ($t < -95.0 or $t > 60.0) { + print "[eto.pl]: t=$t is not in range -95 to +60"; + } else { + my $tmp1 = (17.27 * $t) / ($t + 237.3); + my $tmp2 = 4098 * (0.6108 * exp($tmp1)); + my $delta_es = $tmp2 / ($t + 237.3) ** 2; + return $delta_es; + } +} + +sub ea_from_tmin { + +# Calculates actual vapour pressure, ea [kPa] using equation (48) in +# the FAO paper. This method is to be used where humidity data are +# lacking or are of questionable quality. The method assumes that the +# dewpoint temperature is approximately equal to the minimum temperature +# (T_min), i.e. the air is saturated with water vapour at T_min. +# NOTE: This assumption may not hold in arid/semi-arid areas. +# In these areas it may be better to substract 2 deg C from t_min (see +# Annex 6 in FAO paper). +# +# Arguments: +# tmin - daily minimum temperature [deg C] + + my ($tmin) = @_; + + # Raise exception: + if ($tmin < -95.0 or $tmin > 60.0) { + print "[eto.pl]: tmin=$tmin is not in range -95 to 60 deg C"; + } else { + my $ea = 0.611 * exp((17.27 * $tmin)/($tmin + 237.3)); + return $ea; + } +} + +sub ea_from_rhmin_rhmax { + +# Calculates actual vapour pressure [kPa] from relative humidity data +# using FAO equation (17). +# +# Arguments: +# e_tmin - saturation vapour pressure at daily minimum temperature [kPa] +# e_tmax - saturation vapour pressure at daily maximum temperature [kPa] +# rh_min - minimum relative humidity [%] +# rh_max - maximum relative humidity [%] + + my ($e_tmin, $e_tmax, $rh_min, $rh_max) = @_; + + # Raise exceptions: + if ($rh_min < 0 or $rh_min > 100) { + print "[eto.pl]: RH_min=$rh_min is not in range 0-100"; + } elsif ($rh_max < 0 or $rh_max > 100) { + print "[eto.pl]: RH_max=$rh_max is not in range 0-100"; + } else { + my $tmp1 = $e_tmin * ($rh_max / 100.0); + my $tmp2 = $e_tmax * ($rh_min / 100.0); + my $ea = ($tmp1 + $tmp2) / 2.0; + return $ea; + } +} + +sub ea_from_rhmax { + +# Calculates actual vapour pressure [kPa] from maximum relative humidity +# using FAO equation (18). +# +# Arguments: +# e_tmin - saturation vapour pressure at daily minimum temperature [kPa] +# rh_max - maximum relative humidity [%] + + my ($e_tmin, $rh_max) = @_; + + # Raise exceptions: + if ($rh_max < 0 or $rh_max > 100) { + print "[eto.pl]: RH_max=$rh_max is not in range 0-100"; + } else { + return $e_tmin * ($rh_max / 100.0); + } +} + +sub ea_from_rhmean { + +# Calculates actual vapour pressure, ea [kPa] from mean relative humidity +# (the average of RH min and RH max) using FAO equation (19). +# +# Arguments: +# e_tmin - saturation vapour pressure at daily minimum temperature [kPa] +# e_tmax - saturation vapour pressure at daily maximum temperature [kPa] +# rh_mean - mean relative humidity [%] (average between RH min and RH max) + + my ($e_tmin, $e_tmax, $rh_mean) = @_; + # Raise exceptions: + if ($rh_mean < 0 or $rh_mean > 100) { + print "[eto.pl]: RH_mean=$rh_mean is not in range 0-100"; + } else { + my $ea = ($rh_mean / 100.0) * (($e_tmax + $e_tmin) / 2.0); + return $ea; + } +} + +sub ea_from_tdew { + +# Calculates actual vapour pressure, ea [kPa] from the dewpoint temperature +# using equation (14) in the FAO paper. As the dewpoint temperature is the +# temperature to which air needs to be cooled to make it saturated, the +# actual vapour pressure is the saturation vapour pressure at the dewpoint +# temperature. This method is preferable to calculating vapour pressure from +# minimum temperature. +# +# Arguments: +# tdew - dewpoint temperature [deg C] + + my ($tdew) = @_; + # Raise exception: + if ($tdew < -95.0 or $tdew > 65.0) { + # Are these reasonable bounds? + print "[eto.pl]: tdew=$tdew is not in range -95 to +60 deg C"; + } else { + $tmp = (17.27 * $tdew) / ($tdew + 237.3); + $ea = 0.6108 * exp($tmp); + return $ea; + } +} + +sub ea_from_twet_tdry { + +# Calculates actual vapour pressure, ea [kPa] from the wet and dry bulb +# temperatures using equation (15) in the FAO paper. As the dewpoint temp +# is the temp to which air needs to be cooled to make it saturated, the +# actual vapour pressure is the saturation vapour pressure at the dewpoint +# temperature. This method is preferable to calculating vapour pressure from +# minimum temperature. Values for the psychrometric constant of the +# psychrometer (psy_const) can be calculated using the function +# psyc_const_of_psychrometer(). +# +# Arguments: +# twet - wet bulb temperature [deg C] +# tdry - dry bulb temperature [deg C] +# e_twet - saturated vapour pressure at the wet bulb temperature [kPa] +# psy_const - psychrometric constant of the pyschrometer [kPa deg C-1] + + my ($twet, $tdry, $e_twet, $psy_const) = @_; + # Raise exceptions: + if ($twet < -95.0 or$ twet > 65.0) { + # Are these reasonable bounds? + print "[eto.pl]: T_wet=$twet is not in range -95 to +65 deg C"; + } elsif ($tdry < -95.0 or $tdry > 65.0) { + # Are these reasonable bounds? + print "[eto.pl]: T_dry=$tdry is not in range -95 to +65 deg C"; + } else { + my $ea = $e_twet - ($psy_const * ($tdry - $twet)); + return $ea; + } +} + +sub et_rad { + +# Calculates daily extraterrestrial radiation ('top of the atmosphere +# radiation') [MJ m-2 day-1] using FAO equation 21. If you require a monthly +# mean radiation figure then make sure the solar declination, sunset +# hour angle and inverse relative distance between earth and sun +# provided as function arguments have been calculated using +# the day of the year (doy) that corresponds to the middle of the month. +# +# Arguments: +# lat - latitude [decimal degrees] +# sd - solar declination [rad] +# sha - sunset hour angle [rad] +# irl - inverse relative distance earth-sun [dimensionless] + + my ($lat, $sd, $sha, $irl) = @_; + # Raise exceptions + # TODO: raise exceptions for sd and sha + if ($lat < -90.0 or $lat > 90.0) { + print "[eto.pl]: latitude=$lat is not in range -90 to +90"; + } elsif ($irl < 0.9669 or $irl > 1.0331) { + print "[eto.pl]: irl=$irl is not in range 0.9669-1.0331"; + } else { + my $solar_const = 0.0820; # Solar constant [MJ m-2 min-1] + my $lat_rad = $lat * ($PI / 180); # Convert decimal degrees to radians + + # Calculate daily extraterrestrial radiation based on FAO equation 21 + my $tmp1 = (24 * 60) / $PI; + my $tmp2 = $sha * sin($lat_rad) * sin($sd); + my $tmp3 = cos($lat_rad) * cos($sd) * sin($sha); + my $et_rad = $tmp1 * $solar_const * $irl * ($tmp2 + $tmp3); + return $et_rad; + } +} + +sub hargreaves_ETo { + +# Calcultaes evapotranspiration over grass [mm day-1] using the Hargreaves +# ETo equation. Generally, when solar radiation data, relative humidity data +# and/or wind speed data are missing, they should be estimated using the +# procedures/functions outlined in the comments at the top of this file and +# then ETo calculated using the Penman-Monteith equation. +# As an alternative, ETo can be estimated using the Hargreaves ETo equation. +# +# tmin - minimum daily temperaure [deg C] +# tmax - maximum daily temperaure [deg C] +# tmean - mean daily temperaure [deg C] +# Ra - extraterrestrial radiation as equivalent evaporation [mm day-1] + + my ($tmin, $tmax, $tmean, $Ra) = @_; + my $ETo = 0.0023 * ($tmean + 17.8) * ($tmax - $tmin)**0.5 * $Ra; + return $ETo; +} + +sub inv_rel_dist_earth_sun { + +# Calculates the inverse relative distance between earth and sun from +# day of the year using FAO equation 23. +# +# Arguments: +# doy - day of year [between 1 and 366] + + my ($doy) = @_; + # Raise exception + if ($doy < 1 or $doy > 366) { + print "[eto.pl]: doy=$doy is not in range 1-366"; + } else { + my $inv_rel_dist = 1 + (0.033 * cos((2 * $PI / 365)* $doy)); + return $inv_rel_dist; + } +} + +sub mean_es { + +# Calculates mean saturation vapour pressure, es [kPa] using equations (11) +# and (12) in the FAO paper (see references). Mean saturation vapour +# pressure is calculated as the mean of the saturation vapour pressure at +# tmax (maximum temperature) and tmin (minimum temperature). +# +# Arguments: +# tmin - minimum temperature (deg C) +# tmax - maximum temperature (deg C) + + my ($tmin, $tmax) = @_; + # Raise exceptions + if ($tmin < -95.0 or $tmin > 60.0) { + print "[eto.pl]: tmin=$tmin is not in range -95 to +60"; + } elsif ($tmax < -95.0 or $tmax > 60.0) { + print "[eto.pl]: tmax=$tmax is not in range -95 to +60"; + } else { + # Saturation vapour pressure at minimum daily temp + my $tmp1 = (17.27 * $tmin) / ($tmin + 237.3); + my $es_tmin = 0.6108 * exp($tmp1); + + # Saturation vapour pressure at maximum daily temp + my $tmp1 = (17.27 * $tmax) / ($tmax + 237.3); + my $es_tmax = 0.6108 * exp($tmp1); + my $mean_es = ($es_tmin + $es_tmax) / 2.0; + return $mean_es; + } +} + +sub monthly_soil_heat_flux { + +# Estimates the monthly soil heat flux (Gmonth) [MJ m-2 day-1] +# assuming a grass crop from the mean +# air temperature of the previous month and the next month based on FAO +# equation (43). If the air temperature of the next month is not known use +# function monthly_soil_heat_flux2(). The resluting heat flux can be +# converted to equivalent evaporation [mm day-1] using the equiv_evap() +# function. +# +# Arguments: +# t_month_prev - mean air temperature of previous month [deg C] +# t_month2_next - mean air temperature of next month [deg C] + + my ($t_month_prev, $t_month_next) = @_; + # Raise exceptions + if ($t_month_prev < -95.0 or $t_month_prev > 60.0) { + print "[eto.pl]: t_month_prev=$t_month_prev is not in range -95 to +60"; + } elsif ($t_month_next < -95.0 or $t_month_next > 60.0) { + print "[eto.pl]: t_month_next=$t_month_next is not in range -95 to +60"; + } else { + my $soil_heat_flux = 0.07 * ($t_month_next - $t_month_prev); + return $soil_heat_flux; + } +} + +sub monthly_soil_heat_flux2 { + +# Estimates the monthly soil heat flux (Gmonth) [MJ m-2 day-1] +# assuming a grass crop from the mean +# air temperature of the previous and current month based on FAO +# equation (44). If the air temperature of the next month is available use +# monthly_soil_heat_flux() function instead. The resluting heat flux can be +# converted to equivalent evaporation [mm day-1] using the equiv_evap() +# function. +# +# Arguments: +# t_month_prev - mean air temperature of previous month [deg C] +# t_month2_cur - mean air temperature of current month [deg C] + + my ($t_month_prev, $t_month_cur) = @_; + # Raise exceptions + if ($t_month_prev < -95.0 or $t_month_prev > 60.0) { + print "[eto.pl]: t_month_prev=$t_month_prev is not in range -95 to +60"; + } elsif ($t_month_cur < -95.0 or $t_month_cur > 60.0) { + print "[eto.pl]: t_month_cur=$t_month_cur is not in range -95 to +60"; + } else { + my $soil_heat_flux = 0.14 * ($t_month_cur - $t_month_prev); + return $soil_heat_flux; + } +} + +sub net_out_lw_rad { + +# Calculates net outgoing longwave radiation [MJ m-2 day-1] based on +# FAO equation 39. This is the net longwave energy (net energy flux) leaving +# the earth's surface. It is proportional to the absolute temperature of +# the surface raised to the fourth power according to the Stefan-Boltzmann +# law. However, water vapour, clouds, carbon dioxide and dust are absorbers +# and emitters of longwave radiation. This function corrects the Stefan- +# Boltzmann law for humidty (using actual vapor pressure) and cloudiness +# (using solar radiation and clear sky radiation). The concentrations of all +# other absorbers are assumed to be constant. The output can be converted +# to equivalent evapouration [mm day-1] using the equiv_evap() function. +# +# Arguments: +# tmin - absolute daily minimum temperature [deg C] +# tmax - absolute daily maximum temperature [deg C] +# sol_rad - solar radiation [MJ m-2 day-1] +# clear_sky_rad - clear sky radiation [MJ m-2 day-1] +# ea - actual vapour pressure [kPa] + + my ($tmin, $tmax, $sol_rad, $clear_sky_rad, $ea) = @_; + # Raise exceptions + # TODO: raise exceptions for radiation and avp + if ($tmin < -95.0 or $tmin > 60.0) { + print "[eto.pl]: tmin=$tmin is not in range -95 to +60"; + } elsif ($tmax < -95.0 or $tmax > 60.0) { + print "[eto.pl]: tmax=$tmax is not in range -95 to +60"; + } else { + + # Convert temps in deg C to Kelvin + my $tmin_abs = $tmin + 273.15; + my $tmax_abs = $tmax + 273.15; + + my $sb_const = 0.000000004903; # Stefan-Boltzmann constant [MJ K-4 m-2 day-1] + my $tmp1 = $sb_const * ((($tmax_abs ** 4) + ($tmin_abs ** 4)) / 2); + my $tmp2 = 0.34 - (0.14 * sqrt($ea)); + my $tmp3 = 1.35 * ($sol_rad / $clear_sky_rad) - 0.35; + my $net_out_lw_rad = $tmp1 * $tmp2 * $tmp3; + return $net_out_lw_rad; + } +} + +sub net_rad { + +# Calculates daily net radiation [MJ m-2 day-1] at the crop surface +# based on FAO equations 40 assuming a grass reference crop. +# Net radiation is the difference between the incoming net shortwave (or +# solar) radiation and the outgoing net longwave radiation. Output can be +# converted to equivalent evaporation [mm day-1] using the equiv_evap() +# function. +# +# Arguments: +# ni_sw_rad - net incoming shortwave radiation [MJ m-2 day-1] +# no_lw_rad - net outgoing longwave radiation [MJ m-2 day-1] + + my ($ni_sw_rad, $no_lw_rad) = @_; + # Raise exceptions + # TODO: raise exceptions for radiation arguments + my $net_rad = $ni_sw_rad - $no_lw_rad; + return $net_rad; +} + +sub net_in_sol_rad { + +# Calculates net incoming solar (also known as shortwave) +# radiation [MJ m-2 day-1] +# based on FAO equation 38 for a grass reference crop. This is the net +# shortwave radiation resulting from the balance between incoming and +# reflected solar radiation. The output can be converted to +# equivalent evaporation [mm day-1] using the equiv_evap() function. +# +# Arguments: +# sol_rad - (gross) incoming solar radiation [MJ m-2 day-1] + + my ($sol_rad) = @_; + # Raise exceptions + # TODO: Put in sensible boundaries for solar radiation + #if (sol_rad < ?? or sol_rad > ??): + # print "[eto.pl]: sol_rad=%g is not in range 0-366' %sol_rad + + my $grass_albedo = 0.23; # albedo coefficient for grass [dimensionless] + my $net_in_sw_rad = (1 - $grass_albedo) * $sol_rad; + return $net_in_sw_rad; +} + +sub ETo { + +# Calculates the evapotransporation (ETo) [mm day-1] from a hypothetical +# grass reference surface using the FAO Penman-Monteith equation (equation 6). +# +# Arguments: +# Rn - net radiation at crop surface [MJ m-2 day-1] +# t - air temperature at 2 m height [deg C] +# ws - wind speed at 2 m height [m s-1]. If not measured at 2m, +# convert using wind_speed_at_2m() +# es - saturation vapour pressure [kPa] +# ea - actual vapour pressure [kPa] +# delta_es - slope of vapour pressure curve [kPa deg C] +# psy - psychrometric constant [kPa deg C] +# crop - default 0 = grass, 1 = shrubs/tall grasses +# shf - soil heat flux (MJ m-2 day-1] (default = 0, fine for daily +# time step) + + my ($Rn, $t, $ws, $es, $ea, $delta_es, $psy, $crop, $shf) = @_; + $shf = 0.0 unless ($shf); + # TODO: raise exceptions for radiation and avp/svp etc. + if ($t < -95.0 or $t > 60.0) { + print "[eto.pl]: t=$t is not in range -95 to +60"; + } elsif ($ws < 0.0 or $ws > 150.0) { + print "[eto.pl]: ws=$ws is not in range 0-150"; + } else { + my @grass = (900, 0.34); + my @shrub = (1200, 0.38); + my $Cn = $grass[0]; + my $Cd = $grass[1]; + if ($crop == 1) { + $Cn = $shrub[0]; + $Cd = $shrub[1]; + } + # Convert t in deg C to deg Kelvin + $t += 273.15; + # Calculate evapotranspiration (ET0) + my $a1 = 0.408 * ($Rn - $shf) * $delta_es / ($delta_es + ($psy * (1 + $Cd * $ws))); + my $a2 = $Cn * $ws / $t * ($es - $ea) * $psy / ($delta_es + ($psy * (1 + $Cd * $ws))); + my $ETo = $a1 + $a2; + return $ETo; + } +} + +sub psy_const { + +# Calculates the psychrometric constant (kPa degC-1) using equation (8) +# in the FAO paper (see references below) page 95. This method assumes that +# the air is saturated with water vapour at T_min. This assumption may not +# hold in arid areas. +# +# Arguments: +# atmos_pres - atmospheric pressure [kPa] + + my ($atmos_pres) = @_; + # TODO: raise exception if atmos_press outside sensible bounds + return (0.000665 * $atmos_pres); +} + +sub psy_const_of_psychrometer { + +# Calculates the psychrometric constant [kPa deg C-1] for different +# types of psychrometer at a given atmospheric pressure using FAO equation +# 16. +# +# Arguments: +# psychrometer - integer between 1 and 3 which denotes type of psychrometer +# - 1 = ventilated (Asmann or aspirated type) psychrometer with +# an air movement of approx. 5 m s-1 +# - 2 = natural ventilated psychrometer with an air movement +# of approx. 1 m s-1 +# - 3 = non ventilated psychrometer installed indoors +# atmos_pres - atmospheric pressure [kPa] + + my ($psychrometer, $atmos_pres) = @_; + # TODO: raise exception if atmos_press outside sensible bounds + if ($psychrometer < 1 or $psychrometer > 3) { + print "[eto.pl]: psychrometer=$psychrometer not in range 1-3"; + } else { + # Assign values to coefficient depending on type of ventilation of the + # wet bulb + my $psy_coeff; + if ($psychrometer == 1) { + $psy_coeff = 0.000662; + } elsif ($psychrometer == 2) { + $psy_coeff = 0.000800; + } elsif ($psychrometer == 3) { + $psy_coeff = 0.001200; + } + + my $pys_const = $psy_coeff * $atmos_pres; + return $pys_const; + } +} + +sub rad2equiv_evap { + +# Converts radiation in MJ m-2 day-1 to the equivalent evaporation in +# mm day-1 assuming a grass reference crop using FAO equation 20. +# Energy is converted to equivalent evaporation using a conversion +# factor equal to the inverse of the latent heat of vapourisation +# (1 / lambda = 0.408). +# +# Arguments: +# energy - energy e.g. radiation, heat flux [MJ m-2 day-1] + + my ($energy) = @_; + # Determine the equivalent evaporation [mm day-1] + my $equiv_evap = 0.408 * $energy; + return $equiv_evap; +} + +sub rh_from_ea_es { + +# Calculates relative humidity as the ratio of actual vapour pressure +# to saturation vapour pressure at the same temperature (see FAO paper +# p. 67). +# +# ea - actual vapour pressure [units don't matter as long as same as es] +# es - saturated vapour pressure [units don't matter as long as same as ea] + + my ($ea, $es) = @_; + return 100.0 * $ea / $es; +} + +sub sol_dec { + +# Calculates solar declination [rad] from day of the year based on FAO +# equation 24. +# +# Arguments: +# doy - day of year (between 1 and 366) + + my ($doy) = @_; + # Raise exceptions + if ($doy < 1 or $doy > 366) { + print "[eto.pl]: doy=$doy is not in range 1-366"; + } else { + + # Calculate solar declination [radians] using FAO eq. 24 + $solar_dec = 0.409 * sin(((2 * $PI / 365) * $doy - 1.39)); + return $solar_dec; + } +} + +sub sol_rad_from_sun_hours { + +# Calculates incoming solar (or shortwave) radiation [MJ m-2 day-1] +# (radiation hitting a horizontal plane after scattering by the atmosphere) +# from relative sunshine duration based on FAO equations 34 and 35. +# If measured radiation data are not available this +# method is preferable to calculating solar radiation from temperature . +# If a monthly mean is required then divide the monthly number +# of sunshine hours by number of days in month and ensure that et_rad and +# daylight hours was calculated using the day of the year that +# corresponds to the middle of the month. +# +# Arguments: +# dl_hours - number of daylight hours [hours] +# sun_hours - sunshine duration [hours] +# et_rad - extraterrestrial radiation [MJ m-2 day-1] + + my ($dl_hours, $sun_hours, $et_rad) = @_; + # Raise exceptions + # TODO: Raise exception for et_rad + if ($sun_hours < 0 or $sun_hours > 24) { + print "[eto.pl]: sunshine hours=$sun_hours is not in range 0-24"; + } elsif ($dl_hours < 0 or $dl_hours > 24) { + print "[eto.pl]: daylight hours=$dl_hours is not in range 0-24"; + } else { + + # Use default values of regression constants (Angstrom values) + # recommended by FAO when calibrated values are unavailable. + my $a = 0.25; + my $b = 0.50; + my $solar_rad = ($b * $sun_hours / $dl_hours + $a) * $et_rad; + return $solar_rad; + } +} + +sub sol_rad_from_t { + +# Calculates incoming solar (or shortwave) radiation (Rs) [MJ m-2 day-1] +# (radiation hitting a horizontal plane after scattering by the atmosphere) +# from min and max temperatures together with +# an empirical adjustment coefficient for 'interior' and +# 'coastal' regions. The formula is based on FAO equation 50 which +# is the Hargreaves' radiation formula (Hargreaves and Samani, 1982, 1985). +# This method should be used only when solar radiation or sunshine hours data +# are not available. It is only recommended for locations where it is not +# possible to use radiation data from a regional station (either because +# climate conditions are hetergeneous or data are lacking). +# NOTE: this method is not suitable for island locations +# due to the moderating effects of the surrounding water. +# +# Arguments: +# et_rad - extraterrestrial radiation [MJ m-2 day-1] +# cs_rad - clear sky radiation [MJ m-2 day-1] +# tmin - daily minimum temperature [deg C] +# tmax - daily maximum temperature [deg C] +# coastal - True if site is a coastal location, situated on or adjacent to +# coast of a large land mass and where air masses are influence +# by a nearby water body, False if interior location where land +# mass dominates and air masses are not strongly influenced by a +# large water body. -999 indicates no data. + + my ($et_rad, $cs_rad, $tmin, $tmax, $coastal) = @_; + $costal = -999 unless ($coastal); + # Raise exceptions + # TODO: raise exceptions for cs_rad + if ($tmin < -95.0 or $tmin > 60.0) { + print "[eto.pl]: tmin=$tmin is not in range -95 to +60"; + } elsif ($tmax < -95.0 or $tmax > 60.0) { + print "[eto.pl]: tmax=$tmax is not in range -95 to +60"; + } else { + + # determine value of adjustment coefficient [deg C-0.5] for + # coastal/interior locations + my $adj; + if (lc $coastal eq "true") { + $adj = 0.19; + } elsif (lc $coastal eq "false") { + $adj = 0.16; + } else { + # hedge our bets and give a mean adjustment values and issue a warning + $adj = 0.175; + print "[eto.pl]: WARNING! Location not specified as coastal or interior for calculation of solar radiation. Using defalut adjustment factor."; + } + my $solar_rad = $adj * sqrt($tmax - $tmin) * $et_rad; + + # The solar radiation value is constrained (<=) by the clear sky radiation + $solar_rad = $cs_rad if ($solar_rad > $cs_rad); + + return $solar_rad; + } +} + +sub sol_rad_island { + +# Estimates incoming solar (or shortwave) radiation [MJ m-2 day-1] +# (radiation hitting a horizontal plane after scattering by the atmosphere) +# for an island location using FAO equation 51. An island is defined as a +# land mass with width perpendicular to the coastline <= 20 km. Use this +# method only if radiation data from elsewhere on the island is not +# available. NOTE: This method is only applicable for low altitudes (0-100 m) +# and monthly calculations. +# +# Arguments: +# et_rad - extraterrestrial radiation [MJ m-2 day-1] + + my ($et_rad) = @_; + my $solar_rad = (0.7 * $et_rad) - 4.0; + return $solar_rad; +} + +sub sunset_hour_angle { + +# Calculates sunset hour angle [rad] from latitude and solar +# declination using FAO equation 25. +# +# Arguments: +# lat - latitude [decimal degrees] Note: value should be negative if it is +# degrees south, positive if degrees north +# sd - solar declination [rad] + + my ($lat, $sd) = @_; + # TODO: Raise exception for sd + # Raise exceptions + if ($lat < -90.0 or $lat > 90.0) { + print "[eto.pl]: latitude=$lat is not in range -90 - 906"; + } else { + + # Convert latitude from decimal degrees to radians + my $lat_rad = $lat * ($PI / 180); + + # Calculate sunset hour angle (sha) [radians] from latitude and solar + # declination using FAO equation 25 + my $sha = acos(-1 * tan($lat_rad) * tan($sd)); + return $sha; + } +} + +sub wind_speed_2m { + +# Converts wind speeds measured at different heights above the soil +# surface to wind speed at 2 m above the surface, assuming a short grass +# surface. Formula based on FAO equation 47. +# +# Arguments: +# meas_ws - measured wind speed [m s-1] +# z - height of wind measurement above ground surface [m] + + my ($meas_ws, $z) = @_; + # Raise exceptions + if ($meas_ws < 0.0 or $meas_ws > 150.0) { + print "[eto.pl]: meas_ws=$meas_ws is not in range 0-150 m s-1"; + } elsif ($z < 0.0 or $z > 100.0) { + print "[eto.pl]: z=$z is not in range 0-100 m"; + } else { + + my $tmp1 = (67.8 * $z) - 5.42; + my $ws2m = $meas_ws * (4.87 / log($tmp1)); + return $ws2m; + } +} + +sub acos { + return atan2( sqrt(1 - $_[0] * $_[0]), $_[0] ); +} + +sub tan { + return sin($_[0]) / cos($_[0]); +}