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.
We load two packages, our ‘photobiology’ and ‘lubridate’, as they will be used in the examples.
library(photobiology)
library(lubridate)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.
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 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 x 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
## # ... with 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 x 12
## time tz solartime longitude latitude address azimuth
## <dttm> <chr> <solar_tm> <dbl> <dbl> <chr> <dbl>
## 1 2022-03-24 22:18:37 UTC 22:12:30 0 51.5 Greenwich 328.
## # ... with 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 x 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.
## # ... with 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 x 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.
## # ... with 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 x 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.
## # ... with 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 x 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.
## # ... with 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.710245sun_zenith_angle(time = my.times, geocode = my.geocode)## [1] 142.63935 112.72250 83.28976sun_azimuth(time = my.times, geocode = my.geocode)## [1] 351.04757 96.98377 174.48767Functions sunrise_time(), sunset_time(),
noon_time(), day_length() and
night_length() have all the same parameter signature. An
additional function, day_night returns a data frame
containing all the quantities returned by the other functions. They are
all vectorized for the date and geocode
parameters. As arguments are the same for all these functions, we use
sunrise_time() in the examples below, but these examples
are equally relevant to the other functions listed at the start of this
paragraph.
Both latitude and longitude should be supplied through a
geocode, but be aware that 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. We reuse the geocode data frames created
above, and create a vector of datetime objects differing in their date.
The default time zone of function ymd is UTC, 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"sunrise_time(now("UTC"), tz = "UTC", geocode = my.geocode)## [1] "2022-03-24 04:09:12 UTC"sunrise_time(now("EET"), tz = "EET", geocode = my.geocode)## [1] "2022-03-24 04:09:12 EET"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 = 0))## [1] "2015-02-28 07:01:52 EET" "2015-03-31 05:28:41 EEST"
## [3] "2015-04-30 04:00:50 EEST" "2015-05-31 02:50:49 EEST"
## [5] "2015-06-30 02:41:10 EEST" "2015-07-31 03:38:15 EEST"sunrise_time(dates, geocode = data.frame(lat = -60, lon = 0))## [1] "2015-02-28 05:09:25 EET" "2015-03-31 06:26:18 EEST"
## [3] "2015-04-30 07:38:10 EEST" "2015-05-31 08:44:25 EEST"
## [5] "2015-06-30 09:04:25 EEST" "2015-07-31 08:17:27 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 03:38:58 EET"sunrise_time(ymd("2017-03-21", tz = "EET"),
tz = "EET",
geocode = my.geocode,
twilight = -10)## [1] "2017-03-20 03:05:45 EET"sunrise_time(ymd("2017-03-21", tz = "EET"),
tz = "EET",
geocode = my.geocode,
twilight = +12)## [1] "2017-03-20 06:06:15 EET"Default latitude is zero (the Equator), the default longitude is zero (Greenwich). Be also aware that for summer dates the times are formatted for printing accordingly. In the examples below this can be recognized by the time zone being reported as EEST instead of EET during the summer for Eastern Europe.
The main function is called day_night() and returns a
data frame containing both the input values and the results of the
calculations. See below for additional 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 x 12
## day tz twilight.rise twilight.set longitude latitude address sunrise
## <date> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 2015-02-28 UTC 0 0 24.9 60.2 Helsin~ 5.49
## 2 2015-03-31 UTC 0 0 24.9 60.2 Helsin~ 3.93
## 3 2015-04-30 UTC 0 0 24.9 60.2 Helsin~ 2.47
## # ... with 4 more variables: 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:2],
geocode = my.geocode,
unit.out = "days")## # A tibble: 2 x 12
## day tz twilight.rise twilight.set longitude latitude address sunrise
## <date> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 2015-02-28 UTC 0 0 24.9 60.2 Helsin~ 0.229
## 2 2015-03-31 UTC 0 0 24.9 60.2 Helsin~ 0.164
## # ... with 4 more variables: 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:2],
geocode = my.geocode,
unit.out = "datetime")## # A tibble: 2 x 12
## day tz twilight.rise twilight.set longitude latitude address
## <date> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 2015-02-28 UTC 0 0 24.9 60.2 Helsinki
## 2 2015-03-31 UTC 0 0 24.9 60.2 Helsinki
## # ... with 5 more variables: 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 x 12
## day tz twilight.rise twilight.set longitude latitude address sunrise
## <date> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 2015-02-28 UTC 0 0 24.9 60.2 Helsin~ 5.49
## 2 2015-03-31 UTC 0 0 24.9 60.2 Helsin~ 3.93
## 3 2015-04-30 UTC 0 0 24.9 60.2 Helsin~ 2.47
## 4 2015-02-28 UTC 0 0 25.5 65.0 Oulu 5.68
## 5 2015-03-31 UTC 0 0 25.5 65.0 Oulu 3.78
## 6 2015-04-30 UTC 0 0 25.5 65.0 Oulu 1.96
## # ... with 4 more variables: noon <dbl>, sunset <dbl>, daylength <dbl>,
## # nightlength <dbl>Different convenience functions return the calculated variables individually as vectors.
sunrise_time(date = dates, geocode = my.geocode)## [1] "2015-02-28 05:22:28 EET" "2015-03-31 03:48:44 EEST"
## [3] "2015-04-30 02:20:18 EEST" "2015-05-31 01:09:32 EEST"
## [5] "2015-06-30 00:59:39 EEST" "2015-07-31 01:57:25 EEST"As seen above the default for tz is the time zone of the
argument passed to date. This can be overridden with an
explicit value as argument. In this example, when interpreted in the UTC
time zone, the time instants correspond to the previous calendar day
compared to the EET time zone. We can also see that “summer time”
applies to the EET time zone but not to UTC (universal time
coordinates).
sunrise_time(date = dates, tz = "UTC", geocode = my.geocode)## [1] "2015-02-28 05:22:28 UTC" "2015-03-31 03:48:44 UTC"
## [3] "2015-04-30 02:20:18 UTC" "2015-05-31 01:09:32 UTC"
## [5] "2015-06-30 00:59:39 UTC" "2015-07-31 01:57:25 UTC"noon_time(date = dates, geocode = my.geocode)## [1] "2015-02-28 10:32:51 EET" "2015-03-31 10:24:30 EEST"
## [3] "2015-04-30 10:17:32 EEST" "2015-05-31 10:17:55 EEST"
## [5] "2015-06-30 10:23:53 EEST" "2015-07-31 10:26:41 EEST"sunset_time(date = dates, geocode = my.geocode)## [1] "2015-02-28 15:43:14 EET" "2015-03-31 17:00:17 EEST"
## [3] "2015-04-30 18:14:46 EEST" "2015-05-31 19:26:19 EEST"
## [5] "2015-06-30 19:48:08 EEST" "2015-07-31 18:55:57 EEST"The default for date is the current day, using the
system time zone settings.
noon_time(geocode = my.geocode)## [1] "2022-03-25 10:26:13 UTC"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 04:20:46 EET"sunrise_time(ymd("2017-03-21", tz = "EET"),
tz = "EET",
geocode = my.geocode,
unit.out = "hours")## [1] 6.346365Functions day_length() and night_length()
return by default the length of time in hours.
day_length(dates, geocode = my.geocode)## [1] 10.34596 13.19241 15.90766 18.27962 18.80811 16.97567night_length(dates, geocode = my.geocode)## [1] 13.654040 10.807592 8.092343 5.720384 5.191888 7.024327day_length(dates, geocode = my.geocode, unit.out = "days")## [1] 0.4310817 0.5496837 0.6628190 0.7616507 0.7836713 0.7073197night_length(dates, geocode = my.geocode, unit.out = "days")## [1] 0.5689183 0.4503163 0.3371810 0.2383493 0.2163287 0.2926803In field research it is in many cases preferable to sample or measure, and present and plot data based on local solar time, as this are the time coordinates organisms are exposed to in nature. Two functions are provided. They differ in the value returned, either a time of day in hours, or a datetime.
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] TRUEis.numeric(my.solar.t)## [1] TRUEmy.solar.d <- solar_time(Paris.time, geocode = Paris.geo, unit.out = "datetime")
is.solar_date(my.solar.d)## [1] TRUEis.timepoint(my.solar.d)## [1] TRUEFunction 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.
times <- now(tzone = "UTC") + hours(0:6)
times## [1] "2022-03-24 22:18:39 UTC" "2022-03-24 23:18:39 UTC"
## [3] "2022-03-25 00:18:39 UTC" "2022-03-25 01:18:39 UTC"
## [5] "2022-03-25 02:18:39 UTC" "2022-03-25 03:18:39 UTC"
## [7] "2022-03-25 04:18:39 UTC"as_tod(times)## [1] 22.3108437 23.3108437 0.3108437 1.3108437 2.3108437 3.3108437 4.3108437as_tod(times, unit.out = "minutes")## [1] 1338.65062 1398.65062 18.65062 78.65062 138.65062 198.65062 258.65062Solar 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.83relative_AM(c(80, 60, 40, 20))## [1] 1.01 1.15 1.55 2.90Evapotranspiration 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.8251ET_ref(temperature = 20, # C
water.vp = 1500, # Pa
wind.speed = 4, # m / s
net.irradiance = r_net
)## [1] 0.9448297This 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.