Merge pull request #339 from adokter/develop

new options in sunrise(), sunset() and check_night()

NEWS.md

@@ -2,7 +2,11 @@
 
 * `regularize_vpts()` is now much faster, and chooses more intuitive starting and ending point of the regularized grid, e.g. projecting on half hour grid will have time series start on the nearest half hour (#332)
 
-* new option `keep_timestamp` in `regularize_vpts()`, which allows individual profiles to keep there original timestamp instead of the timestamp of the regularized grid
+* new option `keep_timestamp` in `regularize_vpts()`, which allows individual profiles to keep there original timestamp instead of the timestamp of the regularized grid.
+
+* improved documentation of sunrise/sunset functions (#180) and new option `force_tz` (4968019).
+
+* new option `offset` in `check_night()`, which allows day/night transition to be shifted by a temporal offset (#338). For example, this is useful when selecting night time profiles that start a specific number of hours after sunset.
 
 # bioRad 0.5.1
 

R/check_night.R

@@ -11,7 +11,11 @@
 #' @param lon numeric. Longitude in decimal degrees.
 #' @param lat numeric. Latitude in decimal degrees.
 #' @param tz character. Time zone. Ignored when \code{date} already has an associated time zone
-#' @param elev numeric. Sun elevation in degrees.
+#' @param elev numeric. Sun elevation in degrees defining night time. May also be a numeric vector of
+#' length two, with first element giving sunset elevation, and second element sunrise elevation.
+#' @param offset numeric. Time duration in seconds by which the shift the start and end
+#' of night time. May also be a numeric vector of length two, with first element added to moment
+#' of sunset and second element added to moment of sunrise.
 #' @param ... optional lat,lon arguments.
 #'
 #' @return \code{TRUE} when night, \code{FALSE} when day, \code{NA} if unknown
@@ -30,6 +34,10 @@
 #' Approximate astronomical formula are used, therefore the day/night
 #' transition may be off by a few minutes.
 #'
+#' \code{offset} can be used to shift the moment of sunset and sunrise by a
+#' temporal offset, for example, \code{offset=c(600,-600)} will assume nighttime
+#' starts 600 seconds after sunset (as defined by \code{elev}) and stops 600 seconds before sunrise.
+#'
 #' @examples
 #' # check if it is night at UTC midnight in the Netherlands on January 1st:
 #' check_night("2016-01-01 00:00", 5, 53)
@@ -38,80 +46,113 @@
 #' check_night(example_vp)
 #'
 #' check_night(example_vpts)
-check_night <- function(x, ..., elev = -0.268) {
+#'
+#' # select nighttime profiles that are between 3 hours after sunset
+#' # and 2 hours before sunrise:
+#' index <- check_night(example_vpts, offset=c(3,-2)*3600)
+#' example_vpts[index]
+check_night <- function(x, ..., elev = -0.268, offset = 0) {
   UseMethod("check_night", x)
 }
 
 #' @rdname check_night
 #'
 #' @export
-check_night.default <- function(x, lon, lat, ..., tz = "UTC", elev = -0.268) {
+check_night.default <- function(x, lon, lat, ..., tz = "UTC", elev = -0.268, offset = 0) {
+  # input checks
+  assert_that(is.numeric(elev))
+  assert_that(length(elev)<=2)
+  assert_that(is.numeric(offset))
+  assert_that(length(offset)<=2)
+  assert_that(is.numeric(lon))
+  assert_that(is.numeric(lat))
+  #
   x <- as.POSIXct(x, tz = tz)
+  #
+  elev_sunset = ifelse(length(elev) == 2,elev[1],elev)
+  elev_sunrise = ifelse(length(elev) == 2,elev[2],elev)
   # calculate sunrises
-  trise <- sunrise(x, lon, lat, tz = tz, elev = elev)
+  trise <- sunrise(x, lon, lat, tz = tz, elev = elev_sunrise)
   # calculate sunsets
-  tset <- sunset(x, lon, lat, tz = tz, elev = elev)
-  # for returned rise times on a different day, recalculate for the current day
-  dt <- as.numeric(difftime(as.Date(trise), as.Date(x), units = "days"))
-  change <- which(dt != 0)
-  if (length(change) > 0) {
-    trise[change] <- sunrise(x[change] - dt[change] * 24 * 3600, lon, lat,
-      tz = tz, elev = elev
-    )
-  }
+  tset <- sunset(x, lon, lat, tz = tz, elev = elev_sunset)
 
-  dt <- as.numeric(difftime(as.Date(tset), as.Date(x), units = "days"))
-  change <- which(dt != 0)
-  if (length(change) > 0) {
-    tset[change] <- sunset(x[change] - dt[change] * 24 * 3600, lon, lat,
-      tz = tz, elev = elev
-    )
-  }
+  # make sure the observation is always in between
+  # two subsequent sunrise/sunset events:
+  trise_ordered = trise
+  tset_ordered = tset
+  # sunrise -1 day
+  idx <- which(x < tset  & trise > tset)
+  if(length(idx)>0) trise_ordered[idx] = sunrise(x[idx]-24*3600, lon, lat, tz = tz, elev = elev_sunrise)
+  # sunrise +1 day
+  idx <- which(x > tset  & trise < tset)
+  if(length(idx)>0) trise_ordered[idx] = sunrise(x[idx]+24*3600, lon, lat, tz = tz, elev = elev_sunrise)
+  # sunset -1 day
+  idx <- which(x < trise & trise < tset)
+  if(length(idx)>0) tset_ordered[idx] = sunset(x[idx]-24*3600, lon, lat, tz = tz, elev = elev_sunset)
+  # sunset +1 day
+  idx <- which(x > trise & trise > tset)
+  if(length(idx)>0) tset_ordered[idx] = sunset(x[idx]+24*3600, lon, lat, tz = tz, elev = elev_sunset)
+  # store in original variable
+  trise <- trise_ordered
+  tset <- tset_ordered
+
+  # add offset shifts
+  offset_sunset = ifelse(length(offset) == 2,offset[1],offset)
+  offset_sunrise = ifelse(length(offset) == 2,offset[2],offset)
 
   # prepare output
   output <- rep(NA, length(x))
-  itsday <- (x > trise & x < tset)
+
+  # cases trise < tset
+  itsday <- (x > trise + offset_sunrise & x < tset + offset_sunset)
   output[trise < tset] <- itsday[trise < tset]
-  itsday <- (x < tset | x > trise)
+  # if order of sunrise and sunset switches due to offsets, daytime length is zero
+  output[trise < tset  & (tset + offset_sunset) - (trise + offset_sunrise) <= 0] <- FALSE
+
+  # cases trise >= tset
+  itsday <- (x < tset + offset_sunset | x > trise + offset_sunrise)
   output[trise >= tset] <- itsday[trise >= tset]
+  # if order of sunrise and sunset switches due to offsets, nighttime length is zero
+  output[trise >= tset & (trise + offset_sunrise) - (tset + offset_sunset) <= 0] <- TRUE
+
   !output
 }
 
 #' @rdname check_night
 #'
 #' @export
-check_night.vp <- function(x, ..., elev = -0.268) {
+check_night.vp <- function(x, ..., elev = -0.268, offset = 0) {
   stopifnot(inherits(x, "vp"))
   check_night(x$datetime, x$attributes$where$lon, x$attributes$where$lat,
-    elev = elev
+    elev = elev, offset = offset
   )
 }
 
 #' @rdname check_night
 #'
 #' @export
-check_night.list <- function(x, ..., elev = -0.268) {
+check_night.list <- function(x, ..., elev = -0.268, offset = 0) {
   vptest <- sapply(x, function(y) is(y, "vp"))
   if (FALSE %in% vptest) {
     stop("requires list of vp objects as input")
   }
-  sapply(x, check_night.vp, elev = elev)
+  sapply(x, check_night.vp, elev = elev, offset = offset)
 }
 
 #' @rdname check_night
 #'
 #' @export
-check_night.vpts <- function(x, ..., elev = -0.268) {
+check_night.vpts <- function(x, ..., elev = -0.268, offset = 0) {
   stopifnot(inherits(x, "vpts"))
   check_night(x$datetime, x$attributes$where$lon, x$attributes$where$lat,
-    elev = elev
+    elev = elev, offset = offset
   )
 }
 
 #' @rdname check_night
 #'
 #' @export
-check_night.pvol <- function(x, ..., elev = -0.268) {
+check_night.pvol <- function(x, ..., elev = -0.268, offset = 0) {
   stopifnot(inherits(x, "pvol"))
-  check_night(x$datetime, x$geo$lon, x$geo$lat, elev = elev)
+  check_night(x$datetime, x$geo$lon, x$geo$lat, elev = elev, offset = offset)
 }

R/sunrise_sunset.R

@@ -8,10 +8,20 @@
 #' @param lat Latitude in decimal degrees.
 #' @param elev Sun elevation in degrees.
 #' @param tz output time zone. Ignored if \code{date} has an associated time zone already
+#' @param force_tz whether to convert output to timezone \code{tz}. Default \code{FALSE}.
 #'
-#' @return The moment of sunrise or sunset in UTC time.
+#' @return The moment of sunrise or sunset for the date set by \code{date}and time zone as specified
+#' (by \code{date} and \code{tz}) or in UTC if not specified.
 #'
-#' @details The angular diameter of the sun is about 0.536 degrees,
+#' @details
+#' The day for which sunrise and sunset are calcuated is given by the input date.
+#' Sunrise and sunset are calculated relative to the moment of solar noon for that date,
+#' i.e. the first sunrise before the moment of solar noon, and the first sunset after the
+#' moment of solar noon. Therefore, depending on the timezone provided, it is possible that
+#' the nearest sunrise prior to solar noon occurs a day earlier than the input date.
+#' Similarly, sunset may occur a day later than the input date. See examples for details.
+#'
+#' The angular diameter of the sun is about 0.536 degrees,
 #' therefore the moment of sunrise/sunset corresponds to half that elevation
 #' at -0.268 degrees.
 #'
@@ -20,36 +30,65 @@
 #' Approximate astronomical formula are used, therefore the moment of
 #' sunrise / sunset may be off by a few minutes
 #'
+#' If \code{force_tz} is \code{TRUE}, the output is converted to the timezone
+#' set by \code{tz}
+#'
 #' @examples
 #' # sunrise in the Netherlands
 #' sunrise("2016-01-01", 5, 53)
 #'
 #' # sunset in the Netherlands
 #' sunset("2016-01-01", 5, 53)
 #'
-#' # civil twilight in Ithaca, NY, today
-#' sunrise(Sys.time(), -76.5, 42.4, elev = -6)
+#' # civil twilight in Ithaca, NY
+#' sunrise("2016-01-01", -76.5, 42.4, elev = -6)
+#'
+#' # next sunset in South Dakota, USA
+#' sunset("2016-11-15", -98, 45)
+#'
+#' # Beware that some days have two sunsets, or
+#' # two sunrises! E.g. on 5 Oct (local timezone) at
+#' # this location  sunset is actually on the 6 Oct
+#' # in UTC time zone, i.e. the next day
+#' sunset("2016-10-5", -98, 45)
+#' # One day later, sunset is again on 6 Oct:
+#' sunset("2016-10-6", -98, 45)
+#'
+#' # working in local time zones typically avoids such ambiguities:
+#' sunset(as_datetime("2016-06-05",tz="America/Chicago"), -98, 45)
+#' sunset(as_datetime("2016-06-06",tz="America/Chicago"), -98, 45)
+#'
+#' # use force_tz to force output to a specific time zone, by default UTC:
+#' sunset(as_datetime("2016-06-05",tz="America/Chicago"), -98, 45, force_tz=T)
+#' sunset(as_datetime("2016-06-06",tz="America/Chicago"), -98, 45, force_tz=T)
+#'
+#' # Also beware of jumps in sunrise and sunset date with longitude:
+#' sunrise("2016-11-01", 100, 45)
+#' sunrise("2016-11-01", 102, 45)
+#'
 #' @name sunrise_sunset
 NULL
 
 #' @rdname sunrise_sunset
 #'
 #' @export
-sunrise <- function(date, lon, lat, elev = -0.268, tz = "UTC") {
+sunrise <- function(date, lon, lat, elev = -0.268, tz = "UTC", force_tz = FALSE) {
   locations <- data.frame(lon = lon, lat = lat)
   locations <- SpatialPoints(locations, proj4string = CRS("+proj=longlat +datum=WGS84"))
-  datetime <- as.POSIXct(date, tz = tz)
+  datetime <- as.POSIXct(date, tz = tz) # tz ignored if already set
   suntimes <- crepuscule(locations, datetime, solarDep = -elev, direction = "dawn", POSIXct.out = TRUE)
+  if(force_tz) suntimes$time <- as_datetime(suntimes$time, tz=tz)
   suntimes$time
 }
 
 #' @rdname sunrise_sunset
 #'
 #' @export
-sunset <- function(date, lon, lat, elev = -0.268, tz = "UTC") {
+sunset <- function(date, lon, lat, elev = -0.268, tz = "UTC", force_tz = FALSE) {
   locations <- data.frame(lon = lon, lat = lat)
   locations <- SpatialPoints(locations, proj4string = CRS("+proj=longlat +datum=WGS84"))
-  datetime <- as.POSIXct(date, tz = tz)
+  datetime <- as.POSIXct(date, tz = tz) # tz ignored if already set
   suntimes <- crepuscule(locations, datetime, solarDep = -elev, direction = "dusk", POSIXct.out = TRUE)
+  if(force_tz) suntimes$time <- as_datetime(suntimes$time, tz=tz)
   suntimes$time
 }