---
title: "User Guide"
subtitle: "Package 'SunCalcMeeus' `r packageVersion('SunCalcMeeus')` "
author: "Pedro J. Aphalo"
date: "`r packageDate('SunCalcMeeus')`"
output:
rmarkdown::html_vignette:
toc: yes
vignette: >
%\VignetteIndexEntry{User Guide}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
## Radiation, astronomy and atmosphere
The functions described here used to be part of package 'photobiology',
but as of version 0.11.4 have been moved to this package. To ensure
backwards compatibility, package 'photobiology' will depend on
'SunCalcMeeus' once these functions are removed.
```{r, include=FALSE, echo=FALSE}
knitr::opts_chunk$set(fig.width=8, fig.height=4)
```
```{r, printing-spectra, eval=TRUE, include=FALSE}
# library(tibble)
options(tibble.print_max = 6, tibble.print_min = 4)
```
## Getting started
We load two packages, our '**SunCalcMeeus**' and '**lubridate**', as they will be used in the examples.
```{r, pkg-load, eval=TRUE, message = FALSE}
library(SunCalcMeeus)
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.
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.
```{r}
my.geocode <- data.frame(lat = 60.16, lon = 24.93, address = "Helsinki")
```
Be aware that to obtain correct computed values 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` with one row per time instant or
per geographic location.
The position of the sun at Helsinki, at the given instant in time for
time zone `"Europe/Helsinki"`, which matches _Eastern European Time_.
```{r}
sun_angles(time = ymd_hms("2017-06-20 08:00:00", tz = "Europe/Helsinki"),
geocode = my.geocode)
```
Functions have defaults for their arguments, Greenwhich in U.K. and the
corresponding time zone "UTC". In most cases Greenwich will not be the location
you are interested in. Current UTC time is more likely to be a useful default as
it avoids the difficulty of time shifts in local time coordinates.
```{r}
sun_angles()
```
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.
```{r}
sun_angles(time = ymd_hms("2014-01-01 0:0:0", tz = "Europe/Helsinki") + hours(c(0, 6, 12)),
geocode = my.geocode)
```
```{r}
my.times <- ymd_hms("2014-01-01 0:0:0", tz = "Europe/Helsinki") + hours(c(0, 6, 12))
sun_angles(time = my.times, geocode = my.geocode)
```
Geocodes for several locations can be stored in a data frame with multiple rows.
```{r}
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)
```
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()`_. Thus
if more than one angle is needed, it is more efficient to compute all angles
with function `sun_angles()` and later extract the vectors from the returned
data frame.
```{r}
sun_elevation(time = my.times, geocode = my.geocode)
```
```{r}
sun_zenith_angle(time = my.times, geocode = my.geocode)
```
```{r}
sun_azimuth(time = my.times, geocode = my.geocode)
```
## 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 these
other functions.
These functions are vectorized for their `date` and `geocode` parameters. They use as default location Greenwich in the U.K., and corresponding 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 "Europe/Helsinki" (matching Eastern European Time) and UTC is +2 hours in winter and (matching Eastern European Summer Time) +3 hours in summer.
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 `Europe/Helsinki` to match the geocodes for Finnish cities.
```{r}
dates <- ymd("2015-03-01", tz = "Europe/Helsinki") + months(0:5)
dates
```
As first example we compute the sunrise time for the current day in Helsinki, with the result returned either in UTC or local time coordinates. Time-zone names based on continent and city ("Europe/Helsinki") or continent, country and city ("America/Argentina/Buenos Aires") are supported, while the names "EET" and "CET" and their summer-time versions are no longer supported by R (>= 4.5.0). They have been long deprecated as they do not describe true time zones, and different are within these regions have been in different time zones in the past, making it impossible some computations. Dates and the relationship between time zones and locations have been affected by changes in country boundaries and in national laws.
Use of the Olson time zone names like `"Europe/Helsinki"` is recommended. The list is available in R and can be searched.
```{r}
grep("Argentina", OlsonNames(), value = TRUE)
```
The time zone in use by the computer on which R is running can be found out with the following code.
```{r}
Sys.timezone()
```
At least in R (< 4.5.0) the "EET", "CET", etc. are still used when printing or formatting the output. When not needed, the time zone abbreviation can be disabled in printing and formatting.
```{r}
# defaults to current UTC date and Greenwich, UK as location
sunrise_time()
sunrise_time(date = now(), tz = "Europe/Helsinki", geocode = my.geocode)
sunrise_time(date = now(tzone = "Europe/Helsinki"), geocode = my.geocode)
# time zone abbreviation not shown
print(sunrise_time(date = now("Europe/Helsinki"),
geocode = my.geocode),
usetz = FALSE)
```
----
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 'SunCalcMeeus'. 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 'SunCalcMeeus'. _An instant in time and time zone strictly define a corresponding date at any location on Earth, even though the date is not the same at all these locations._
The time zone used by default for the returned value is that of the `POSIXct` value passed as argument to parameter `date`. This behaviour can be overridden by an argument passed to `tz`. However, to obtain a correct value expressed in local time We must make sure that the time zone matches that at the `geocode` location.
```{r}
## Using date times as POSIXct
sunrise_time(geocode = my.geocode)
sunrise_time(date = now("UTC"), geocode = my.geocode)
sunrise_time(date = now("UTC"), tz = "UTC", geocode = my.geocode)
sunrise_time(date = now("Europe/Helsinki"), geocode = my.geocode)
sunrise_time(date = now(""), tz = "Europe/Helsinki", geocode = my.geocode)
## Using Date
# correct always as time zones match
sunrise_time(today("Europe/Helsinki"), tz = "Europe/Helsinki", geocode = my.geocode)
# sometimes the value returned will be correct and sometimes off by 1 d at Helsinki
sunrise_time(today("Australia/Canberra"), tz = "Europe/Helsinki", geocode = my.geocode)
```
----
We can also compute the time at solar noon and at sunset.
```{r}
noon_time(now("UTC"), geocode = my.geocode)
noon_time(now("Europe/Helsinki"), geocode = my.geocode)
```
By default, sunset and sunrise are defined as the time when the upper rim of the solar disk is at the horizon. How to override this default to account for twilight and or obstacles that occlude the sun will be shown later.
```{r}
sunset_time(now("UTC"), geocode = my.geocode)
sunset_time(now("Europe/Helsinki"), geocode = my.geocode)
sunrise_time(now("Europe/Helsinki"), geocode = my.geocode)
```
Functions `day_length()` and `night_length()` return the length of time between sunrise and sunset and between sunset and sunrise, respectively, by default expressed in hours and fraction.
```{r}
day_length(dates, geocode = my.geocode)
night_length(dates, geocode = my.geocode)
day_length(dates, geocode = my.geocode, unit.out = "day")
```
Southern hemisphere latitudes as well as longitudes to the West of the Greenwich
meridian should be supplied as negative numbers. In this case the time zone is
abbreviated as a time difference from UTC.
```{r}
sunrise_time(dates, tz = "America/Argentina/Buenos_Aires",
geocode = data.frame(lat = -34.6, lon = -58.3))
noon_time(dates, tz = "America/Argentina/Buenos_Aires",
geocode = data.frame(lat = -34.6, lon = -58.3))
```
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.
```{r}
sunrise_time(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = "none") # center of the sun disk
sunrise_time(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = "sunlight") # upper rim of the sun disk
sunrise_time(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = "rim") # lower rim of the sun disk
sunrise_time(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = "civil") # civil twilight = -6 degrees
sunrise_time(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = -10) # 10 degrees below the horizon
sunrise_time(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = +12) # 12 degrees above the horizon
```
Twilight is also relevant to the computation of day length and night length. The default is to use the centre of the sun disk, but this can be changed. For the values returned by `day_length()` and `night_length()` to add to 24 h they must be computed using the same twilight definition.
```{r}
day_length(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode)
day_length(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = 0)
day_length(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = "rim")
day_length(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = "civil")
day_length(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = -10)
```
Say if there is a mountain blocking the view above the Western horizon, we can set different twilight angles for sunrise and sunset.
```{r}
day_length(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = c(0, 12))
night_length(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
twilight = c(0, 12))
```
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`.
```{r}
sunrise_time(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode)
sunrise_time(ymd("2017-03-21", tz = "Europe/Helsinki"),
tz = "Europe/Helsinki",
geocode = my.geocode,
unit.out = "hours")
```
Similarly, the units can also be selected for the values returned by `day_length()` and `night_length()`.
```{r}
day_length(dates, geocode = my.geocode, unit.out = "days")
night_length(dates, geocode = my.geocode, unit.out = "days")
```
----
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.
```{r}
day_night(dates[1:3],
geocode = my.geocode)
```
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.
```{r}
day_night(dates[1:3],
geocode = my.geocode,
unit.out = "days")
```
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.
```{r}
day_night(dates[1:3],
geocode = my.geocode,
unit.out = "datetime")
```
When multiple times and locations are supplied as arguments, the values returned
are for all combinations of locations and times.
```{r}
day_night(dates[1:3],
geocode = two.geocodes)
```
## 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 'SunCalcMeeus', with `print()` and `format()` 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.
```{r}
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)
solar_time(Paris.time, geocode = Paris.geo, unit.out = "datetime")
```
```{r}
my.solar.t <- solar_time(Paris.time, geocode = Paris.geo)
is.solar_time(my.solar.t)
is.numeric(my.solar.t)
```
```{r}
my.solar.d <- solar_time(Paris.time, geocode = Paris.geo, unit.out = "datetime")
is.solar_date(my.solar.d)
is.timepoint(my.solar.d)
```
## Time of day
When analysing data as a time series the usual way to represent time is as a
date plus time value, i.e., as an instant in time. In contrast, when data need
to summarised or plotted as a function of time of day, the date portion of a
data time representation of time becomes a nuisance.
Function `as_tod()` facilitates conversion of R's time date objects into values
representing the time of day as numerical value giving the time elapsed since
the most recent past midnight. This value can be represented as a numeric value
using `"day"`, `"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.
```{r}
times <- now(tzone = "UTC") + hours(0:6)
times
as_tod(times)
as_tod(times, unit.out = "minutes")
```
## 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
in the UV-B region. Function `relative_AM()` can be used to calculate an
empirical _estimate_ of this quantity from elevation angles in degrees. This
function is vectorised. As seem above the apparent position of the sun for an
observer on Earth differs from the true astronomical position as a result of
atmospheric refraction. The functions presented above, can apply a correction
based on an estimate of atmospheric refraction. If the correction is applied we
obtain apparent sun elevation angles suitable for estimating the relative air
mass with Kasten and Youngs' (1989) equation as follows.
```{r}
relative_AM(33)
```
```{r}
relative_AM(c(90, 60, 40, 20, 10, 5, 2, 1, 0.5))
```
Young's (1994) equation can be used to estimate AM from the true sun elevation.
At high solar elevations the impact of atmospheric refraction is very small
but at low elevations is is large enough to make a clear difference in the
AM estimates.
```{r}
relative_AMt(33)
```
```{r}
relative_AMt(c(90, 60, 40, 20, 10, 5, 2, 1, 0.5))
```
Two additional functions make it possible to compute the relative AM from time and geographic coordinates. They both apply a correction for atmospheric refraction, but using two different algorithms in the sun elevation computation or in Young's (1994) equation. The difference is small. These convenience functions are wrappers on `relative_AM()` and `relative_AMt()` that call function `sun_elevation()` to obtain the input to pass to the wrapped functions.
```{r}
january.times <- ymd_h("2020-01-01 12", tz = "Europe/Helsinki") + hours(-2:+2)
relative_AM_geotime(january.times, my.geocode, tz = "Europe/Helsinki")
```
```{r}
relative_AMt_geotime(january.times, my.geocode, tz = "Europe/Helsinki")
```