---
title: "User Guide: 2 Astronomy and Atmosphere"
subtitle: "Package 'photobiology' `r packageVersion('photobiology')` "
author: "Pedro J. Aphalo"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: yes
vignette: >
  %\VignetteIndexEntry{User Guide: 2 Astronomy and Atmosphere}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

## 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.

```{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 '**photobiology**' and '**lubridate**', as they will be used in the examples.

```{r, pkg-load, eval=TRUE, message = FALSE}
library(photobiology)
library(lubridate)
```

## 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 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.

```{r}
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_.

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

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.

```{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 = "EET") + hours(c(0, 6, 12)), 
           geocode = my.geocode)
```

```{r}
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)
```

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)
```

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

```{r}
sun_angles(time = getWhenMeasured(sun.spct), geocode = getWhereMeasured(sun.spct))
```

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()`)

```{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

Functions `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.

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

```{r}
sunrise_time(now("UTC"), tz = "UTC", geocode = my.geocode)
sunrise_time(now("EET"), tz = "EET", geocode = my.geocode)
```

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

```{r}
sunrise_time(dates, geocode = data.frame(lat = 60, lon = 0))
sunrise_time(dates, geocode = data.frame(lat = -60, lon = 0))
```

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 = "EET"), 
             tz = "EET", 
             geocode = my.geocode,
             twilight = "civil")
sunrise_time(ymd("2017-03-21", tz = "EET"), 
             tz = "EET", 
             geocode = my.geocode,
             twilight = -10)
sunrise_time(ymd("2017-03-21", tz = "EET"), 
             tz = "EET", 
             geocode = my.geocode,
             twilight = +12)
```

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.

```{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:2], 
          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:2], 
          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)
```

Different convenience functions return the calculated variables individually as
vectors.

```{r}
sunrise_time(date = dates, geocode = my.geocode)
```

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).

```{r}
sunrise_time(date = dates, tz = "UTC", geocode = my.geocode)
```

```{r}
noon_time(date = dates, geocode = my.geocode)
```

```{r}
sunset_time(date = dates, geocode = my.geocode)
```

The default for `date` is the current day, using the system time zone settings.

```{r}
noon_time(geocode = my.geocode)
```

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 = "EET"), 
             tz = "EET", 
             geocode = my.geocode)
sunrise_time(ymd("2017-03-21", tz = "EET"), 
             tz = "EET", 
             geocode = my.geocode,
             unit.out = "hours")
```

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

```{r}
day_length(dates, geocode = my.geocode)
night_length(dates, geocode = my.geocode)
```

```{r}
day_length(dates, geocode = my.geocode, unit.out = "days")
night_length(dates, geocode = my.geocode, unit.out = "days")
```

## 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, 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.

```{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

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.

```{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
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.

```{r}
relative_AM(33)
```

```{r}
relative_AM(c(80, 60, 40, 20))
```

## 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.

| Function | parameters | returned value | methods implemented |
|:--------|:-------|:----------|:---------------------------|
| water_vp_sat() | temperature, C | vapour pressure, Pa | Tetens, Magnus, Wexler, Goff-Gratch |
| water_dp() | vapour pressure, Pa | temperature, C | Tetens, Magnus, Wexler |
| water_fp() | vapour pressure, Pa | temperature, C | Tetens, Magnus, Wexler |
| water_vp2RH() | vapour pressure, Pa; temperature, C | relative humidity, % | Tetens, Magnus, Wexler, Goff-Gratch |
| water_RH2vp() | relative humidity, %; temperature, C | vapour pressure, Pa | Tetens, Magnus, Wexler, Goff-Gratch |
| water_vp2mvc() | vapour pressure, Pa | mass/vol. concentration, g m$^{-3}$ | |
| water_mvc2vp() | mass/vol. concentration, g m$^{-3}$ | vapour pressure, Pa | |


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

```{r}
water_vp_sat(c(0, 10, 20, 30))
```

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

```{r}
water_vp2mvc(water_vp_sat(20), 20)
```

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

```{r}
water_RH2vp(50, 20)
```

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

```{r}
water_vp2RH(500, 40)
```

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

```{r}
water_dp(3000)
```

## 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}
r_net <- net_irradiance(temperature = 20, # C
                        sw.down.irradiance = 800, # W / m2
                        water.vp = 1500) # Pa
r_net # W / m2
```
```{r}
ET_ref(temperature = 20, # C
       water.vp = 1500, # Pa
       wind.speed = 4, # m / s
       net.irradiance = r_net
       )
```
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.


