R/noaa.sunrise.code.R
In karlropkins/grey.area: This is a bit of a grey area
Defines functions
hour.offset
cutDaylight
calcNOAASunrise
Documented in
calcNOAASunrise
cutDaylight
hour.offset
####################################################
#' @title NOAA Sunrise Calculations
####################################################
#'
#' @name noaa.sunrise.code
#' @aliases noaa.sunrise calcNOAASunrise cutDaylight
#' hour.offset
#' @description Methods for the calculation of sunrise
#' and sunset times based on NOAA methods.
#' @param time.stamp (POSIX) A vector of date/time of
#' \code{POSIX*} class.
#' @param latitude,longitude (Numerics) The latitude and
#' longitude of the location where the associated \code{time.stamp}
#' measurement was logged. (Note: \code{latitude} + to North; - to
#' South. \code{longitude} + to East; - to West.)
#' @param local.hour.offset (Numeric or character) The associated
#' local time offset relative to GMT in hours (numeric) or a reference
#' term that \code{\link{hour.offset}} can convert into a numeric time
#' offset (character). (Note: + to East; - to West).
#' @param output (Character) Either the default \code{"all"}, which
#' returns all inputs and calculated parameters as a data.frame, or a
#' character vector of field/column names specifying which fields/columns
#' to return. See Value below for summary of calculated parameters
#' generated by this function.
#' @param mydata (Data.frame). A date frame which includes the
#' \code{time.stamp} vector as a field/column called \code{date}.
#' @param local (Character) The location (or reference) for required
#' hour offset.
#' @param pmatch,ignore.case (Logicals) Supplied \code{local} handling.
#' \code{pmatch} applies partial (leading edge) matching of reference
#' \code{local} terms. \code{ignore.case} resets \code{local} to upper
#' case before comparing with reference \code{local} terms. Both defaults
#' \code{TRUE}.
#' @details \code{calcNOAASunrise} is based on NOAA methods and calculates
#' a range of parameters based on the relative position of the Sun and a
#' receptor point of the Earth`s surface (see value below). The method is
#' valid for dates between 1901 and 2099.\code{cutDaylight} is a cut-down
#' version of the same method that was used in the \code{openair} package.
#'
#' \code{hour.offset} uses \code{\link{ref.hour.offset}} as a look-up
#' table to match supplied \code{local} terms and return
#' \code{hour.offset} values for use in \code{\link{calcNOAASunrise}}
#' or \code{\link{cutDaylight}}.
#' @return \code{calcNOAASunrise} returns a data frame which, depending
#' on the call \code{output} value, may include the call arguments and
#' any of the following calculated parameters:
#' \describe{
#' \item{julian.day }{(Numeric) The Julian day.}
#' \item{julian.century }{(Numeric) The Julian century.}
#' \item{geom.mean.long.sun.deg }{(Numeric) The geometric mean longitude
#' of the Sun in degrees.}
#' \item{geom.mean.anom.sun.deg }{(Numeric) The geometric mean anomaly of
#' the Sun in degrees.}
#' \item{eccent.earth.orbit }{(Numeric) The eccentricity of Earth`s orbit.}
#' \item{sun.eq.of.ctr }{(Numeric) The estimated center for the Sun, in
#' degrees.}
#' \item{sun.true.long.deg }{(Numeric) The true longitude of the Sun, in
#' degrees.}
#' \item{sun.true.anom.deg }{(Numeric) The true anamoly of the Sun, in
#' degrees.}
#' \item{sun.rad.vector.aus }{(Numeric) The distance to the Sun, in AUs
#' (radial).}
#' \item{sun.app.long.deg }{(Numeric) The apparent longitude of the Sun,
#' in degrees.}
#' \item{mean.obliq.ecliptic.deg }{(Numeric) The mean obliquity of the
#' ecliptic, in degrees.}
#' \item{obliq.corr.deg }{(Numeric) The corrected obliquity of the
#' ecliptic, in degrees.}
#' \item{sun.rt.ascen.deg }{(Numeric) The right ascension of the Sun,
#' in degress.}
#' \item{sun.declin.deg }{(Numeric) The declination of the Sun, in
#' degrees.}
#' \item{vary }{(Numeric) TO BE CONFIRMED}
#' \item{eq.of.time.minutes }{(Numeric) TO BE CONFIRMED}
#' \item{ha.sunrise.deg }{(Numeric) TO BE CONFIRMED}
#' \item{solar.noon.lst }{(Numeric) TO BE CONFIRMED}
#' \item{sunrise.time.lst }{(Numeric) The estimated sunrise time as
#' fraction of day.}
#' \item{sunset.time.lst }{(Numeric) The estimated sunset time as
#' fraction of day.}
#' \item{sunlight.duration.minutes }{(Numeric) TO BE CONFIRMED}
#' \item{true.solar.time.min }{(Numeric) TO BE CONFIRMED}
#' \item{hour.angle.deg }{(Numeric) TO BE CONFIRMED}
#' \item{solar.zenith.angle.deg }{(Numeric) TO BE CONFIRMED}
#' \item{solar.elevation.angle.deg }{(Numeric) TO BE CONFIRMED}
#' \item{approx.atmospheric.refraction.deg }{(Numeric) TO BE CONFIRMED}
#' \item{solar.elevation.corrected.for.atm.refraction.deg }{(Numeric) TO
#' BE CONFIRMED}
#' \item{solar.azimuth.angle.deg.cw.from.n }{(Numeric) TO BE CONFIRMED}
#' \item{daylight }{(Factor) A factor specifying if the \code{time.stamp},
#' \code{latitude}, \code{longitude} and \code{local.hour.offset} is
#' estimated to be a \code{daylight} or \code{nighttime} instance.}
#' }
#' \code{cutDaylight} is an alternative to \code{calcNOAASunrise} which
#' accepts the time stamp as part of a data frame and returns that data
#' frame with the \code{calcNOAACalc} output component \code{daylight}
#' as an addition field/column.
#'
#' \code{hour.offset} returns the hour.offset as a numeric.
#' @references Methods based on:
#' \url{http://www.srrb.noaa.gov/highlights/sunrise/sunrise.html}
#' @author Karl Ropkins
#' @note The NOAA Sunrise methods were originally recoded for R as part
#' of work on the NERC penair project:
#' \url{http://www.openair-project.org/}.
#' \code{cutDaylight} is an alternative version of the function. It is
#' highly similar to code used as part of the function \code{cutData} in
#' the package \code{openair} to provide \code{daylight} conditioning for
#' many plot types within that package.
#'
#' Character handling for \code{local.hour.offset} is by
#' \code{hour.offset} and \code{\link{ref.hour.offset}}.
#' @seealso \code{\link{ref.hour.offset}}
#' @examples
#' #example 1
#' #calcNOAASunrise
#' d <- seq(ISOdatetime(2000,1,31,0,0,0), by = "hour", length.out = 24)
#' ans <- calcNOAASunrise(d)
#' #solar elevation vs hour of day at (0,0) on 2000-01-31
#' plot(solar.elevation.corrected.for.atm.refraction.deg~time.stamp,
#' main="2000-01-31; latitude 0, longitude 0",
#' data=ans, type="b")
#'
#' #anotate further, add dusk and dawn
#' #NB: sunrise and sunset are given as fractions of day
#' dusk <- ISOdatetime(2000,1,31,0,0,0) + (ans$sunset.time.lst*86400)
#' dawn <- ISOdatetime(2000,1,31,0,0,0) + (ans$sunrise.time.lst*86400)
#'
#' abline(v=dusk, lty=2)
#' abline(v=dawn, lty=2)
#' abline(h=0, lty=3)
#' arrows(dusk, 25, dusk, 0, col="blue")
#' text(dusk, 25, "dusk", col="blue")
#' arrows(dawn, 25, dawn, 0, col="blue")
#' text(dawn, 25, "dawn", col="blue")
############################
##calcNOAASunrise
############################
#daylight parameters
#kr v 0.3 2013/08/28
##################
#suggestions:
##################
#local hour offset could be a lookup table linked to tz
#### #' @importFrom devtools as.package
#' @rdname noaa.sunrise.code
#' @export
#to do
##########################
#comments
##########################
calcNOAASunrise <- function(time.stamp,
latitude = 0, longitude = 0, local.hour.offset = 0,
output = "all"){
##################
#check OK to start
##################
if(!"POSIXt" %in% class(time.stamp))
stop(paste("In noaaCalc(...) 'time.stamp' not required 'POSIXt' class",
sep=""), call.=FALSE)
if(!is.numeric(latitude) | !is.numeric(longitude))
stop(paste("In noaaCalc(...) 'latitude' and longitude' must both be numeric",
sep=""), call.=FALSE)
if(is.character(local.hour.offset))
local.hour.offset <- hour.offset(local.hour.offset)
if(!is.numeric(local.hour.offset))
stop(paste("In noaaCalc(...) 'local.hour.offset' must both be numeric",
sep=""), call.=FALSE)
##################
#local functions
##################
rad <- function(x) x * pi / 180
degrees <- function(x) x * (180 / pi)
#################
#make julian.refs
#################
#ref Gregorian calendar back extrapolated.
#assumed good for years between 1800 and 2100
p.day <- (as.numeric(format(time.stamp, "%H")) * 3600) +
(as.numeric(format(time.stamp, "%M")) * 60) +
as.numeric(format(time.stamp, "%S"))
p.day <- p.day/86400
julian.day <- as.numeric(as.Date(time.stamp, format= "%m/%d/%Y")) +
2440587.5 + p.day - (local.hour.offset/24)
julian.century <- (julian.day-2451545)/36525
################
#main calculations
################
geom.mean.long.sun.deg <- (280.46646 + julian.century * (36000.76983 + julian.century * 0.0003032)) %% 360
geom.mean.anom.sun.deg <- 357.52911 + julian.century * (35999.05029 - 0.0001537 * julian.century)
eccent.earth.orbit <- 0.016708634 - julian.century * (0.000042037 + 0.0001537 * julian.century)
sun.eq.of.ctr <- sin(rad(geom.mean.anom.sun.deg)) *
(1.914602 - julian.century *
(0.004817 + 0.000014*julian.century)) +
sin(rad(2 * geom.mean.anom.sun.deg)) *
(0.019993 - 0.000101*julian.century) +
sin(rad(3 * geom.mean.anom.sun.deg)) * 0.000289
sun.true.long.deg <- sun.eq.of.ctr + geom.mean.long.sun.deg
sun.true.anom.deg <- sun.eq.of.ctr + geom.mean.anom.sun.deg
sun.rad.vector.aus <- (1.000001018 *
(1 - eccent.earth.orbit * eccent.earth.orbit)) /
(1 + eccent.earth.orbit *
cos(rad(sun.true.anom.deg)))
sun.app.long.deg <- sun.true.long.deg - 0.00569 - 0.00478 *
sin(rad(125.04 - 1934.136 * julian.century))
mean.obliq.ecliptic.deg <- 23 + (26 + ((21.448 - julian.century *
(46.815 + julian.century *
(0.00059 - julian.century
* 0.001813)))) / 60) / 60
obliq.corr.deg <- mean.obliq.ecliptic.deg +
0.00256 * cos(rad(125.04 - 1934.136 * julian.century))
sun.rt.ascen.deg <- degrees(atan2(cos(rad(sun.app.long.deg)),
cos(rad(obliq.corr.deg)) * sin(rad(sun.app.long.deg))))
sun.declin.deg <- degrees(asin(sin(rad(obliq.corr.deg)) *
sin(rad(sun.app.long.deg))))
vary <- tan(rad(obliq.corr.deg / 2)) * tan(rad(obliq.corr.deg/2))
eq.of.time.minutes <- 4 * degrees(vary * sin(2 * rad(geom.mean.long.sun.deg)) -
2 * eccent.earth.orbit * sin(rad(geom.mean.anom.sun.deg)) +
4 * eccent.earth.orbit * vary * sin(rad(geom.mean.anom.sun.deg)) *
cos(2 * rad(geom.mean.long.sun.deg)) - 0.5 * vary * vary *
sin(4 * rad(geom.mean.long.sun.deg)) - 1.25 * eccent.earth.orbit *
eccent.earth.orbit * sin(2 * rad(geom.mean.anom.sun.deg)))
#####################
#modified selection
#####################
#original nooa code
#####################
#ha.sunrise.deg <- degrees(acos(cos(rad(90.833)) /
# (cos(rad(latitude)) * cos(rad(sun.declin.deg))) -
# tan(rad(latitude)) * tan(rad(sun.declin.deg))))
#####################
#R error catcher
#for long nights (> 24hours) and short nights (<0)
#####################
#need to test this at extremes
#####################
ha.sunrise.deg <- cos(rad(90.833)) /
(cos(rad(latitude)) * cos(rad(sun.declin.deg))) -
tan(rad(latitude)) * tan(rad(sun.declin.deg))
ha.sunrise.deg <- ifelse(ha.sunrise.deg > 1, 1, ha.sunrise.deg)
ha.sunrise.deg <- ifelse(ha.sunrise.deg < -1, -1, ha.sunrise.deg)
ha.sunrise.deg <- degrees(acos(ha.sunrise.deg))
solar.noon.lst <- (720 - 4 * longitude - eq.of.time.minutes + local.hour.offset * 60) / 1440
sunrise.time.lst <- solar.noon.lst - ha.sunrise.deg * 4 / 1440
sunset.time.lst <- solar.noon.lst + ha.sunrise.deg * 4 / 1440
sunlight.duration.minutes <- 8 * ha.sunrise.deg
true.solar.time.min <- (p.day * 1440 + eq.of.time.minutes + 4 * longitude - 60 * local.hour.offset) %% 1440
hour.angle.deg <- ifelse(true.solar.time.min/4<0,
true.solar.time.min/4 + 180,
true.solar.time.min/4 - 180)
solar.zenith.angle.deg <- degrees(acos(sin(rad(latitude)) *
sin(rad(sun.declin.deg)) + cos(rad(latitude)) *
cos(rad(sun.declin.deg)) * cos(rad(hour.angle.deg))))
solar.elevation.angle.deg <- 90 - solar.zenith.angle.deg
approx.atmospheric.refraction.deg <- ifelse(solar.elevation.angle.deg>85, 0, ifelse(solar.elevation.angle.deg>5,
58.1 / tan(rad(solar.elevation.angle.deg)) - 0.07 /
(tan(rad(solar.elevation.angle.deg)))^3 +
0.000086 / (tan(rad(solar.elevation.angle.deg)))^5,
ifelse(solar.elevation.angle.deg> -0.575,
1735 + solar.elevation.angle.deg *
(-518.2 + solar.elevation.angle.deg *
(103.4 + solar.elevation.angle.deg *
(-12.79 + solar.elevation.angle.deg * 0.711))),
-20.772 / tan(rad(solar.elevation.angle.deg)))
))/3600
solar.elevation.corrected.for.atm.refraction.deg <- solar.elevation.angle.deg + approx.atmospheric.refraction.deg
solar.azimuth.angle.deg.cw.from.n <- ifelse(hour.angle.deg>0,
#if modolo
(degrees(acos(((sin(rad(latitude)) *
cos(rad(solar.zenith.angle.deg))) -
sin(rad(sun.declin.deg))) / (cos(rad(latitude)) *
sin(rad(solar.zenith.angle.deg))))) + 180) %% 360,
#else modolo
#could rationalise?
(540 - degrees(acos(((sin(rad(latitude)) *
cos(rad(solar.zenith.angle.deg))) -
sin(rad(sun.declin.deg))) / (cos(rad(latitude)) *
sin(rad(solar.zenith.angle.deg)))))) %% 360)
#################################
#need to confirm dusk/dawn handing
#################################
daylight <- ifelse(sunlight.duration.minutes==0, FALSE,
ifelse(sunlight.duration.minutes==1440, TRUE,
ifelse(sunrise.time.lst<sunset.time.lst,
ifelse(p.day < sunset.time.lst & p.day > sunrise.time.lst, TRUE, FALSE),
ifelse(p.day <= sunrise.time.lst & p.day >= sunset.time.lst, FALSE, TRUE)
)))
#as ordered factor
daylight <- factor(daylight, levels=c(TRUE, FALSE), labels=c("daylight", "nighttime"))
###################
#NOTES
###################
#solar.noon.lst is a faction of day
#seconds into that day = p.time * 86400
#so for example sunset time is
#as.POSIXct(sunset.time.lst * 86400, origin = format(time.stamp, "%Y-%m-%d"))
###################
#currently unsure about extremes
#long nights and days at poles
####################
#output
####################
ans <- (data.frame(time.stamp = time.stamp, latitude = latitude,
longitude = longitude, local.hour.offset = local.hour.offset,
julian.day = julian.day, julian.century = julian.century,
geom.mean.long.sun.deg = geom.mean.long.sun.deg,
geom.mean.anom.sun.deg = geom.mean.anom.sun.deg,
eccent.earth.orbit = eccent.earth.orbit, sun.eq.of.ctr = sun.eq.of.ctr,
sun.true.long.deg = sun.true.long.deg, sun.true.anom.deg = sun.true.anom.deg,
sun.rad.vector.aus = sun.rad.vector.aus, sun.app.long.deg = sun.app.long.deg,
mean.obliq.ecliptic.deg = mean.obliq.ecliptic.deg,
obliq.corr.deg = obliq.corr.deg, sun.rt.ascen.deg = sun.rt.ascen.deg,
sun.declin.deg = sun.declin.deg, vary = vary,
eq.of.time.minutes = eq.of.time.minutes, ha.sunrise.deg = ha.sunrise.deg,
solar.noon.lst = solar.noon.lst, sunrise.time.lst = sunrise.time.lst,
sunset.time.lst = sunset.time.lst,
sunlight.duration.minutes = sunlight.duration.minutes,
true.solar.time.min = true.solar.time.min, hour.angle.deg = hour.angle.deg,
solar.zenith.angle.deg = solar.zenith.angle.deg,
solar.elevation.angle.deg = solar.elevation.angle.deg,
approx.atmospheric.refraction.deg = approx.atmospheric.refraction.deg,
solar.elevation.corrected.for.atm.refraction.deg = solar.elevation.corrected.for.atm.refraction.deg,
solar.azimuth.angle.deg.cw.from.n = solar.azimuth.angle.deg.cw.from.n,
daylight=daylight
))
if(is.character(output))
if(length(output)==1 && output=="all")
return(ans) else {
temp <- try(ans[output], silent=TRUE)
if(!is(temp)[1]=="try-error")
return(temp)
}
warning(paste("In noaaCalc(...) 'output' ignored because value not understood",
sep=""), call.=FALSE)
ans
}
############################
##cutDaylight
############################
#condition openair data by daylight
#using date (POSIXt)
#kr v 0.2
#based on noaa methods
#http://www.srrb.noaa.gov/highlights/sunrise/calcdetails.html
#by Chris Cornwall, Aaron Horiuchi and Chris Lehman
######################
#notes
######################
#calculations use
#(lat, long) position relative to sun
#to estimate if daylight or nighttime hour
######################
#solar.noon.lst, etc are factions of day
#seconds into that day = p.time * 86400
#so for example sunset time is
#as.POSIXct(sunset.time.lst * 86400, origin = format(mydata$date, "%Y-%m-%d"))
#(assuming you do not run into next day!)
######################
#currently unsure about extremes
#long nights and days at poles need checking
#
##################
#suggestions:
##################
#local hour offset could be a lookup table linked to tz
#
#### #' @importFrom devtools as.package
#' @rdname noaa.sunrise.code
#' @export
cutDaylight <- function(mydata, latitude = 0, longitude = 0,
local.hour.offset = 0){
##################
#check OK to start
##################
if(!"POSIXt" %in% class(mydata$date))
stop("required field 'date' missing or not POSIXt\n", call. = FALSE)
if(is.character(local.hour.offset))
local.hour.offset <- hour.offset(local.hour.offset)
###################
#temp functions
###################
rad <- function(x) x * pi / 180
degrees <- function(x) x * (180 / pi)
###############
#get local time
###############
temp <- mydata$date
#################
#make julian.refs
#################
#ref Gregorian calendar back extrapolated.
#assumed good for years between 1800 and 2100
p.day <- (as.numeric(format(temp, "%H")) * 3600) +
(as.numeric(format(temp, "%M")) * 60) +
as.numeric(format(temp, "%S"))
p.day <- p.day/86400
#julian century (via julian day)
julian.century <- as.numeric(as.Date(temp, format= "%m/%d/%Y")) + 2440587.5 + p.day - (local.hour.offset/24)
julian.century <- (julian.century-2451545)/36525
##################
#main calcs
##################
#as of noaa
geom.mean.long.sun.deg <- (280.46646 + julian.century *(36000.76983 + julian.century * 0.0003032)) %% 360
geom.mean.anom.sun.deg <- 357.52911 + julian.century * (35999.05029 - 0.0001537 * julian.century)
eccent.earth.orbit <- 0.016708634 - julian.century * (0.000042037 + 0.0001537 * julian.century)
sun.eq.of.ctr <- sin(rad(geom.mean.anom.sun.deg)) *
(1.914602 - julian.century * (0.004817 + 0.000014*julian.century)) +
sin(rad(2 * geom.mean.anom.sun.deg)) *
(0.019993 - 0.000101*julian.century) +
sin(rad(3 * geom.mean.anom.sun.deg)) * 0.000289
sun.true.long.deg <- sun.eq.of.ctr + geom.mean.long.sun.deg
sun.app.long.deg <- sun.true.long.deg - 0.00569 - 0.00478 *
sin(rad(125.04 - 1934.136 * julian.century))
mean.obliq.ecliptic.deg <- 23 + (26 + ((21.448 - julian.century *
(46.815 + julian.century *
(0.00059 - julian.century
* 0.001813)))) / 60) / 60
obliq.corr.deg <- mean.obliq.ecliptic.deg +
0.00256 * cos(rad(125.04 - 1934.136 * julian.century))
sun.declin.deg <- degrees(asin(sin(rad(obliq.corr.deg)) *
sin(rad(sun.app.long.deg))))
vary <- tan(rad(obliq.corr.deg / 2)) * tan(rad(obliq.corr.deg/2))
eq.of.time.minutes <- 4 * degrees(vary * sin(2 * rad(geom.mean.long.sun.deg)) -
2 * eccent.earth.orbit * sin(rad(geom.mean.anom.sun.deg)) +
4 * eccent.earth.orbit * vary * sin(rad(geom.mean.anom.sun.deg)) *
cos(2 * rad(geom.mean.long.sun.deg)) - 0.5 * vary * vary *
sin(4 * rad(geom.mean.long.sun.deg)) - 1.25 * eccent.earth.orbit *
eccent.earth.orbit * sin(2 * rad(geom.mean.anom.sun.deg)))
#original nooa code
##
#ha.sunrise.deg <- degrees(acos(cos(rad(90.833)) /
# (cos(rad(latitude)) * cos(rad(sun.declin.deg))) -
# tan(rad(latitude)) * tan(rad(sun.declin.deg))))
##
#R error catcher added
#for long nights>24hours/short nights<0
ha.sunrise.deg <- cos(rad(90.833)) /
(cos(rad(latitude)) * cos(rad(sun.declin.deg))) -
tan(rad(latitude)) * tan(rad(sun.declin.deg))
ha.sunrise.deg <- ifelse(ha.sunrise.deg > 1, 1, ha.sunrise.deg)
ha.sunrise.deg <- ifelse(ha.sunrise.deg < -1, -1, ha.sunrise.deg)
ha.sunrise.deg <- degrees(acos(ha.sunrise.deg))
solar.noon.lst <- (720 - 4 * longitude - eq.of.time.minutes + local.hour.offset * 60) / 1440
sunrise.time.lst <- solar.noon.lst - ha.sunrise.deg * 4 / 1440
sunset.time.lst <- solar.noon.lst + ha.sunrise.deg * 4 / 1440
sunlight.duration.minutes <- 8 * ha.sunrise.deg
#################################
#daylight factor
#################################
#need to confirm dusk/dawn handing
daylight <- ifelse(sunlight.duration.minutes==0, FALSE,
ifelse(sunlight.duration.minutes==1440, TRUE,
ifelse(sunrise.time.lst<sunset.time.lst,
ifelse(p.day < sunset.time.lst & p.day > sunrise.time.lst, TRUE, FALSE),
ifelse(p.day <= sunrise.time.lst & p.day >= sunset.time.lst, FALSE, TRUE)
)))
#as ordered factor
daylight <- factor(daylight, levels=c(TRUE, FALSE),
labels=c("daylight", "nighttime"))
###############################
#output
###############################
mydata <- cbind(mydata, daylight = daylight)
}
############################
##hour.offset
############################
#simple lookup using ref.hour.offset
#kr v 0.3 2013/08/28
#################
#notes
#################
#################
#pmatch actual grep(local, front of ref)
#for leading edge matching
#### #' @importFrom devtools as.package
#' @rdname noaa.sunrise.code
#' @export
hour.offset <- function(local, pmatch=TRUE, ignore.case=TRUE){
if(missing(local))
stop("No 'local' set. \n\t[?hour.offset for details]",
call.=FALSE)
if(!is.character(local))
stop("invalid 'local' class. \n\t[numeric or character. see ?hour.offset]",
call.=FALSE)
if(ignore.case)
local <- toupper(local)
temp <- if(pmatch)
grep(local, substr(ref.hour.offset$local,1,nchar(local))) else
match(local, ref.hour.offset$local)
if(length(temp)<1)
stop("'local' not known. \n\t[see ?hour.offset]", call.=FALSE)
if(length(temp)>1)
stop("'local' not unique.\n\t[suggest one of: "
,paste(ref.hour.offset$local[temp], collapse=", "),
"]\n\t[or see ?hour.offset]", call.=FALSE)
ref.hour.offset$hour.offset[temp]
}
karlropkins/grey.area documentation built on Dec. 27, 2024, 7:10 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions hour.offset cutDaylight calcNOAASunrise
Documented in calcNOAASunrise cutDaylight hour.offset
####################################################
#' @title NOAA Sunrise Calculations
####################################################
#'
#' @name noaa.sunrise.code
#' @aliases noaa.sunrise calcNOAASunrise cutDaylight
#' hour.offset
#' @description Methods for the calculation of sunrise
#' and sunset times based on NOAA methods.
#' @param time.stamp (POSIX) A vector of date/time of
#' \code{POSIX*} class.
#' @param latitude,longitude (Numerics) The latitude and
#' longitude of the location where the associated \code{time.stamp}
#' measurement was logged. (Note: \code{latitude} + to North; - to
#' South. \code{longitude} + to East; - to West.)
#' @param local.hour.offset (Numeric or character) The associated
#' local time offset relative to GMT in hours (numeric) or a reference
#' term that \code{\link{hour.offset}} can convert into a numeric time
#' offset (character). (Note: + to East; - to West).
#' @param output (Character) Either the default \code{"all"}, which
#' returns all inputs and calculated parameters as a data.frame, or a
#' character vector of field/column names specifying which fields/columns
#' to return. See Value below for summary of calculated parameters
#' generated by this function.
#' @param mydata (Data.frame). A date frame which includes the
#' \code{time.stamp} vector as a field/column called \code{date}.
#' @param local (Character) The location (or reference) for required
#' hour offset.
#' @param pmatch,ignore.case (Logicals) Supplied \code{local} handling.
#' \code{pmatch} applies partial (leading edge) matching of reference
#' \code{local} terms. \code{ignore.case} resets \code{local} to upper
#' case before comparing with reference \code{local} terms. Both defaults
#' \code{TRUE}.
#' @details \code{calcNOAASunrise} is based on NOAA methods and calculates
#' a range of parameters based on the relative position of the Sun and a
#' receptor point of the Earth`s surface (see value below). The method is
#' valid for dates between 1901 and 2099.\code{cutDaylight} is a cut-down
#' version of the same method that was used in the \code{openair} package.
#'
#' \code{hour.offset} uses \code{\link{ref.hour.offset}} as a look-up
#' table to match supplied \code{local} terms and return
#' \code{hour.offset} values for use in \code{\link{calcNOAASunrise}}
#' or \code{\link{cutDaylight}}.
#' @return \code{calcNOAASunrise} returns a data frame which, depending
#' on the call \code{output} value, may include the call arguments and
#' any of the following calculated parameters:
#' \describe{
#' \item{julian.day }{(Numeric) The Julian day.}
#' \item{julian.century }{(Numeric) The Julian century.}
#' \item{geom.mean.long.sun.deg }{(Numeric) The geometric mean longitude
#' of the Sun in degrees.}
#' \item{geom.mean.anom.sun.deg }{(Numeric) The geometric mean anomaly of
#' the Sun in degrees.}
#' \item{eccent.earth.orbit }{(Numeric) The eccentricity of Earth`s orbit.}
#' \item{sun.eq.of.ctr }{(Numeric) The estimated center for the Sun, in
#' degrees.}
#' \item{sun.true.long.deg }{(Numeric) The true longitude of the Sun, in
#' degrees.}
#' \item{sun.true.anom.deg }{(Numeric) The true anamoly of the Sun, in
#' degrees.}
#' \item{sun.rad.vector.aus }{(Numeric) The distance to the Sun, in AUs
#' (radial).}
#' \item{sun.app.long.deg }{(Numeric) The apparent longitude of the Sun,
#' in degrees.}
#' \item{mean.obliq.ecliptic.deg }{(Numeric) The mean obliquity of the
#' ecliptic, in degrees.}
#' \item{obliq.corr.deg }{(Numeric) The corrected obliquity of the
#' ecliptic, in degrees.}
#' \item{sun.rt.ascen.deg }{(Numeric) The right ascension of the Sun,
#' in degress.}
#' \item{sun.declin.deg }{(Numeric) The declination of the Sun, in
#' degrees.}
#' \item{vary }{(Numeric) TO BE CONFIRMED}
#' \item{eq.of.time.minutes }{(Numeric) TO BE CONFIRMED}
#' \item{ha.sunrise.deg }{(Numeric) TO BE CONFIRMED}
#' \item{solar.noon.lst }{(Numeric) TO BE CONFIRMED}
#' \item{sunrise.time.lst }{(Numeric) The estimated sunrise time as
#' fraction of day.}
#' \item{sunset.time.lst }{(Numeric) The estimated sunset time as
#' fraction of day.}
#' \item{sunlight.duration.minutes }{(Numeric) TO BE CONFIRMED}
#' \item{true.solar.time.min }{(Numeric) TO BE CONFIRMED}
#' \item{hour.angle.deg }{(Numeric) TO BE CONFIRMED}
#' \item{solar.zenith.angle.deg }{(Numeric) TO BE CONFIRMED}
#' \item{solar.elevation.angle.deg }{(Numeric) TO BE CONFIRMED}
#' \item{approx.atmospheric.refraction.deg }{(Numeric) TO BE CONFIRMED}
#' \item{solar.elevation.corrected.for.atm.refraction.deg }{(Numeric) TO
#' BE CONFIRMED}
#' \item{solar.azimuth.angle.deg.cw.from.n }{(Numeric) TO BE CONFIRMED}
#' \item{daylight }{(Factor) A factor specifying if the \code{time.stamp},
#' \code{latitude}, \code{longitude} and \code{local.hour.offset} is
#' estimated to be a \code{daylight} or \code{nighttime} instance.}
#' }
#' \code{cutDaylight} is an alternative to \code{calcNOAASunrise} which
#' accepts the time stamp as part of a data frame and returns that data
#' frame with the \code{calcNOAACalc} output component \code{daylight}
#' as an addition field/column.
#'
#' \code{hour.offset} returns the hour.offset as a numeric.
#' @references Methods based on:
#' \url{http://www.srrb.noaa.gov/highlights/sunrise/sunrise.html}
#' @author Karl Ropkins
#' @note The NOAA Sunrise methods were originally recoded for R as part
#' of work on the NERC penair project:
#' \url{http://www.openair-project.org/}.
#' \code{cutDaylight} is an alternative version of the function. It is
#' highly similar to code used as part of the function \code{cutData} in
#' the package \code{openair} to provide \code{daylight} conditioning for
#' many plot types within that package.
#'
#' Character handling for \code{local.hour.offset} is by
#' \code{hour.offset} and \code{\link{ref.hour.offset}}.
#' @seealso \code{\link{ref.hour.offset}}
#' @examples
#' #example 1
#' #calcNOAASunrise
#' d <- seq(ISOdatetime(2000,1,31,0,0,0), by = "hour", length.out = 24)
#' ans <- calcNOAASunrise(d)
#' #solar elevation vs hour of day at (0,0) on 2000-01-31
#' plot(solar.elevation.corrected.for.atm.refraction.deg~time.stamp,
#' main="2000-01-31; latitude 0, longitude 0",
#' data=ans, type="b")
#'
#' #anotate further, add dusk and dawn
#' #NB: sunrise and sunset are given as fractions of day
#' dusk <- ISOdatetime(2000,1,31,0,0,0) + (ans$sunset.time.lst*86400)
#' dawn <- ISOdatetime(2000,1,31,0,0,0) + (ans$sunrise.time.lst*86400)
#'
#' abline(v=dusk, lty=2)
#' abline(v=dawn, lty=2)
#' abline(h=0, lty=3)
#' arrows(dusk, 25, dusk, 0, col="blue")
#' text(dusk, 25, "dusk", col="blue")
#' arrows(dawn, 25, dawn, 0, col="blue")
#' text(dawn, 25, "dawn", col="blue")
############################
##calcNOAASunrise
############################
#daylight parameters
#kr v 0.3 2013/08/28
##################
#suggestions:
##################
#local hour offset could be a lookup table linked to tz
#### #' @importFrom devtools as.package
#' @rdname noaa.sunrise.code
#' @export
#to do
##########################
#comments
##########################
calcNOAASunrise <- function(time.stamp,
latitude = 0, longitude = 0, local.hour.offset = 0,
output = "all"){
##################
#check OK to start
##################
if(!"POSIXt" %in% class(time.stamp))
stop(paste("In noaaCalc(...) 'time.stamp' not required 'POSIXt' class",
sep=""), call.=FALSE)
if(!is.numeric(latitude) | !is.numeric(longitude))
stop(paste("In noaaCalc(...) 'latitude' and longitude' must both be numeric",
sep=""), call.=FALSE)
if(is.character(local.hour.offset))
local.hour.offset <- hour.offset(local.hour.offset)
if(!is.numeric(local.hour.offset))
stop(paste("In noaaCalc(...) 'local.hour.offset' must both be numeric",
sep=""), call.=FALSE)
##################
#local functions
##################
rad <- function(x) x * pi / 180
degrees <- function(x) x * (180 / pi)
#################
#make julian.refs
#################
#ref Gregorian calendar back extrapolated.
#assumed good for years between 1800 and 2100
p.day <- (as.numeric(format(time.stamp, "%H")) * 3600) +
(as.numeric(format(time.stamp, "%M")) * 60) +
as.numeric(format(time.stamp, "%S"))
p.day <- p.day/86400
julian.day <- as.numeric(as.Date(time.stamp, format= "%m/%d/%Y")) +
2440587.5 + p.day - (local.hour.offset/24)
julian.century <- (julian.day-2451545)/36525
################
#main calculations
################
geom.mean.long.sun.deg <- (280.46646 + julian.century * (36000.76983 + julian.century * 0.0003032)) %% 360
geom.mean.anom.sun.deg <- 357.52911 + julian.century * (35999.05029 - 0.0001537 * julian.century)
eccent.earth.orbit <- 0.016708634 - julian.century * (0.000042037 + 0.0001537 * julian.century)
sun.eq.of.ctr <- sin(rad(geom.mean.anom.sun.deg)) *
(1.914602 - julian.century *
(0.004817 + 0.000014*julian.century)) +
sin(rad(2 * geom.mean.anom.sun.deg)) *
(0.019993 - 0.000101*julian.century) +
sin(rad(3 * geom.mean.anom.sun.deg)) * 0.000289
sun.true.long.deg <- sun.eq.of.ctr + geom.mean.long.sun.deg
sun.true.anom.deg <- sun.eq.of.ctr + geom.mean.anom.sun.deg
sun.rad.vector.aus <- (1.000001018 *
(1 - eccent.earth.orbit * eccent.earth.orbit)) /
(1 + eccent.earth.orbit *
cos(rad(sun.true.anom.deg)))
sun.app.long.deg <- sun.true.long.deg - 0.00569 - 0.00478 *
sin(rad(125.04 - 1934.136 * julian.century))
mean.obliq.ecliptic.deg <- 23 + (26 + ((21.448 - julian.century *
(46.815 + julian.century *
(0.00059 - julian.century
* 0.001813)))) / 60) / 60
obliq.corr.deg <- mean.obliq.ecliptic.deg +
0.00256 * cos(rad(125.04 - 1934.136 * julian.century))
sun.rt.ascen.deg <- degrees(atan2(cos(rad(sun.app.long.deg)),
cos(rad(obliq.corr.deg)) * sin(rad(sun.app.long.deg))))
sun.declin.deg <- degrees(asin(sin(rad(obliq.corr.deg)) *
sin(rad(sun.app.long.deg))))
vary <- tan(rad(obliq.corr.deg / 2)) * tan(rad(obliq.corr.deg/2))
eq.of.time.minutes <- 4 * degrees(vary * sin(2 * rad(geom.mean.long.sun.deg)) -
2 * eccent.earth.orbit * sin(rad(geom.mean.anom.sun.deg)) +
4 * eccent.earth.orbit * vary * sin(rad(geom.mean.anom.sun.deg)) *
cos(2 * rad(geom.mean.long.sun.deg)) - 0.5 * vary * vary *
sin(4 * rad(geom.mean.long.sun.deg)) - 1.25 * eccent.earth.orbit *
eccent.earth.orbit * sin(2 * rad(geom.mean.anom.sun.deg)))
#####################
#modified selection
#####################
#original nooa code
#####################
#ha.sunrise.deg <- degrees(acos(cos(rad(90.833)) /
# (cos(rad(latitude)) * cos(rad(sun.declin.deg))) -
# tan(rad(latitude)) * tan(rad(sun.declin.deg))))
#####################
#R error catcher
#for long nights (> 24hours) and short nights (<0)
#####################
#need to test this at extremes
#####################
ha.sunrise.deg <- cos(rad(90.833)) /
(cos(rad(latitude)) * cos(rad(sun.declin.deg))) -
tan(rad(latitude)) * tan(rad(sun.declin.deg))
ha.sunrise.deg <- ifelse(ha.sunrise.deg > 1, 1, ha.sunrise.deg)
ha.sunrise.deg <- ifelse(ha.sunrise.deg < -1, -1, ha.sunrise.deg)
ha.sunrise.deg <- degrees(acos(ha.sunrise.deg))
solar.noon.lst <- (720 - 4 * longitude - eq.of.time.minutes + local.hour.offset * 60) / 1440
sunrise.time.lst <- solar.noon.lst - ha.sunrise.deg * 4 / 1440
sunset.time.lst <- solar.noon.lst + ha.sunrise.deg * 4 / 1440
sunlight.duration.minutes <- 8 * ha.sunrise.deg
true.solar.time.min <- (p.day * 1440 + eq.of.time.minutes + 4 * longitude - 60 * local.hour.offset) %% 1440
hour.angle.deg <- ifelse(true.solar.time.min/4<0,
true.solar.time.min/4 + 180,
true.solar.time.min/4 - 180)
solar.zenith.angle.deg <- degrees(acos(sin(rad(latitude)) *
sin(rad(sun.declin.deg)) + cos(rad(latitude)) *
cos(rad(sun.declin.deg)) * cos(rad(hour.angle.deg))))
solar.elevation.angle.deg <- 90 - solar.zenith.angle.deg
approx.atmospheric.refraction.deg <- ifelse(solar.elevation.angle.deg>85, 0, ifelse(solar.elevation.angle.deg>5,
58.1 / tan(rad(solar.elevation.angle.deg)) - 0.07 /
(tan(rad(solar.elevation.angle.deg)))^3 +
0.000086 / (tan(rad(solar.elevation.angle.deg)))^5,
ifelse(solar.elevation.angle.deg> -0.575,
1735 + solar.elevation.angle.deg *
(-518.2 + solar.elevation.angle.deg *
(103.4 + solar.elevation.angle.deg *
(-12.79 + solar.elevation.angle.deg * 0.711))),
-20.772 / tan(rad(solar.elevation.angle.deg)))
))/3600
solar.elevation.corrected.for.atm.refraction.deg <- solar.elevation.angle.deg + approx.atmospheric.refraction.deg
solar.azimuth.angle.deg.cw.from.n <- ifelse(hour.angle.deg>0,
#if modolo
(degrees(acos(((sin(rad(latitude)) *
cos(rad(solar.zenith.angle.deg))) -
sin(rad(sun.declin.deg))) / (cos(rad(latitude)) *
sin(rad(solar.zenith.angle.deg))))) + 180) %% 360,
#else modolo
#could rationalise?
(540 - degrees(acos(((sin(rad(latitude)) *
cos(rad(solar.zenith.angle.deg))) -
sin(rad(sun.declin.deg))) / (cos(rad(latitude)) *
sin(rad(solar.zenith.angle.deg)))))) %% 360)
#################################
#need to confirm dusk/dawn handing
#################################
daylight <- ifelse(sunlight.duration.minutes==0, FALSE,
ifelse(sunlight.duration.minutes==1440, TRUE,
ifelse(sunrise.time.lst<sunset.time.lst,
ifelse(p.day < sunset.time.lst & p.day > sunrise.time.lst, TRUE, FALSE),
ifelse(p.day <= sunrise.time.lst & p.day >= sunset.time.lst, FALSE, TRUE)
)))
#as ordered factor
daylight <- factor(daylight, levels=c(TRUE, FALSE), labels=c("daylight", "nighttime"))
###################
#NOTES
###################
#solar.noon.lst is a faction of day
#seconds into that day = p.time * 86400
#so for example sunset time is
#as.POSIXct(sunset.time.lst * 86400, origin = format(time.stamp, "%Y-%m-%d"))
###################
#currently unsure about extremes
#long nights and days at poles
####################
#output
####################
ans <- (data.frame(time.stamp = time.stamp, latitude = latitude,
longitude = longitude, local.hour.offset = local.hour.offset,
julian.day = julian.day, julian.century = julian.century,
geom.mean.long.sun.deg = geom.mean.long.sun.deg,
geom.mean.anom.sun.deg = geom.mean.anom.sun.deg,
eccent.earth.orbit = eccent.earth.orbit, sun.eq.of.ctr = sun.eq.of.ctr,
sun.true.long.deg = sun.true.long.deg, sun.true.anom.deg = sun.true.anom.deg,
sun.rad.vector.aus = sun.rad.vector.aus, sun.app.long.deg = sun.app.long.deg,
mean.obliq.ecliptic.deg = mean.obliq.ecliptic.deg,
obliq.corr.deg = obliq.corr.deg, sun.rt.ascen.deg = sun.rt.ascen.deg,
sun.declin.deg = sun.declin.deg, vary = vary,
eq.of.time.minutes = eq.of.time.minutes, ha.sunrise.deg = ha.sunrise.deg,
solar.noon.lst = solar.noon.lst, sunrise.time.lst = sunrise.time.lst,
sunset.time.lst = sunset.time.lst,
sunlight.duration.minutes = sunlight.duration.minutes,
true.solar.time.min = true.solar.time.min, hour.angle.deg = hour.angle.deg,
solar.zenith.angle.deg = solar.zenith.angle.deg,
solar.elevation.angle.deg = solar.elevation.angle.deg,
approx.atmospheric.refraction.deg = approx.atmospheric.refraction.deg,
solar.elevation.corrected.for.atm.refraction.deg = solar.elevation.corrected.for.atm.refraction.deg,
solar.azimuth.angle.deg.cw.from.n = solar.azimuth.angle.deg.cw.from.n,
daylight=daylight
))
if(is.character(output))
if(length(output)==1 && output=="all")
return(ans) else {
temp <- try(ans[output], silent=TRUE)
if(!is(temp)[1]=="try-error")
return(temp)
}
warning(paste("In noaaCalc(...) 'output' ignored because value not understood",
sep=""), call.=FALSE)
ans
}
############################
##cutDaylight
############################
#condition openair data by daylight
#using date (POSIXt)
#kr v 0.2
#based on noaa methods
#http://www.srrb.noaa.gov/highlights/sunrise/calcdetails.html
#by Chris Cornwall, Aaron Horiuchi and Chris Lehman
######################
#notes
######################
#calculations use
#(lat, long) position relative to sun
#to estimate if daylight or nighttime hour
######################
#solar.noon.lst, etc are factions of day
#seconds into that day = p.time * 86400
#so for example sunset time is
#as.POSIXct(sunset.time.lst * 86400, origin = format(mydata$date, "%Y-%m-%d"))
#(assuming you do not run into next day!)
######################
#currently unsure about extremes
#long nights and days at poles need checking
#
##################
#suggestions:
##################
#local hour offset could be a lookup table linked to tz
#
#### #' @importFrom devtools as.package
#' @rdname noaa.sunrise.code
#' @export
cutDaylight <- function(mydata, latitude = 0, longitude = 0,
local.hour.offset = 0){
##################
#check OK to start
##################
if(!"POSIXt" %in% class(mydata$date))
stop("required field 'date' missing or not POSIXt\n", call. = FALSE)
if(is.character(local.hour.offset))
local.hour.offset <- hour.offset(local.hour.offset)
###################
#temp functions
###################
rad <- function(x) x * pi / 180
degrees <- function(x) x * (180 / pi)
###############
#get local time
###############
temp <- mydata$date
#################
#make julian.refs
#################
#ref Gregorian calendar back extrapolated.
#assumed good for years between 1800 and 2100
p.day <- (as.numeric(format(temp, "%H")) * 3600) +
(as.numeric(format(temp, "%M")) * 60) +
as.numeric(format(temp, "%S"))
p.day <- p.day/86400
#julian century (via julian day)
julian.century <- as.numeric(as.Date(temp, format= "%m/%d/%Y")) + 2440587.5 + p.day - (local.hour.offset/24)
julian.century <- (julian.century-2451545)/36525
##################
#main calcs
##################
#as of noaa
geom.mean.long.sun.deg <- (280.46646 + julian.century *(36000.76983 + julian.century * 0.0003032)) %% 360
geom.mean.anom.sun.deg <- 357.52911 + julian.century * (35999.05029 - 0.0001537 * julian.century)
eccent.earth.orbit <- 0.016708634 - julian.century * (0.000042037 + 0.0001537 * julian.century)
sun.eq.of.ctr <- sin(rad(geom.mean.anom.sun.deg)) *
(1.914602 - julian.century * (0.004817 + 0.000014*julian.century)) +
sin(rad(2 * geom.mean.anom.sun.deg)) *
(0.019993 - 0.000101*julian.century) +
sin(rad(3 * geom.mean.anom.sun.deg)) * 0.000289
sun.true.long.deg <- sun.eq.of.ctr + geom.mean.long.sun.deg
sun.app.long.deg <- sun.true.long.deg - 0.00569 - 0.00478 *
sin(rad(125.04 - 1934.136 * julian.century))
mean.obliq.ecliptic.deg <- 23 + (26 + ((21.448 - julian.century *
(46.815 + julian.century *
(0.00059 - julian.century
* 0.001813)))) / 60) / 60
obliq.corr.deg <- mean.obliq.ecliptic.deg +
0.00256 * cos(rad(125.04 - 1934.136 * julian.century))
sun.declin.deg <- degrees(asin(sin(rad(obliq.corr.deg)) *
sin(rad(sun.app.long.deg))))
vary <- tan(rad(obliq.corr.deg / 2)) * tan(rad(obliq.corr.deg/2))
eq.of.time.minutes <- 4 * degrees(vary * sin(2 * rad(geom.mean.long.sun.deg)) -
2 * eccent.earth.orbit * sin(rad(geom.mean.anom.sun.deg)) +
4 * eccent.earth.orbit * vary * sin(rad(geom.mean.anom.sun.deg)) *
cos(2 * rad(geom.mean.long.sun.deg)) - 0.5 * vary * vary *
sin(4 * rad(geom.mean.long.sun.deg)) - 1.25 * eccent.earth.orbit *
eccent.earth.orbit * sin(2 * rad(geom.mean.anom.sun.deg)))
#original nooa code
##
#ha.sunrise.deg <- degrees(acos(cos(rad(90.833)) /
# (cos(rad(latitude)) * cos(rad(sun.declin.deg))) -
# tan(rad(latitude)) * tan(rad(sun.declin.deg))))
##
#R error catcher added
#for long nights>24hours/short nights<0
ha.sunrise.deg <- cos(rad(90.833)) /
(cos(rad(latitude)) * cos(rad(sun.declin.deg))) -
tan(rad(latitude)) * tan(rad(sun.declin.deg))
ha.sunrise.deg <- ifelse(ha.sunrise.deg > 1, 1, ha.sunrise.deg)
ha.sunrise.deg <- ifelse(ha.sunrise.deg < -1, -1, ha.sunrise.deg)
ha.sunrise.deg <- degrees(acos(ha.sunrise.deg))
solar.noon.lst <- (720 - 4 * longitude - eq.of.time.minutes + local.hour.offset * 60) / 1440
sunrise.time.lst <- solar.noon.lst - ha.sunrise.deg * 4 / 1440
sunset.time.lst <- solar.noon.lst + ha.sunrise.deg * 4 / 1440
sunlight.duration.minutes <- 8 * ha.sunrise.deg
#################################
#daylight factor
#################################
#need to confirm dusk/dawn handing
daylight <- ifelse(sunlight.duration.minutes==0, FALSE,
ifelse(sunlight.duration.minutes==1440, TRUE,
ifelse(sunrise.time.lst<sunset.time.lst,
ifelse(p.day < sunset.time.lst & p.day > sunrise.time.lst, TRUE, FALSE),
ifelse(p.day <= sunrise.time.lst & p.day >= sunset.time.lst, FALSE, TRUE)
)))
#as ordered factor
daylight <- factor(daylight, levels=c(TRUE, FALSE),
labels=c("daylight", "nighttime"))
###############################
#output
###############################
mydata <- cbind(mydata, daylight = daylight)
}
############################
##hour.offset
############################
#simple lookup using ref.hour.offset
#kr v 0.3 2013/08/28
#################
#notes
#################
#################
#pmatch actual grep(local, front of ref)
#for leading edge matching
#### #' @importFrom devtools as.package
#' @rdname noaa.sunrise.code
#' @export
hour.offset <- function(local, pmatch=TRUE, ignore.case=TRUE){
if(missing(local))
stop("No 'local' set. \n\t[?hour.offset for details]",
call.=FALSE)
if(!is.character(local))
stop("invalid 'local' class. \n\t[numeric or character. see ?hour.offset]",
call.=FALSE)
if(ignore.case)
local <- toupper(local)
temp <- if(pmatch)
grep(local, substr(ref.hour.offset$local,1,nchar(local))) else
match(local, ref.hour.offset$local)
if(length(temp)<1)
stop("'local' not known. \n\t[see ?hour.offset]", call.=FALSE)
if(length(temp)>1)
stop("'local' not unique.\n\t[suggest one of: "
,paste(ref.hour.offset$local[temp], collapse=", "),
"]\n\t[or see ?hour.offset]", call.=FALSE)
ref.hour.offset$hour.offset[temp]
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.