# User Guide: 3 Astronomy

### Package ‘photobiology’ 0.10.0

#### Pedro J. Aphalo

#### 2020-04-02

## Getting started

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

    library(photobiology)
    library(lubridate)

## Introduction

The functions and methods described in this section 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.

## Position of the sun

In photobiology research we sometimes need to calculate the position on the sun at arbitrary geographic locations and times. The function `sun_angles` returns the azimuth in degrees eastwards, altitude in degrees above the horizon, solar disk diameter in degrees and sun to earth distance in astronomical units. The time should be a `POSIXct` vector, possibly of length one. The easiest way create date and time constant values is to use package `lubridate`.

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 object to be supplied as argument for `geocode` is a data frame with variables `lon` and `lat` giving the location of Earth. This matches the return value of function `ggmap::geocode()`, function that can be used to find the coordinates using any address entered as a character string understood by the Google maps API. We use the “geocode” for Helsinki.

Be aware that to obtain results based on local time, the correct time zone needs to be set. 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`.

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

The position of the sun at Helsinki, at the given time instant decoded 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 11
    ##   time                tz    solartime longitude latitude address azimuth
    ##   <dttm>              <chr> <solar_t>     <dbl>    <dbl> <chr>     <dbl>
    ## 1 2017-06-20 08:00:00 EET   06:38:09       24.9     60.2 Helsin~    85.8
    ## # ... with 4 more variables: elevation <dbl>, declination <dbl>,
    ## #   eq.of.time <dbl>, hour.angle <dbl>

Functions have defaults for their arguments, but rarely Greenwich will be the location you are interested in.

    sun_angles()

    ## # A tibble: 1 x 11
    ##   time                tz    solartime longitude latitude address azimuth
    ##   <dttm>              <chr> <solar_t>     <dbl>    <dbl> <chr>     <dbl>
    ## 1 2020-04-02 09:54:58 UTC   09:51:32          0     51.5 Greenw~    139.
    ## # ... with 4 more variables: elevation <dbl>, declination <dbl>,
    ## #   eq.of.time <dbl>, hour.angle <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 11
    ##   time                tz    solartime longitude latitude address azimuth
    ##   <dttm>              <chr> <solar_t>     <dbl>    <dbl> <chr>     <dbl>
    ## 1 2014-01-01 00:00:00 EET   23:36:26       24.9     60.2 Helsin~   351.
    ## 2 2014-01-01 06:00:00 EET   05:36:19       24.9     60.2 Helsin~    97.0
    ## 3 2014-01-01 12:00:00 EET   11:36:12       24.9     60.2 Helsin~   174.
    ## # ... with 4 more variables: elevation <dbl>, declination <dbl>,
    ## #   eq.of.time <dbl>, hour.angle <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 11
    ##   time                tz    solartime longitude latitude address azimuth
    ##   <dttm>              <chr> <solar_t>     <dbl>    <dbl> <chr>     <dbl>
    ## 1 2014-01-01 00:00:00 EET   23:36:26       24.9     60.2 Helsin~   351.
    ## 2 2014-01-01 06:00:00 EET   05:36:19       24.9     60.2 Helsin~    97.0
    ## 3 2014-01-01 12:00:00 EET   11:36:12       24.9     60.2 Helsin~   174.
    ## # ... with 4 more variables: elevation <dbl>, declination <dbl>,
    ## #   eq.of.time <dbl>, hour.angle <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 11
    ##   time                tz    solartime longitude latitude address azimuth
    ##   <dttm>              <chr> <solar_t>     <dbl>    <dbl> <chr>     <dbl>
    ## 1 2014-01-01 00:00:00 EET   23:36:26       24.9     60.2 Helsin~   351.
    ## 2 2014-01-01 06:00:00 EET   05:36:19       24.9     60.2 Helsin~    97.0
    ## 3 2014-01-01 12:00:00 EET   11:36:12       24.9     60.2 Helsin~   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 4 more variables: elevation <dbl>, declination <dbl>,
    ## #   eq.of.time <dbl>, hour.angle <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 11
    ##   time                tz    solartime longitude latitude address azimuth
    ##   <dttm>              <chr> <solar_t>     <dbl>    <dbl> <chr>     <dbl>
    ## 1 2010-06-22 09:51:00 UTC   11:28:54       25.0     60.2 Kumpul~    168.
    ## # ... with 4 more variables: elevation <dbl>, declination <dbl>,
    ## #   eq.of.time <dbl>, hour.angle <dbl>

One what is needed is only one of the solar angles, functions returning vectors instead of data frames can be useful.

    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

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 they apply to the other functions described in this section.

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] "2020-04-02 03:40:15 UTC"

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

    ## [1] "2020-04-02 06:40:15 EEST"

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 09:01:52 EET"  "2015-03-31 08:28:41 EEST"
    ## [3] "2015-04-30 07:00:50 EEST" "2015-05-31 05:50:49 EEST"
    ## [5] "2015-06-30 05:41:10 EEST" "2015-07-31 06:38:15 EEST"

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

    ## [1] "2015-02-28 07:09:25 EET"  "2015-03-31 09:26:18 EEST"
    ## [3] "2015-04-30 10:38:10 EEST" "2015-05-31 11:44:25 EEST"
    ## [5] "2015-06-30 12:04:25 EEST" "2015-07-31 11: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 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"

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 Helsin~
    ## 2 2015-03-31 UTC               0            0      24.9     60.2 Helsin~
    ## # ... 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 07:22:28 EET"  "2015-03-31 06:48:44 EEST"
    ## [3] "2015-04-30 05:20:18 EEST" "2015-05-31 04:09:32 EEST"
    ## [5] "2015-06-30 03:59:39 EEST" "2015-07-31 04: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 12:32:51 EET"  "2015-03-31 13:24:30 EEST"
    ## [3] "2015-04-30 13:17:32 EEST" "2015-05-31 13:17:55 EEST"
    ## [5] "2015-06-30 13:23:53 EEST" "2015-07-31 13:26:41 EEST"

    sunset_time(date = dates, geocode = my.geocode)

    ## [1] "2015-02-28 17:43:14 EET"  "2015-03-31 20:00:17 EEST"
    ## [3] "2015-04-30 21:14:46 EEST" "2015-05-31 22:26:19 EEST"
    ## [5] "2015-06-30 22:48:08 EEST" "2015-07-31 21:55:57 EEST"

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

    noon_time(geocode = my.geocode)

    ## [1] "2020-04-02 10:23:41 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 06:20:46 EET"

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

    ## [1] 6.346365

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

    night_length(dates, geocode = my.geocode)

    ## [1] 13.654040 10.807592  8.092343  5.720384  5.191888  7.024327

    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

## 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. 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] 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 in one of hour, minute or second as unit of expression.

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

    ## [1] "2020-04-02 09:54:59 UTC" "2020-04-02 10:54:59 UTC"
    ## [3] "2020-04-02 11:54:59 UTC" "2020-04-02 12:54:59 UTC"
    ## [5] "2020-04-02 13:54:59 UTC" "2020-04-02 14:54:59 UTC"
    ## [7] "2020-04-02 15:54:59 UTC"

    as_tod(times)

    ## [1]  9.916555 10.916555 11.916555 12.916555 13.916555 14.916555 15.916555

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

    ## [1] 594.9933 654.9933 714.9933 774.9933 834.9933 894.9933 954.9933

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

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

    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