Nothing
#' Solar angles
#'
#' Function \code{sun_angles()} returns the solar angles and Sun to Earth
#' relative distance for given times and locations using a very precise
#' algorithm. Convenience functions \code{sun_azimuth()},
#' \code{sun_elevation()}, \code{sun_zenith_angle()} and
#' \code{distance_to_sun()} are wrappers on \code{sun_angles()} that return
#' individual vectors.
#'
#' @param time A "vector" of POSIXct Time, with any valid time zone (TZ) is
#' allowed, default is current time.
#' @param tz character string indicating time zone to be used in output.
#' @param geocode data frame with variables lon and lat as numeric values
#' (degrees), nrow > 1, allowed.
#' @param use.refraction logical Flag indicating whether to correct for
#' fraction in the atmosphere.
#'
#' @return A \code{data.frame} with variables \code{time} (in same TZ as input),
#' \code{TZ}, \code{solartime}, \code{longitude}, \code{latitude},
#' \code{address}, \code{azimuth}, \code{elevation}, \code{declination},
#' \code{eq.of.time}, \code{hour.angle}, and \code{distance}. If a data frame
#' with multiple rows is passed to \code{geocode} and a vector of times longer
#' than one is passed to \code{time}, sun position for all combinations of
#' locations and times are returned by \code{sun_angles}. Angles are expressed
#' in degrees, \code{solartime} is a vector of class \code{"solar.time"},
#' \code{distance} is expressed in relative sun units.
#'
#' @family astronomy related functions
#'
#' @details This function is an implementation of Meeus equations as used in
#' NOAAs on-line web calculator, which are precise and valid for a very broad
#' range of dates (years -1000 to 3000 at least). The apparent solar
#' elevations near sunrise and sunset are affected by refraction in the
#' atmosphere, which does in turn depend on weather conditions. The effect of
#' refraction on the apparent position of the sun is only an estimate based on
#' "typical" conditions for the atmosphere. The computation is not defined for
#' latitudes 90 and -90 degrees, i.e. exactly at the poles. The function is
#' vectorized and in particular passing a vector of times for a single geocode
#' enhances performance very much as the equation of time, the most time
#' consuming step, is computed only once.
#'
#' For improved performance, if more than one angle is needed it
#' is preferable to directly call \code{sun_angles} instead of the wrapper
#' functions as this avoids the unnecesary recalculation.
#'
#' @section Important!: Given an instant in time and a time zone, the date is
#' computed from these, and may differ by one day to that at the location
#' pointed by \code{geocode} at the same instant in time, unless the argument
#' passed to \code{tz} matches the time zone at this location.
#'
#' @note There exists a different R implementation of the same algorithms called
#' "AstroCalcPureR" available as function \code{astrocalc4r} in package
#' 'fishmethods'. Although the equations used are almost all the same, the
#' function signatures and which values are returned differ. In particular,
#' the present implementation splits the calculation into two separate
#' functions, one returning angles at given instants in time, and a separate
#' one returning the timing of events for given dates.
#'
#' @references
#' The primary source for the algorithm used is the book:
#' Meeus, J. (1998) Astronomical Algorithms, 2 ed., Willmann-Bell, Richmond,
#' VA, USA. ISBN 978-0943396613.
#'
#' A different implementation is available at
#' \url{https://github.com/NEFSC/READ-PDB-AstroCalc4R/}.
#'
#' An interactive web page using the same algorithms is available at
#' \url{https://gml.noaa.gov/grad/solcalc/}. There are small
#' differences in the returned times compared to our function that seem to be
#' related to the estimation of atmospheric refraction (about 0.1 degrees).
#'
#' @export
#'
#' @examples
#' library(lubridate)
#' sun_angles()
#' sun_azimuth()
#' sun_elevation()
#' sun_zenith_angle()
#' sun_angles(ymd_hms("2014-09-23 12:00:00"))
#' sun_angles(ymd_hms("2014-09-23 12:00:00"),
#' geocode = data.frame(lat=60, lon=0))
#' sun_angles(ymd_hms("2014-09-23 12:00:00") + minutes((0:6) * 10))
#'
sun_angles <- function(time = lubridate::now(tzone = "UTC"),
tz = lubridate::tz(time),
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
use.refraction = FALSE)
{
geocode <- validate_geocode(geocode)
stopifnot(lubridate::is.POSIXct(time))
stopifnot(length(tz) == 1)
z <- list(nrow(geocode))
for (i in seq_len(nrow(geocode))) {
temp <- sun_angles_fast(time = time,
tz = tz,
geocode = dplyr::slice(geocode, i),
use.refraction = use.refraction)
z[[i]] <- temp # needed so that class attribute is retained
}
# we supress warning of dropped attributes and restore them
z <- suppressWarnings(dplyr::bind_rows(z))
class(z[["solartime"]]) <- class(temp[["solartime"]])
# we could use rbind instead of dplyr::bind_rows as the second drops the class attribute
# for solartime.
# first.iter <- TRUE
# for (i in 1:nrow(geocode)) {
# temp <- sun_angles_fast(time = time,
# tz = tz,
# geocode = geocode[i, ],
# use.refraction = use.refraction)
# if (first.iter) {
# z <- temp
# first.iter <- FALSE
# } else {
# z <- rbind(z, temp,
# make.row.names = FALSE)
# }
# }
# assertion
if (any(z[["elevation"]] < (-90)) || any(z[["elevation"]] > 90))
warning("Returned 'elevation' value(s) off range")
if (any(!is.na(z[["azimuth"]]) & (z[["azimuth"]] < -1e-10 |
z[["azimuth"]] > (360 + 1e-10)))) {
warning("Returned 'azimuth' values(s) off range")
}
z
}
#' @rdname sun_angles
#'
# Internal function, called by sun_angles()
#
sun_angles_fast <- function(time,
tz,
geocode,
use.refraction)
{
# Input validation done in sun_angles() before calling this function.
# stopifnot(!anyNA(time))
# stopifnot(is.data.frame(geocode))
# stopifnot(nrow(geocode == 1) && length(tz == 1))
# We have a single geocode and all times are expressed in the same time zone!
# If time is a vector we can vectorize the whole calculation, and do the
# expensive calculations only once.
lon <- geocode[["lon"]]
lat <- geocode[["lat"]]
address <- geocode[["address"]]
cent <- julian_century(time)
sun.lon.mean <- geom_mean_lon_sun(cent)
sun.anom.mean <- geom_mean_anom_sun(cent)
eccent.earth <- eccent_earth_orbit(cent)
delta <- sun_eq_of_ctr(cent, sun.anom.mean)
sun.lon <- sun.lon.mean + delta
sun.anom <- sun.anom.mean + delta
sun.dist <- sun_rad_vector(eccent.earth, sun.anom)
sun.app.lon <- sun_app_lon(cent, sun.lon)
sun.ecliptic <- mean_obliq_eclip(cent)
obliq.corr <- obliq_corr(cent, sun.ecliptic)
# rt.ascen <- sun_rt_ascen(sun.app.lon, obliq.corr)
var.y <- var_y(obliq.corr)
eq.of.time <- eq_of_time(mean.lon = sun.lon.mean,
eccent.earth = eccent.earth,
anom.mean = sun.anom.mean,
var.y = var.y)
solar.time <- solar_tod(time, lat, lon, eq.of.time)
hour.angle <- hour_angle(solar.time)
sun.declin <- sun_decline(sun.app.lon, obliq.corr)
zenith.angle <- zenith_angle(lat, hour.angle, sun.declin)
elevation.angle <- 90 - zenith.angle
if (use.refraction) {
elevation.angle <-
elevation.angle + atm_refraction_approx(elevation.angle)
}
azimuth.angle <- azimuth_angle(lat, hour.angle, zenith.angle, sun.declin)
solar.time <- solar.time / 60 # hours
solar.time <- solar.time %% 24 # needed for DST
class(solar.time) <- c("solar_time", class(solar.time))
tibble::tibble(time = lubridate::with_tz(time, tz),
tz = rep(tz, length(time)),
solartime = solar.time,
longitude = rep(lon, length(time)),
latitude = rep(lat, length(time)),
address = rep(address, length(time)),
azimuth = azimuth.angle,
elevation = elevation.angle,
declination = sun.declin,
eq.of.time = eq.of.time,
hour.angle = hour.angle,
distance = sun.dist,
.name_repair = "minimal")
}
#' @rdname sun_angles
#'
#' @export
#'
sun_elevation <- function(time = lubridate::now(),
tz = lubridate::tz(time),
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
use.refraction = FALSE)
{
stopifnot(length(time) == 1 || nrow(geocode) == 1)
sun_angles(time = time,
tz = tz,
geocode = geocode,
use.refraction = use.refraction)[["elevation"]]
}
#' @rdname sun_angles
#'
#' @export
#'
sun_zenith_angle <- function(time = lubridate::now(),
tz = lubridate::tz(time),
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
use.refraction = FALSE)
{
stopifnot(length(time) == 1 || nrow(geocode) == 1)
90 - sun_angles(time = time,
tz = tz,
geocode = geocode,
use.refraction = use.refraction)[["elevation"]]
}
#' @rdname sun_angles
#'
#' @export
#'
sun_azimuth <- function(time = lubridate::now(),
tz = lubridate::tz(time),
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
use.refraction = FALSE)
{
stopifnot(length(time) == 1 || nrow(geocode) == 1)
sun_angles(time = time,
tz = tz,
geocode = geocode,
use.refraction = use.refraction)[["azimuth"]]
}
#' @rdname sun_angles
#'
#' @export
#'
distance_to_sun <- function(time = lubridate::now(),
tz = lubridate::tz(time),
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
use.refraction = FALSE)
{
stopifnot(length(time) == 1 || nrow(geocode) == 1)
90 - sun_angles(time = time,
tz = tz,
geocode = geocode,
use.refraction = use.refraction)[["distance"]]
}
#' Time difference between two time zones
#'
#' Returns the difference in local time expressed in hours between two time
#' zones at a given instant in time. The difference due to daylight saving time
#' or Summer and Winter time as well as historical changes in time zones are
#' taken into account.
#'
#' @note This function is implemented using functions from package 'lubridate'.
#' For details on the handling of time zones, please, consult the
#' documentation for \code{\link{Sys.timezone}} about system differences in
#' time zone names and handling.
#'
#' @param when datetime A time instant
#' @param tz.target,tz.reference character Two time zones using names
#' recognized by functions from package 'lubridate'
#'
#' @return A \code{numeric} value.
#'
#' @export
#'
tz_time_diff <- function(when = lubridate::now(),
tz.target = lubridate::tz(when),
tz.reference = "UTC") {
if (lubridate::is.Date(when)) {
when <- lubridate::as_datetime(when, tz = tz.target)
}
(as.numeric(lubridate::force_tz(when, tz.reference)) -
as.numeric(lubridate::force_tz(when, tz.target))) / 3600
}
#' Times for sun positions
#'
#' Functions for calculating the timing of solar positions, given geographical
#' coordinates and dates. They can be also used to find the time for an
#' arbitrary solar elevation between 90 and -90 degrees by supplying "twilight"
#' angle(s) as argument.
#'
#' @param date "vector" of \code{POSIXct} times or\code{Date} objects, any valid
#' TZ is allowed, default is current date at Greenwich matching the default
#' for \code{geocode}.
#' @param tz character vector indicating time zone to be used in output and to
#' interpret \code{Date} values passed as argument to \code{date}.
#' @param geocode data frame with one or more rows and variables lon and lat as
#' numeric values (degrees). If present, address will be copied to the output.
#' @param twilight character string, one of "none", "rim", "refraction",
#' "sunlight", "civil", "nautical", "astronomical", or a \code{numeric} vector
#' of length one, or two, giving solar elevation angle(s) in degrees (negative
#' if below the horizon).
#' @param unit.out character string, One of "datetime", "day", "hour", "minute",
#' or "second".
#'
#' @return A tibble with variables day, tz, twilight.rise, twilight.set,
#' longitude, latitude, address, sunrise, noon, sunset, daylength,
#' nightlength or the corresponding individual vectors.
#'
#' @family astronomy related functions
#'
#' @details Twilight names are interpreted as follows. "none": solar elevation =
#' 0 degrees. "rim": upper rim of solar disk at the horizon or solar elevation
#' = -0.53 / 2. "refraction": solar elevation = 0 degrees + refraction
#' correction. "sunlight": upper rim of solar disk corrected for refraction,
#' which is close to the value used by the online NOAA Solar Calculator.
#' "civil": -6 degrees, "naval": -12 degrees, and "astronomical": -18 degrees.
#' Unit names for output are as follows: "day", "hours", "minutes" and
#' "seconds" times for sunrise and sunset are returned as times-of-day since
#' midnight expressed in the chosen unit. "date" or "datetime" return the same
#' times as datetime objects with TZ set (this is much slower than "hours").
#' Day length and night length are returned as numeric values expressed in
#' hours when `"datetime"' is passed as argument to \code{unit.out}. If
#' twilight is a numeric vector of length two, the element with index 1 is
#' used for sunrise and that with index 2 for sunset.
#'
#' \code{is_daytime()} supports twilight specifications by name, a test
#' like \code{sun_elevation() > 0} may be used directly for a numeric angle.
#'
#' @seealso \code{\link{sun_angles}}.
#'
#' @note Function \code{day_night()} is an implementation of Meeus equations as
#' used in NOAAs on-line web calculator, which are very precise and valid for
#' a very broad range of dates. For sunrise and sunset the times are affected
#' by refraction in the atmosphere, which does in turn depend on weather
#' conditions. The effect of refraction on the apparent position of the sun is
#' only an estimate based on "typical" conditions. The more tangential to the
#' horizon is the path of the sun, the larger the effect of refraction is on
#' the times of visual occlusion of the sun behind the horizon---i.e. the
#' largest timing errors occur at high latitudes. The computation is not
#' defined for latitudes 90 and -90 degrees, i.e. at the poles.
#'
#' There exists a different R implementation of the same algorithms called
#' "AstroCalcPureR" available as function \code{astrocalc4r} in package
#' 'fishmethods'. Although the equations used are almost all the same, the
#' function signatures and which values are returned differ. In particular,
#' the implementation in 'photobiology' splits the calculation into two
#' separate functions, one returning angles at given instants in time, and a
#' separate one returning the timing of events for given dates. In
#' 'fishmethods' (= 1.11-0) there is a bug in function astrocalc4r() that
#' affects sunrise and sunset times. The times returned by the functions in
#' package 'photobiology' have been validated against the NOAA base
#' implementation.
#'
#' In the current implementation functions \code{sunrise_time},
#' \code{noon_time}, \code{sunset_time}, \code{day_length},
#' \code{night_length} and \code{is_daytime} are all wrappers
#' on \code{day_night}, so if more than one quantity is needed it is
#' preferable to directly call \code{day_night} and extract the different
#' components from the returned list.
#'
#' @section Warning: Be aware that R's \code{Date} class does not save time zone
#' metadata. This can lead to ambiguities in the current implementation
#' based on time instants. The argument passed to \code{date} should be
#' of class \code{POSIXct}, in other words an instant in time, from which
#' the correct date will be computed based on the \code{tz} argument.
#'
#' The time zone in which times passed to \code{date} as argument are
#' expressed does not need to be the local one or match the geocode, however,
#' the returned values will be in the same time zone as the input.
#'
#' @references
#' The primary source for the algorithm used is the book:
#' Meeus, J. (1998) Astronomical Algorithms, 2 ed., Willmann-Bell, Richmond,
#' VA, USA. ISBN 978-0943396613.
#'
#' A different implementation is available at
#' \url{https://github.com/NEFSC/READ-PDB-AstroCalc4R/} and in R paclage
#' 'fishmethods'. In 'fishmethods' (= 1.11-0) there is a bug in function
#' astrocalc4r() that affects sunrise and sunset times.
#'
#' An interactive web page using the same algorithms is available at
#' \url{https://gml.noaa.gov/grad/solcalc/}. There are small
#' differences in the returned times compared to our function that seem to be
#' related to the estimation of atmospheric refraction (about 0.1 degrees).
#'
#' @export
#' @examples
#' library(lubridate)
#'
#' my.geocode <- data.frame(lon = 24.93838,
#' lat = 60.16986,
#' address = "Helsinki, Finland")
#'
#' day_night(ymd("2015-05-30", tz = "EET"),
#' geocode = my.geocode)
#' day_night(ymd("2015-05-30", tz = "EET") + days(1:10),
#' geocode = my.geocode,
#' twilight = "civil")
#' sunrise_time(ymd("2015-05-30", tz = "EET"),
#' geocode = my.geocode)
#' noon_time(ymd("2015-05-30", tz = "EET"),
#' geocode = my.geocode)
#' sunset_time(ymd("2015-05-30", tz = "EET"),
#' geocode = my.geocode)
#' day_length(ymd("2015-05-30", tz = "EET"),
#' geocode = my.geocode)
#' day_length(ymd("2015-05-30", tz = "EET"),
#' geocode = my.geocode,
#' unit.out = "day")
#' is_daytime(ymd("2015-05-30", tz = "EET") + hours(c(0, 6, 12, 18, 24)),
#' geocode = my.geocode)
#' is_daytime(ymd_hms("2015-05-30 03:00:00", tz = "EET"),
#' geocode = my.geocode)
#' is_daytime(ymd_hms("2015-05-30 00:00:00", tz = "UTC"),
#' geocode = my.geocode)
#' is_daytime(ymd_hms("2015-05-30 03:00:00", tz = "EET"),
#' geocode = my.geocode,
#' twilight = "civil")
#' is_daytime(ymd_hms("2015-05-30 00:00:00", tz = "UTC"),
#' geocode = my.geocode,
#' twilight = "civil")
#'
day_night <- function(date = lubridate::now(tzone = "UTC"),
tz = ifelse(lubridate::is.Date(date),
"UTC",
lubridate::tz(date)),
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
twilight = "none",
unit.out = "hours") {
stopifnot(!anyNA(date)) # NAs could be propagated instead
tz <- unique(tz)
if (length(tz) > 1L) {
tz <- tz[1]
warning("'tz' is a heterogeneous vector, using only: ", tz)
}
geocode <- validate_geocode(geocode)
if (any(lubridate::is.Date(date))) {
date <- as.POSIXct(date, tz = tz, origin = lubridate::origin)
}
# as 'date' is not a Date but a time, we find the corresponding date in UTC
# as calculations are done in UTC time
date <- lubridate::with_tz(date, tzone = "UTC")
date <- lubridate::floor_date(date, unit = "days")
if (unit.out == "date") {
unit.out <- "datetime"
} else if (unit.out == "days") {
unit.out <- "day"
} else if (unit.out == "hours") {
unit.out <- "hour"
} else if (unit.out == "minutes") {
unit.out <- "minute"
} else if (unit.out == "seconds") {
unit.out <- "second"
}
z <- list()
for (i in seq_len(nrow(geocode))) {
temp <- day_night_fast(date = date,
tz = tz, # used for returned times
geocode = dplyr::slice(geocode, i),
twilight = twilight,
unit.out = unit.out)
z[[i]] <- temp
}
z <- do.call(rbind, z)
# we use rbind instead of dplyr::bind_rows as the second drops the class attribute
# for solartime.
# assertion
if (any(!is.na(z[["daylength"]]) & z[["daylength"]] < 0) ||
any(!is.na(z[["nightlength"]]) & z[["nightlength"]] < 0)) {
warning("Returned 'daylength/nightlength' value(s) off range")
}
attr(z, "unit.out") <- unit.out
z
}
#' @rdname day_night
#'
day_night_fast <- function(date,
tz,
geocode,
twilight,
unit.out) {
# Input validation done in day_night() before calling this function.
# stopifnot(!anyNA(time))
# date should be always a POSIXct object with tz set to "UTC" at hms all set to zero.
# stopifnot(is.data.frame(geocode))
# stopifnot(nrow(geocode == 1) && length(tz == 1))
# We have a single geocode and all dates are expressed in the same time zone!
# As date is a "vector" we can vectorize the whole calculation, and do the
# expensive calculations only once.
tz <- tz[1]
multiplier <- switch(unit.out,
hour = 1,
minute = 60,
second = 3600,
day = 1/24,
1) # default
# not vectorized, but possibly different angle for sunset and sunrise
# twilight.angles is always of length 2
twilight.angles <- twilight2angle(twilight)
# geocode passed by day_night() has always one row
if (!exists("address", geocode)) {
geocode[["address"]] <- NA_character_
}
lon <- geocode[["lon"]]
if (lon > 180 || lon < -180) {
stop("Longitude is off-range.")
}
lat <- geocode[["lat"]]
if (lat > 89.99 || lat < -89.99) {
stop("Latitude is off-range.")
}
address <- geocode[["address"]]
noon.of.date <- lubridate::with_tz(date, tzone = "UTC") + 43200 # faster
cent <- julian_century(noon.of.date)
tz.diff <- tz_time_diff(noon.of.date, tz.target = tz)
sun.lon.mean <- geom_mean_lon_sun(cent)
sun.anom.mean <- geom_mean_anom_sun(cent)
eccent.earth <- eccent_earth_orbit(cent)
delta <- sun_eq_of_ctr(cent, sun.anom.mean)
sun.lon <- sun.lon.mean + delta
# sun.anom <- sun.anom.mean + delta
# sun.dist <- sun_rad_vector(eccent.earth, sun.anom)
sun.app.lon <- sun_app_lon(cent, sun.lon)
sun.ecliptic <- mean_obliq_eclip(cent)
obliq.corr <- obliq_corr(cent, sun.ecliptic)
# rt.ascen <- sun_rt_ascen(sun.app.lon, obliq.corr)
sun.declin <- sun_decline(sun.app.lon, obliq.corr)
var.y <- var_y(obliq.corr)
eq.of.time <- eq_of_time(mean.lon = sun.lon.mean,
eccent.earth = eccent.earth,
anom.mean = sun.anom.mean,
var.y = var.y)
solar.noon <- solar_noon(lon, eq.of.time)
# We need to test for 24 h and 0 h days
sun.noon.elevation <- elevation_angle(lat, 0, sun.declin)
low.sunrise.sun <- sun.noon.elevation < max(twilight.angles[1])
low.sunset.sun <- sun.noon.elevation < max(twilight.angles[2])
sun.midnight.elevation <- elevation_angle(lat, 180, sun.declin)
high.sunrise.sun <- sun.midnight.elevation > min(twilight.angles[1])
high.sunset.sun <- sun.midnight.elevation > min(twilight.angles[2])
# Vectorized
day.start <-
ifelse(!low.sunrise.sun & !high.sunrise.sun,
# normal case
sunrise(solar.noon,
ha_sunrise(lat, sun.declin, nag = -twilight.angles[1])),
0
)
sunrise <- ifelse(day.start == 0, NA_real_, day.start)
day.end <-
ifelse(!low.sunset.sun & !high.sunset.sun,
# normal case
sunset(solar.noon,
ha_sunrise(lat, sun.declin, nag = -twilight.angles[2])),
1
)
sunset <- ifelse(day.end == 1, NA_real_, day.end)
daylength.hours <- ifelse(low.sunrise.sun & low.sunset.sun,
0,
(day.end - day.start) * 24)
# we assemble the data frame for one geographic location
# the slowest part of calculations are the calls to lubridate
# so we try to minimize them according to output format.
if (unit.out == "datetime") {
sunrise.time <- date +
lubridate::seconds(sunrise * 86400)
noon.time <- date +
lubridate::seconds(solar.noon * 86400)
sunset.time <- date +
lubridate::seconds(sunset * 86400)
tibble::tibble(day = date,
tz = rep(!!tz, length(date)),
twilight.rise = rep(twilight.angles[1], length(date)),
twilight.set = rep(twilight.angles[2], length(date)),
longitude = rep(lon, length(date)),
latitude = rep(lat, length(date)),
address = rep(address, length(date)),
sunrise = lubridate::with_tz(sunrise.time, tzone = !!tz),
noon = lubridate::with_tz(noon.time, tzone = !!tz),
sunset = lubridate::with_tz(sunset.time, tzone = !!tz),
daylength = daylength.hours,
nightlength = 24 - daylength.hours,
.name_repair = "minimal"
)
} else if (unit.out %in% c("day", "hour", "minute", "second")) {
sunrise.tod <- (sunrise * 24 + tz.diff) %% 24
noon.tod <- (solar.noon * 24 + tz.diff) %% 24
sunset.tod <- (sunset * 24 + tz.diff) %% 24
tibble::tibble(day = date,
tz = rep(!!tz, length(date)),
twilight.rise = rep(twilight.angles[1], length(date)),
twilight.set = rep(twilight.angles[2], length(date)),
longitude = rep(lon, length(date)),
latitude = rep(lat, length(date)),
address = rep(address, length(date)),
sunrise = sunrise.tod * multiplier,
noon = noon.tod * multiplier,
sunset = sunset.tod * multiplier,
daylength = daylength.hours * multiplier,
nightlength = (24 - daylength.hours) * multiplier,
.name_repair = "minimal"
)
} else {
stop("Unit out '", unit.out, "' not recognized")
}
}
#' @rdname day_night
#'
#' @return \code{is_daytime()} returns a logical vector, with \code{TRUE} for
#' day time and \code{FALSE} for night time.
#'
#' @export
#'
is_daytime <- function(date = lubridate::now(tzone = "UTC"),
tz = ifelse(lubridate::is.Date(date),
"UTC",
lubridate::tz(date)),
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
twilight = "none",
unit.out = "hours") {
if (!lubridate::is.POSIXct(date)) {
warning("'date' must be a 'POSIXct' vector")
return(rep(NA, length(date)))
}
z <- day_night(date = date,
tz = tz,
geocode = geocode,
twilight = twilight,
unit.out = "datetime")
date > z[["sunrise"]] & date < z[["sunset"]]
}
#' twilight argument check and conversion
#'
#' @return numeric Solar elevation angle at sunrise or sunset
#' @keywords internal
#'
twilight2angle <- function(twilight) {
# refraction depends on the elevation angle (ha), atmospheric pressure and temperature
# we assume atmospheric pressure 1010 hPa and temperature 10 C
# constants below are approximate
# function atm_refraction_approx() can be used to estimate refraction for other angles
if (length(twilight) == 0) {
warning("'twilight' too short or NULL")
twilight_angle <- c(NA_real_, NA_real_)
} else if (is.character(twilight)) {
if (length(twilight) != 1) {
warning("'twilight' is character but of length > 1, using first value")
twilight <- twilight[1L]
}
twilight_angle <-
switch(twilight,
none = c(0, 0), # center of solar disk
rim = c(-0.53/2, -0.53/2), # upper rim of solar disk
refraction = c(-0.4819444, -0.4819444), # center of solar disk, refraction corrected
sunlight = c(-0.833, -0.833),
civil = c(-6, -6),
nautical = c(-12, -12),
astronomical = c(-18, -18),
{warning("Unrecognised 'twilight' value: ", twilight); c(NA_real_, NA_real_)}
)
} else if (is.numeric(twilight)) {
if (length(twilight) == 1) {
twilight_angle <- rep(twilight, 2)
} else if (length(twilight) == 2) {
twilight_angle <- twilight
} else {
stop("'twilight' vector longer than 2")
}
if (any(!is.na(twilight_angle) &
twilight_angle < -90 | twilight_angle > 90 ) ) {
warning("Off range twilight angle(s) replaced with NAs.")
twilight_angle <-
ifelse(!is.na(twilight_angle) &
twilight_angle > 90 | twilight_angle < -90,
NA_real_, twilight_angle)
}
} else {
warning("'twilight' must be numeric or character, but is ", class(twilight))
twilight_angle <- c(NA_real_, NA_real_)
}
twilight_angle
}
#' @rdname day_night
#' @export
#' @return \code{noon_time}, \code{sunrise_time} and \code{sunset_time} return a
#' vector of POSIXct times
#'
noon_time <- function(date = lubridate::now(tzone = "UTC"),
tz = lubridate::tz(date),
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
twilight = "none",
unit.out = "datetime") {
stopifnot(length(date) == 1 || nrow(geocode) == 1)
day_night(date = date,
tz = tz,
geocode = geocode,
twilight = twilight,
unit.out = unit.out)[["noon"]]
}
#' @rdname day_night
#'
#' @export
#'
sunrise_time <- function(date = lubridate::now(tzone = "UTC"),
tz = lubridate::tz(date),
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
twilight = "sunlight",
unit.out = "datetime") {
# stopifnot(length(date) == 1L || nrow(geocode) == 1L)
day_night(date = date,
tz = tz,
geocode = geocode,
twilight = twilight,
unit.out = unit.out)[["sunrise"]]
}
#' @rdname day_night
#'
#' @export
#'
sunset_time <- function(date = lubridate::now(tzone = "UTC"),
tz = lubridate::tz(date),
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
twilight = "sunlight",
unit.out = "datetime") {
# stopifnot(length(date) == 1L || nrow(geocode) == 1L)
day_night(date = date,
tz = tz,
geocode = geocode,
twilight = twilight,
unit.out = unit.out)[["sunset"]]
}
#' @rdname day_night
#'
#' @export
#' @return \code{day_length} and \code{night_length} return numeric a vector
#' giving the length in hours
day_length <- function(date = lubridate::now(tzone = "UTC"),
tz = "UTC",
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
twilight = "sunlight", unit.out = "hours") {
# stopifnot(length(date) == 1L || nrow(geocode) == 1L)
day_night(date = date,
tz = tz,
geocode = geocode,
twilight = twilight,
unit.out = unit.out)[["daylength"]]
}
#' @rdname day_night
#'
#' @export
#' @note \code{night_length} returns the length of night-time conditions in one
#' day (00:00:00 to 23:59:59), rather than the length of the night between two
#' consecutive days.
night_length <- function(date = lubridate::now(tzone = "UTC"),
tz = "UTC",
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
twilight = "sunlight", unit.out = "hours") {
# stopifnot(length(date) == 1L || nrow(geocode) == 1L)
day_night(date = date,
tz = tz,
geocode = geocode,
twilight = twilight,
unit.out = unit.out)[["nightlength"]]
}
#' Convert datetime to time-of-day
#'
#' Convert a datetime into a time of day expressed in hours, minutes or seconds
#' from midnight in local time for a time zone. This conversion is useful when
#' time-series data for different days needs to be compared or plotted based on
#' the local time-of-day.
#'
#' @param x a datetime object accepted by lubridate functions
#' @param unit.out character string, One of "tod_time", "hours", "minutes", or "seconds".
#' @param tz character string indicating time zone to be used in output.
#'
#' @return A numeric vector of the same length as \code{x}. If
#' \code{unit.out = "tod_time"} an object of class \code{"tod_time"} which
#' the same as for \code{unit.out = "hours"} but with the class attribute
#' set, which dispatches to special \code{format()} nad \code{print()}
#' methods.
#'
#' @seealso \code{\link{solar_time}}
#'
#' @family Time of day functions
#'
#' @export
#'
#'
#' @examples
#' library(lubridate)
#' my_instants <- ymd_hms("2020-05-17 12:05:03") + days(c(0, 30))
#' my_instants
#' as_tod(my_instants)
#' as_tod(my_instants, unit.out = "tod_time")
#'
as_tod <- function(x, unit.out = "hours", tz = NULL) {
stopifnot(lubridate::is.timepoint(x))
if (!is.null(tz)) {
x <- lubridate::with_tz(x, tzone = tz[1])
}
if (unit.out == "tod_time") {
tod <- lubridate::hour(x) + lubridate::minute(x) / 60 + lubridate::second(x) / 3600
class(tod) <- c("tod_time", class(tod))
tod
} else if (unit.out == "hours") {
lubridate::hour(x) + lubridate::minute(x) / 60 + lubridate::second(x) / 3600
} else if (unit.out == "minutes") {
lubridate::hour(x) * 60 + lubridate::minute(x) + lubridate::second(x) / 60
} else if (unit.out == "seconds") {
lubridate::hour(x) * 3600 + lubridate::minute(x) * 60 + lubridate::second(x)
} else {
stop("Unrecognized 'unit.out': ", unit.out)
}
}
#' Encode in a Common Format
#'
#' Format a \code{tod_time} object for pretty printing
#'
#' @param x an R object
#' @param ... ignored
#' @param sep character used as separator
#'
#' @family Time of day functions
#'
#' @export
#'
format.tod_time <- function(x, ..., sep = ":") {
hours <- as.integer(trunc(x))
minutes <- as.integer((x * 60) %% 60)
seconds <- as.integer((x * 3600) %% 60)
fmt <- paste(rep("%02d", 3), collapse = sep)
time_string <-
sprintf(fmt = fmt, hours, minutes, seconds)
time_string
}
#' Print time-of-day objects
#'
#' @param x an R object
#' @param ... passed to \code{format} method
#'
#' @family Time of day functions
#'
#' @note Default is to print the underlying \code{numeric} vector as a solar time.
#'
#' @export
#'
print.tod_time <- function(x, ...) {
print(format(x, ...))
invisible(x)
}
#' Local solar time
#'
#' \code{solar_time()} computes the time of day expressed in seconds since the
#' astronomical midnight using and instant in time and a geocode as input. Solar
#' time is useful when we want to plot data according to the local solar time
#' rather than the local time in use at a time zone. How the returned instant in
#' time is expressed depends on the argument passed to \code{unit.out}.
#'
#' @details Solar time is determined by the position of the sun in the sky and
#' it almost always differs from the time expressed in the local time
#' coordinates in use. The differences can vary from a few minutes up to a
#' couple of hours depending on the exact location within the time zone and the
#' use or not of daylight saving time.
#'
#' @param time POSIXct Time, any valid time zone (TZ) is allowed, default is
#' current time
#' @param geocode data frame with variables lon and lat as numeric values
#' (degrees).
#' @param unit.out character string, One of "datetime", "time", "hour", "minute", or
#' "second".
#'
#' @seealso \code{\link{as_tod}}
#'
#' @family Local solar time functions
#'
#' @note The algorithm is approximate, it calculates the difference between
#' local solar noon and noon in the time zone of \code{time} and uses this
#' value for the whole day when converting times into solar time. Days are not
#' exactly 24 h long. Between successive days the shift is only a few seconds,
#' and this leads to a small jump at midnight.
#'
#' @section Warning!: Returned values are computed based on the time zone of the
#' argument for parameter time. In the case of solar time, this timezone does
#' not affect the result. However, in the case of solar dates the date part
#' may be off by one day, if the time zone does not match the coordinates of
#' the geocode value provided as argument.
#'
#' @return In all cases solar time is expressed as time since local astronomical
#' midnight and, thus, lacks date information. If \code{unit.out = "time"}, a
#' numeric value in seconds with an additional class attribute
#' "solar_time"; if \code{unit.out = "datetime"}, a "POSIXct" value in seconds
#' from midnight but with an additional class attribute "solar_date"; if
#' \code{unit.out = "hour"} or \code{unit.out = "minute"} or \code{unit.out =
#' "second"}, a numeric value.
#'
#' @export
#'
#' @examples
#' BA.geocode <-
#' data.frame(lon = -58.38156, lat = -34.60368, address = "Buenos Aires, Argentina")
#' sol_t <- solar_time(lubridate::dmy_hms("21/06/2016 10:00:00", tz = "UTC"),
#' BA.geocode)
#' sol_t
#' class(sol_t)
#'
#' sol_d <- solar_time(lubridate::dmy_hms("21/06/2016 10:00:00", tz = "UTC"),
#' BA.geocode,
#' unit.out = "datetime")
#' sol_d
#' class(sol_d)
#'
solar_time <- function(time = lubridate::now(),
geocode = tibble::tibble(lon = 0,
lat = 51.5,
address = "Greenwich"),
unit.out = "time")
{
# solar time in hours from midnight
solar.time <- sun_angles(time = time,
tz = lubridate::tz(time),
geocode = geocode)[["solartime"]]
switch(unit.out,
"date" = as.solar_date(solar.time, time),
"datetime" = as.solar_date(solar.time, time),
"time" = solar.time,
"days" = as.numeric(solar.time) / 24,
"hours" = as.numeric(solar.time),
"minutes" = as.numeric(solar.time * 60),
"seconds" = as.numeric(solar.time * 3600)
)
}
#' Convert a solar_time object into solar_date object
#'
#' @param x solar_time object.
#' @param time an R date time object
#'
#' @return For method \code{as.solar_date()} a date-time object with the class attr
#' set to "solar.time". This is needed only for unambiguous formatting and
#' printing.
#'
#' @family Local solar time functions
#'
#' @export
#'
as.solar_date <- function(x, time)
{
stopifnot(is.solar_time(x))
stopifnot(lubridate::is.timepoint(time))
solar.date <-
lubridate::floor_date(time, unit = "days") +
lubridate::seconds(as.numeric(x) * 3600)
class(solar.date) <- c("solar_date", class(solar.date))
solar.date
}
#' Query class
#'
#' @param x an R object.
#'
#' @family Local solar time functions
#'
#' @export
#'
is.solar_time <- function(x) {
inherits(x, "solar_time")
}
#' @rdname is.solar_time
#'
#' @export
#'
is.solar_date <- function(x) {
inherits(x, "solar_date")
}
#' Encode in a Common Format
#'
#' Format a \code{solar_time} object for pretty printing
#'
#' @param x an R object
#' @param ... ignored
#' @param sep character used as separator
#'
#' @family astronomy related functions
#'
#' @export
#'
format.solar_time <- function(x, ..., sep = ":") {
hours <- as.integer(trunc(x))
minutes <- as.integer((x * 60) %% 60)
seconds <- as.integer((x * 3600) %% 60)
fmt <- paste(rep("%02d", 3), collapse = sep)
time_string <-
sprintf(fmt = fmt, hours, minutes, seconds)
time_string
}
#' Print solar time and solar date objects
#'
#' @param x an R object
#' @param ... passed to \code{format} method
#'
#' @family Local solar time functions
#'
#' @note Default is to print the underlying POSIXct as a solar time.
#'
#' @export
#'
print.solar_time <- function(x, ...) {
print(format(x, ...))
invisible(x)
}
#' @rdname print.solar_time
#'
#' @export
#'
print.solar_date <- function(x, ...) {
print(paste(format(x, ...), "solar"))
invisible(x)
}
#' Extraterrestrial irradiance
#'
#' Estimate of down-welling solar (short wave) irradiance at the top of the
#' atmosphere above a location on Earth, computed based on angles, Sun-Earth
#' distance and the solar constant. Astronomical computations are done with
#' function \code{sun_angles()}.
#'
#' @param time A "vector" of POSIXct Time, with any valid time zone (TZ) is
#' allowed, default is current time.
#' @param tz character string indicating time zone to be used in output.
#' @param geocode data frame with variables lon and lat as numeric values
#' (degrees), nrow > 1, allowed.
#' @param solar.constant numeric or character If character, "WMO" or "NASA", if
#' numeric, an irradiance value in the same units as the value to be returned.
#'
#' @return Numeric vector of extraterrestrial irradiance (in W / m2 if solar
#' constant is a character value).
#'
#' @seealso Function \code{\link{sun_angles}}.
#'
#' @examples
#' library(lubridate)
#'
#' irrad_extraterrestrial(ymd_hm("2021-06-21 12:00", tz = "UTC"))
#'
#' irrad_extraterrestrial(ymd_hm("2021-12-21 20:00", tz = "UTC"))
#'
#' irrad_extraterrestrial(ymd_hm("2021-06-21 00:00", tz = "UTC") + hours(1:23))
#'
#' @export
#'
irrad_extraterrestrial <-
function(time = lubridate::now(tzone = "UTC"),
tz = lubridate::tz(time),
geocode = tibble::tibble(lon = 0, lat = 51.5, address = "Greenwich"),
solar.constant = "NASA") {
solar.cnst.map <- c(NASA = 1360, WMO = 1367) # W / m2
if (is.character(solar.constant)) {
solar.constant <- solar.cnst.map[solar.constant]
}
angles <- sun_angles(time = time,
tz = tz,
geocode = geocode,
use.refraction = FALSE) # no atmosphere!
rel.distance <- angles[["distance"]]
sun.elevation <- angles[["elevation"]]
ifelse(sun.elevation <= 0,
0,
solar.constant *
cos((90 - sun.elevation) / 180 * pi) / # degrees -> radians
rel.distance^2)
}
#' Validate a geocode
#'
#' Test validity of a geocode or ensure that a geocode is valid.
#'
#' @details
#' \code{validate_geocode} Converts to tibble, checks data bounds, converts
#' address to character if it is not already a character vector, or add
#' character NAs if the address column is missing.
#'
#' \code{is_valid_geocode} Checks if a geocode is valid, returning 0L if not,
#' and the number of row otherwise.
#'
#' @param geocode data.frame with geocode data in columns \code{"lat"},
#' \code{"lon"}, and possibly also \code{"address"}.
#'
#' @return A valid geocode stored in a tibble.
#'
#' @examples
#'
#' validate_geocode(NA)
#' validate_geocode(data.frame(lon = -25, lat = 66))
#'
#' is_valid_geocode(NA)
#' is_valid_geocode(1L)
#' is_valid_geocode(data.frame(lon = -25, lat = 66))
#'
#' na_geocode()
#'
#' @export
#'
validate_geocode <- function(geocode) {
if (is.atomic(geocode) && (length(geocode) == 1L) && is.na(geocode)) {
geocode <- na_geocode()
} else if (is.data.frame(geocode)) {
geocode <- tibble::as_tibble(geocode, .name_repair = "minimal")
} else {
stop("Bad geocode: ", format(geocode))
}
stopifnot(nrow(geocode) >= 1) # needs to be replace by generation of no output in all functions
stopifnot(exists("lon", geocode), exists("lat", geocode))
geocode[["lon"]] <- as.numeric(geocode[["lon"]]) # convert logical NA
geocode[["lat"]] <- as.numeric(geocode[["lat"]]) # convert logical NA
if (any(na.omit(geocode[["lon"]]) > 180 | na.omit(geocode[["lon"]]) < -180)) {
stop("Longitude is off-range.")
}
if (any(na.omit(geocode[["lat"]]) > 89.99 | na.omit(geocode[["lat"]]) < -89.99)) {
stop("Latitude is off-range.")
}
if (!exists("address", geocode)) {
geocode[["address"]] <- NA_character_
} else if (!is.character(geocode[["address"]])) {
geocode[["address"]] <- as.character(geocode[["address"]])
}
geocode
}
#' @rdname validate_geocode
#'
#' @return FALSE for invalid, TRUE for valid.
#'
#' @export
#'
is_valid_geocode <- function(geocode) {
if (!is.list(geocode)) return(FALSE)
if (!is.data.frame(geocode)) {
# walk list of geocodes using recursion
is_valid <- all(sapply(geocode, is_valid_geocode))
} else {
is_valid <- nrow(geocode) >= 1L &&
all(c("lon", "lat") %in% names(geocode)) &&
all(c(is.numeric(geocode[["lon"]]), is.numeric(geocode[["lat"]]))) &&
if ("address" %in% names(geocode)) is.character(geocode[["address"]]) else TRUE
if (!is_valid && "address" %in% names(geocode) && is.factor(geocode[["address"]])) {
warning("'address' is a factor instead of a character vector.")
}
}
is_valid
}
#' @rdname validate_geocode
#'
#' @return FALSE for invalid, number of rows for valid.
#'
#' @export
#'
length_geocode <- function(geocode) {
if (!is_valid_geocode(geocode)) return(NA_integer_)
nrow(geocode)
}
#' @rdname validate_geocode
#'
#' @return A geo_code tibble with all fields set to suitable NAs.
#'
#' @export
#'
na_geocode <- function() {
tibble::tibble(lon = NA_real_,
lat = NA_real_,
address = NA_character_)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.