User Guide: 2 Astronomy and Atmosphere

Radiation, astronomy and atmosphere

The User Guide is divided into two volumes: 1. Radiation and 2. Astronomy. The second volume also describes functions related to water vapour in the atmosphere.

Getting started

We load two packages, our ‘photobiology’ and ‘lubridate’, as they will be used in the examples.

library(photobiology)
library(lubridate)
# if installed, we use 'lutz' to lookup time zones from geocodes
eval_lutz <- requireNamespace("lutz", quietly = TRUE)
if (eval_lutz) {library(lutz)}

Introduction

Most organisms, including plants and animals, have circadian internal clocks. These clocks are entrained to the day-night cycle through perception of light. For example, night length informs plants about seasons of the year. This information allows the synchronization of developmental events like flowering to take place at the “right” time of the year.

From the point of view of interactions between light and vegetation, the direction of the direct light beam is of interest. Hence, the position of the sun in the sky is also important for photobiology. This is the reason for the inclusion of astronomical calculations about the sun in this package. On the other hand, such calculations are also highly relevant to other fields including solar energy.

The functions and methods described in this volume return either values that represent angles or times. In the current version angles are always expressed in degrees. In the case of times, the unit of expression, can be changed through parameter unit.out which accepts the following arguments "datetime", "hours", "minutes", "seconds". For backwards compatibility "date" is also accepted as equivalent to "datetime" but deprecated.

All astronomical computations rely on the algorithms of Meuss (1998), and consequently returned values are very precise. However, these algorithms are computationally rather costly. Contrary to other faster algorithms, they give reliable estimates even for latitudes near the poles.

However, at high latitudes due to the almost tangential path of the sun to the horizon, atmospheric effects that slightly alter the apparent elevation of the sun have much larger effects on the apparent timing of events given that the solar elevation angle changes at a slower rate than at lower latitudes.

Position of the sun

The position of the sun at an arbitrary geographic locations and time instant can be described with two angles: elevation above the horizon and azimuth angle relative to the geographic North. If expressed in degrees, solar elevation (\(h\)) varies between -90 and 90 degrees, while being visible when angles are positive and otherwise occluded below the horizon. Azimuth angle (\(\alpha\)) varies clockwise from North between 0 and 360 degrees. Zenith angle (\(z\)), provides the same information as the elevation angle but using the zenith as starting point, consequently taking values between 0 and 180 degrees, such that \(z = 90 - h\). Declination angle describes the angle between the plane of the Equator and the plane of the Earth’s orbit around the sun, and varies with the seasons of the year.

The function sun_angles returns location, civil time, local solar time, the azimuth in degrees eastwards, elevation in degrees above the horizon, declination, the equation of time and the hour angle.

In versions up to 0.9.11 in addition to parameter geocode most functions also had the redundant formal parameters lon and lat which were removed in version 0.9.12. This change may require users’ scripts to be revised.

In versions 0.9.16 and later the code has been optimized for performance with time vectors, making massive calculations such as the sun position for every minute of the year quite fast (couple of seconds). We keep, however, examples with rather short vectors.

For calculation of the position of the sun we need to supply geographic coordinates and a time instant. The time instant passed as argument should be a POSIXct vector, possibly of length one. The easiest way create date and time constant values is to use package ‘lubridate’.

The object to be supplied as argument for geocode is a data frame with variables lon and lat giving the location on Earth. This matches the return value of functions tidygeocoder::geo_osm(), tidygeocoder::geo_google() and ggmap::geocode(), functions that can be used to find the coordinates using an address entered as a character string understood by the OSM or Google maps APIs (Google requires an API key and registration, while OSM is open). We use the “geocode” for Helsinki, defined explicitly rather than searched for.

my.geocode <- data.frame(lat = 60.16, lon = 24.93, address = "Helsinki")

Be aware that to obtain correct the time zone must be correctly set for the argument passed to time. To obtain results based on local time, this time zone needs to be set in the POSIXct object or passed as a argument to tz. In the examples we use functions from package ‘lubridate’ for working with times and dates. The argument passed to parameter time can be a “vector” of POSIXct values. The returned value is a data.frame.

The position of the sun at Helsinki, at the given instant in time for time zone Eastern European Time.

sun_angles(time = ymd_hms("2017-06-20 08:00:00", tz = "EET"), geocode = my.geocode)

## # A tibble: 1 × 12
##   time                tz    solartime  longitude latitude address  azimuth
##   <dttm>              <chr> <solar_tm>     <dbl>    <dbl> <chr>      <dbl>
## 1 2017-06-20 08:00:00 EET   06:38:09        24.9     60.2 Helsinki    85.8
## # ℹ 5 more variables: elevation <dbl>, declination <dbl>, eq.of.time <dbl>,
## #   hour.angle <dbl>, distance <dbl>

Functions have defaults for their arguments, but rarely Greenwich will be the location you are interested in. Current UTC time is probably a useful default in some cases.

sun_angles()

## # A tibble: 1 × 12
##   time                tz    solartime  longitude latitude address   azimuth
##   <dttm>              <chr> <solar_tm>     <dbl>    <dbl> <chr>       <dbl>
## 1 2024-08-25 19:11:36 UTC   19:09:43           0     51.5 Greenwich    290.
## # ℹ 5 more variables: elevation <dbl>, declination <dbl>, eq.of.time <dbl>,
## #   hour.angle <dbl>, distance <dbl>

A vector of times is accepted as argument, and as performance is optimized for this case, vectorization will markedly improve performance compared to multiple calls to the function. The vector of times can be created on the fly, or stored beforehand.

sun_angles(time = ymd_hms("2014-01-01 0:0:0", tz = "EET") + hours(c(0, 6, 12)), 
           geocode = my.geocode)

## # A tibble: 3 × 12
##   time                tz    solartime  longitude latitude address  azimuth
##   <dttm>              <chr> <solar_tm>     <dbl>    <dbl> <chr>      <dbl>
## 1 2014-01-01 00:00:00 EET   23:36:26        24.9     60.2 Helsinki   351. 
## 2 2014-01-01 06:00:00 EET   05:36:19        24.9     60.2 Helsinki    97.0
## 3 2014-01-01 12:00:00 EET   11:36:12        24.9     60.2 Helsinki   174. 
## # ℹ 5 more variables: elevation <dbl>, declination <dbl>, eq.of.time <dbl>,
## #   hour.angle <dbl>, distance <dbl>

my.times <- ymd_hms("2014-01-01 0:0:0", tz = "EET") + hours(c(0, 6, 12))
sun_angles(time = my.times, geocode = my.geocode)

## # A tibble: 3 × 12
##   time                tz    solartime  longitude latitude address  azimuth
##   <dttm>              <chr> <solar_tm>     <dbl>    <dbl> <chr>      <dbl>
## 1 2014-01-01 00:00:00 EET   23:36:26        24.9     60.2 Helsinki   351. 
## 2 2014-01-01 06:00:00 EET   05:36:19        24.9     60.2 Helsinki    97.0
## 3 2014-01-01 12:00:00 EET   11:36:12        24.9     60.2 Helsinki   174. 
## # ℹ 5 more variables: elevation <dbl>, declination <dbl>, eq.of.time <dbl>,
## #   hour.angle <dbl>, distance <dbl>

Geocodes for several locations can be stored in a data frame with multiple rows.

two.geocodes <- data.frame(lat = c(60.16, 65.02), 
                           lon = c(24.93, 25.47),
                           address = c("Helsinki", "Oulu"))
sun_angles(time = my.times, geocode = two.geocodes)

## # A tibble: 6 × 12
##   time                tz    solartime  longitude latitude address  azimuth
##   <dttm>              <chr> <solar_tm>     <dbl>    <dbl> <chr>      <dbl>
## 1 2014-01-01 00:00:00 EET   23:36:26        24.9     60.2 Helsinki   351. 
## 2 2014-01-01 06:00:00 EET   05:36:19        24.9     60.2 Helsinki    97.0
## 3 2014-01-01 12:00:00 EET   11:36:12        24.9     60.2 Helsinki   174. 
## 4 2014-01-01 00:00:00 EET   23:38:36        25.5     65.0 Oulu       353. 
## 5 2014-01-01 06:00:00 EET   05:38:29        25.5     65.0 Oulu        95.4
## 6 2014-01-01 12:00:00 EET   11:38:22        25.5     65.0 Oulu       175. 
## # ℹ 5 more variables: elevation <dbl>, declination <dbl>, eq.of.time <dbl>,
## #   hour.angle <dbl>, distance <dbl>

When spectra contain suitable metadata, the position of the sun for the spectral irradiance data measurement can be easily obtained.

sun_angles(time = getWhenMeasured(sun.spct), geocode = getWhereMeasured(sun.spct))

## # A tibble: 1 × 12
##   time                tz    solartime  longitude latitude address        azimuth
##   <dttm>              <chr> <solar_tm>     <dbl>    <dbl> <chr>            <dbl>
## 1 2010-06-22 09:51:00 UTC   11:28:54        25.0     60.2 Kumpula, Hels…    168.
## # ℹ 5 more variables: elevation <dbl>, declination <dbl>, eq.of.time <dbl>,
## #   hour.angle <dbl>, distance <dbl>

If what is needed is only one of the solar angles, functions returning vectors instead of data frames can be useful. (In their current implementation these functions do not have improved performance compared to sun_angles())

sun_elevation(time = my.times, geocode = my.geocode)

## [1] -52.639345 -22.722495   6.710245

sun_zenith_angle(time = my.times, geocode = my.geocode)

## [1] 142.63935 112.72250  83.28976

sun_azimuth(time = my.times, geocode = my.geocode)

## [1] 351.04757  96.98377 174.48767

Times of sunrise, solar noon and sunset

Convenience functions sunrise_time(), sunset_time(), noon_time(), day_length() and night_length() have all the same parameter signature and are wrappers on function day_night(). Function day_night returns a data frame containing all the quantities returned by the other functions.

These functions are vectorized for their date and geocode parameters. They use as default location Greenwich in the U.K., and default time zone “UTC”. The date is given by default by the current date based on “UTC”. Universal Time Coordinate (“UTC”) is the reference used to describe differences among time zones and is never modified by daylight saving time or summer time. The difference between EET (Eastern European Time) and UTC is +2 hours in winter and with EEST (Eastern European Summer Time) +3 hours.

Latitude and longitude are passed through a geocode (a data frame). If the returned value is desired in the local time coordinates of the argument passed to geocode, the time zone should match the geographic coordinates. If geocodes contain a variable "address" it will be copied to the output.

In some of the examples below we reuse the geocode data frames created above, and we here create a vector of datetime objects differing in their date. The default time zone of function ymd() is NULL, so we override it with EET to match the geocodes for Finnish cities.

dates <- ymd("2015-03-01", tz = "EET") + months(0:5)
dates

## [1] "2015-03-01 EET"  "2015-04-01 EEST" "2015-05-01 EEST" "2015-06-01 EEST"
## [5] "2015-07-01 EEST" "2015-08-01 EEST"

As first example we compute the sunrise time for the current day in Helsinki, with the result returned eithe in UTC or EET coordinates. Time-zone names based on continent and city are also supported, and are to be preferred for past dates and the relationship between time zones and locations are affected by changes in country boundaries and in national laws.

sunrise_time(now("UTC"), geocode = my.geocode)

## [1] "2024-08-25 02:59:44 UTC"

sunrise_time(now("EET"), geocode = my.geocode)

## [1] "2024-08-25 05:59:44 EEST"

sunrise_time(now("Europe/Helsinki"), geocode = my.geocode)

## [1] "2024-08-25 05:59:44 EEST"

Be aware of the behaviour of functions ymd(), dmy(), ym() and my() from package ‘lubridate’. A function call like ymd("2015-03-01", tz = "UTC") returns a POSIXct object while a call like ymd("2015-03-01") is equivalent to ymd("2015-03-01", tz = NULL) and returns a Date object. Date objects do not carry time zone information in the way POSIXct objects do, and consequently Dates do not uniquely match a period between two absolute instants in time, but rather between two instants in local time. Given these features, it is safer to use POSIXct objects as arguments to the functions from package ‘photobiology’. Function today() always returns a Date while function now() always returns a POSIXct, independently of the argument passed to their tzone parameter. Consequently it is preferable to use now(), but if you do use today() make sure to path the same value as argument to parameter tzone of today() and to parameter tz of the functions from package ‘photobiology’. An instant in time and time zone strictly define a corresponding date at any location on Earth.

The default for time zone is that of the POSIXct value passed as argument to parameter date, but this can be overridden.

sunrise_time(geocode = my.geocode)

## [1] "2024-08-25 02:59:44 UTC"

sunrise_time(date = now("UTC"), geocode = my.geocode)

## [1] "2024-08-25 02:59:44 UTC"

sunrise_time(date = now("UTC"), tz = "UTC", geocode = my.geocode)

## [1] "2024-08-25 02:59:44 UTC"

sunrise_time(date = now("EET"), geocode = my.geocode)

## [1] "2024-08-25 05:59:44 EEST"

sunrise_time(date = now("EET"), tz = "EET", geocode = my.geocode)

## [1] "2024-08-25 05:59:44 EEST"

sunrise_time(today("EET"), tz = "EET", geocode = my.geocode)

## [1] "2024-08-25 05:59:44 EEST"

We can also compute the time of solar noon and of sunset.

noon_time(now("UTC"), geocode = my.geocode)

## [1] "2024-08-25 10:22:14 UTC"

noon_time(now("EET"), geocode = my.geocode)

## [1] "2024-08-25 13:22:14 EEST"

sunset_time(now("UTC"), geocode = my.geocode)

## [1] "2024-08-25 17:44:45 UTC"

sunset_time(now("EET"), geocode = my.geocode)

## [1] "2024-08-25 20:44:45 EEST"

Functions day_length() and night_length() return the length of time, by default in hours and fraction.

day_length(dates, geocode = my.geocode)

## [1] 10.34596 13.19241 15.90766 18.27962 18.80811 16.97567

night_length(dates, geocode = my.geocode)

## [1] 13.654040 10.807592  8.092343  5.720384  5.191888  7.024327

Southern hemisphere latitudes as well as longitudes to the West of the Greenwich meridian should be supplied as negative numbers.

sunrise_time(dates, geocode = data.frame(lat = 60, lon = 25))

## [1] "2015-02-28 07:21:52 EET"  "2015-03-31 06:48:41 EEST"
## [3] "2015-04-30 05:20:50 EEST" "2015-05-31 04:10:49 EEST"
## [5] "2015-06-30 04:01:10 EEST" "2015-07-31 04:58:15 EEST"

noon_time(dates, geocode = data.frame(lat = 60, lon = 25))

## [1] "2015-02-28 12:32:34 EET"  "2015-03-31 13:24:14 EEST"
## [3] "2015-04-30 13:17:15 EEST" "2015-05-31 13:17:38 EEST"
## [5] "2015-06-30 13:23:36 EEST" "2015-07-31 13:26:24 EEST"

sunrise_time(dates, geocode = data.frame(lat = -60, lon = 25))

## [1] "2015-02-28 05:29:25 EET"  "2015-03-31 07:46:18 EEST"
## [3] "2015-04-30 08:58:10 EEST" "2015-05-31 10:04:25 EEST"
## [5] "2015-06-30 10:24:25 EEST" "2015-07-31 09:37:27 EEST"

noon_time(dates, geocode = data.frame(lat = -60, lon = 25))

## [1] "2015-02-28 12:32:34 EET"  "2015-03-31 13:24:14 EEST"
## [3] "2015-04-30 13:17:15 EEST" "2015-05-31 13:17:38 EEST"
## [5] "2015-06-30 13:23:36 EEST" "2015-07-31 13:26:24 EEST"

The angle used in the twilight calculation can be supplied, either as the name of a standard definition, or as an angle in degrees (negative for sun positions below the horizon). Positive angles can be used when the time of sun occlusion behind a building, mountain, or other obstacle needs to be calculated. The default for twilight is "none" meaning that times correspond to the occlusion of the upper rim of the sun disk below the theoretical horizon.

sunrise_time(ymd("2017-03-21", tz = "EET"), 
             tz = "EET", 
             geocode = my.geocode,
             twilight = "civil")

## [1] "2017-03-20 05:38:58 EET"

sunrise_time(ymd("2017-03-21", tz = "EET"), 
             tz = "EET", 
             geocode = my.geocode,
             twilight = -10)

## [1] "2017-03-20 05:05:45 EET"

sunrise_time(ymd("2017-03-21", tz = "EET"), 
             tz = "EET", 
             geocode = my.geocode,
             twilight = +12)

## [1] "2017-03-20 08:06:15 EET"

Parameter unit.out can be used to obtain the returned value expressed as time-of-day in hours, minutes, or seconds since midnight, instead of the default datetime.

sunrise_time(ymd("2017-03-21", tz = "EET"), 
             tz = "EET", 
             geocode = my.geocode)

## [1] "2017-03-20 06:20:46 EET"

sunrise_time(ymd("2017-03-21", tz = "EET"), 
             tz = "EET", 
             geocode = my.geocode,
             unit.out = "hours")

## [1] 6.346365

Similarly, the units can also be selected for the values returned by day_length() and night_length().

day_length(dates, geocode = my.geocode, unit.out = "days")

## [1] 0.4310817 0.5496837 0.6628190 0.7616507 0.7836713 0.7073197

night_length(dates, geocode = my.geocode, unit.out = "days")

## [1] 0.5689183 0.4503163 0.3371810 0.2383493 0.2163287 0.2926803

The core function is called day_night() and returns a data frame containing both the input values and the results of the calculations. See above for convenience functions useful in the case one needs only one of the calculated variables. In other cases it is more efficient to compute the whole data frame and later select the columns of interest.

day_night(dates[1:3], 
          geocode = my.geocode)

## # A tibble: 3 × 12
##   day                 tz    twilight.rise twilight.set longitude latitude
##   <dttm>              <chr>         <dbl>        <dbl>     <dbl>    <dbl>
## 1 2015-02-28 00:00:00 EET               0            0      24.9     60.2
## 2 2015-03-31 00:00:00 EET               0            0      24.9     60.2
## 3 2015-04-30 00:00:00 EET               0            0      24.9     60.2
## # ℹ 6 more variables: address <chr>, sunrise <dbl>, noon <dbl>, sunset <dbl>,
## #   daylength <dbl>, nightlength <dbl>

The default for unit.out is "hours" with decimal fractions, as seen above. However other useful units like "seconds", "minutes", and "days" can be useful.

day_night(dates[1:3], 
          geocode = my.geocode, 
          unit.out = "days")

## # A tibble: 3 × 12
##   day                 tz    twilight.rise twilight.set longitude latitude
##   <dttm>              <chr>         <dbl>        <dbl>     <dbl>    <dbl>
## 1 2015-02-28 00:00:00 EET               0            0      24.9     60.2
## 2 2015-03-31 00:00:00 EET               0            0      24.9     60.2
## 3 2015-04-30 00:00:00 EET               0            0      24.9     60.2
## # ℹ 6 more variables: address <chr>, sunrise <dbl>, noon <dbl>, sunset <dbl>,
## #   daylength <dbl>, nightlength <dbl>

Finally it is also possible to have the timing of solar events returned as POSIXct time values, in which case lengths of time are still returned as fractional hours.

day_night(dates[1:3], 
          geocode = my.geocode, 
          unit.out = "datetime")

## # A tibble: 3 × 12
##   day                 tz    twilight.rise twilight.set longitude latitude
##   <dttm>              <chr>         <dbl>        <dbl>     <dbl>    <dbl>
## 1 2015-02-28 00:00:00 EET               0            0      24.9     60.2
## 2 2015-03-31 00:00:00 EET               0            0      24.9     60.2
## 3 2015-04-30 00:00:00 EET               0            0      24.9     60.2
## # ℹ 6 more variables: address <chr>, sunrise <dttm>, noon <dttm>,
## #   sunset <dttm>, daylength <dbl>, nightlength <dbl>

When multiple times and locations are supplied as arguments, the values returned are for all combinations of locations and times.

day_night(dates[1:3], 
          geocode = two.geocodes)

## # A tibble: 6 × 12
##   day                 tz    twilight.rise twilight.set longitude latitude
##   <dttm>              <chr>         <dbl>        <dbl>     <dbl>    <dbl>
## 1 2015-02-28 00:00:00 EET               0            0      24.9     60.2
## 2 2015-03-31 00:00:00 EET               0            0      24.9     60.2
## 3 2015-04-30 00:00:00 EET               0            0      24.9     60.2
## 4 2015-02-28 00:00:00 EET               0            0      25.5     65.0
## 5 2015-03-31 00:00:00 EET               0            0      25.5     65.0
## 6 2015-04-30 00:00:00 EET               0            0      25.5     65.0
## # ℹ 6 more variables: address <chr>, sunrise <dbl>, noon <dbl>, sunset <dbl>,
## #   daylength <dbl>, nightlength <dbl>

Solar time

In field research it is in many cases preferable to sample or measure, and present and plot data based on local solar time. A new class is defined in package ‘photobiology’, with a print method, a constructor, a conversion function and a class query function.

The constructor takes as arguments a POSIXct object describing and instant in time and a geocode describing the geographic coordinates.

Paris.geo <- data.frame(lon = 2.352222, lat = 48.85661, address = "Paris")
Paris.time <- ymd_hms("2016-09-30 06:00:00", tz = "UTC")
solar_time(Paris.time, geocode = Paris.geo)

## [1] "06:19:28"

solar_time(Paris.time, geocode = Paris.geo, unit.out = "datetime")

## [1] "2016-09-30 06:19:28 solar"

my.solar.t <- solar_time(Paris.time, geocode = Paris.geo)
is.solar_time(my.solar.t)

## [1] TRUE

is.numeric(my.solar.t)

## [1] TRUE

my.solar.d <- solar_time(Paris.time, geocode = Paris.geo, unit.out = "datetime")
is.solar_date(my.solar.d)

## [1] TRUE

is.timepoint(my.solar.d)

## [1] TRUE

Time of day

Function as_tod() facilitates conversion of R’s time date objects into a numeric value representing the time of day as numerical value with a decimal fraction using one of hour, minute or second as unit of expression. While solar time is based on the astronomical position of the sun, time of day is based on the time coordinates for a time zone.

times <- now(tzone = "UTC") + hours(0:6)
times

## [1] "2024-08-25 19:11:38 UTC" "2024-08-25 20:11:38 UTC"
## [3] "2024-08-25 21:11:38 UTC" "2024-08-25 22:11:38 UTC"
## [5] "2024-08-25 23:11:38 UTC" "2024-08-26 00:11:38 UTC"
## [7] "2024-08-26 01:11:38 UTC"

as_tod(times)

## [1] 19.1939852 20.1939852 21.1939852 22.1939852 23.1939852  0.1939852  1.1939852

as_tod(times, unit.out = "minutes")

## [1] 1151.63911 1211.63911 1271.63911 1331.63911 1391.63911   11.63911   71.63911

Relative air mass

Solar elevation determines the length of the path of the sun beam through the Earth’s atmosphere. This affects the solar spectrum at ground level, specially the UVB region. Function relative_AM() can be used to calculate an empirical estimate of this quantity from elevation angles in degrees. This function is vectorised.

relative_AM(33)

## [1] 1.83

relative_AM(c(80, 60, 40, 20))

## [1] 1.01 1.15 1.55 2.90

Water vapour-related calculations

There are several different empirical equations in common use to estimate water vapour pressure in air over liquid water or ice. This package implements some of them as well as their inverses, which are useful to compute the dew point. Some ancillary functions help with unit conversion and computation of relative humidity.

Saturation vapour pressure in Pa at four air temperatures in C, over water, computed using Tetens equation (the default method).

water_vp_sat(c(0, 10, 20, 30))

## [1]  610.780 1227.873 2338.023 4242.450

Saturation water vapour mass/volume concentration for air at 20 C, expressed in grammes per cubic meter.

water_vp2mvc(water_vp_sat(20), 20)

## [1] 17.27438

Vapour pressure for 50% relative humidity in air at a temperature of 20 C.

water_RH2vp(50, 20)

## [1] 1169.011

And the inverse computation, RH from vapour pressure in Pa and air temperature in C.

water_vp2RH(500, 40)

## [1] 6.780297

Dew point for air with a water vapour pressure of 3 kPa.

water_dp(3000)

## [1] 24.09167

Reference evapotranspiration

Evapotranspiration is a measure of the water flux between land and atmosphere expressed as millimeters of water per unit ground area and unit time. One millimeter of evapotranspiration or precipitation is equivalent to one litre per square meter. Evapotranspiration is composed of the evaporation from the soil surface and other wet surfaces plus transpiration from plants (water evaporated inside leaves and diffusing as vapour to the outside of leaves). The condensation of dew represents a negative component.

When not measured, evapotranspiration can be estimated based on energy balance and resistances to mass and heat transport. The Penman-Monteith equation is widely used. One special idealized condition is used to compute reference evapotranspiration (similar to potential evapotranspiration but with a standardized formulation by FAO that is widely used in agriculture).

As implemented the computation of the energy balance is done with function net_irradiance() and using the value returned by this function as one argument, reference evapotranspiration is computed with function ET_ref().

r_net <- net_irradiance(temperature = 20, # C
                        sw.down.irradiance = 800, # W / m2
                        water.vp = 1500) # Pa
r_net # W / m2

## [1] 546.8251

ET_ref(temperature = 20, # C
       water.vp = 1500, # Pa
       wind.speed = 4, # m / s
       net.irradiance = r_net
       )

## [1] 0.9448297

This function returns an estimate of reference evapotranspiration expressed in mm / h. Can be used with hourly means for the inputs or more frequent observations. In all cases the result is expressed as an instantaneous water flux rate expressed per hour.

Function ET_ref_day() should be used when the input data are daily meeans or less frequent.