Nothing
R/earthtide_class.R
In earthtide: Parallel Implementation of 'ETERNA 3.40' for Prediction and Analysis of Earth Tides
#' @title Earthtide class
#'
#' @description Class to generate synthetic earthtide signals.
#'
#' @rdname Earthtide_class
#'
#' @section Usage:
#' \preformatted{
#' et <- Earthtide$new(
#' utc = as.POSIXct("2017-01-01", tz = "UTC") + 0:(24 * 7) * 3600,
#' latitude = 52.3868,
#' longitude = 9.7144,
#' catalog = "ksm04",
#' wave_groups = data.frame(start = 0.0, end = 6.0))
#'
#' et$predict(method = "gravity", n_thread = 1)
#' et$analyze(method = "gravity", n_thread = 1)
#' et$lod_tide()
#' et$pole_tide()
#' et$tide()
#' et$print()
#' }
#'
#' @section Arguments:
#' \code{Earthtide$new}
#' \describe{
#' \item{et: }{An \code{Earthtide} object.}
#' \item{utc: }{The date-time in UTC (POSIXct vector).}
#' \item{latitude: }{The station latitude (WGS84) (degree) (numeric)
#' defaults to 0.0}
#' \item{longitude: }{The station longitude (WGS84) (degree) (numeric)
#' defaults to 0.0}
#' \item{elevation: }{The station ellipsoidal height (WGS84) (m) (numeric)
#' defaults to 0.0}
#' \item{azimuth: }{Earth azimuth (numeric) defaults to 0 (degrees)}
#' \item{gravity: }{Gravity at the station (m/s^2) (numeric) 0 to
#' estimate gravity from elevation and latitude.}
#' \item{earth_radius: }{Radius of earth (m) (numeric) defaults to 6378136.3 }
#' \item{earth_eccen: }{Eccentricity of earth (numeric)
#' defaults to 6.69439795140e-3}
#' \item{cutoff: }{Cutoff amplitude for constituents (numeric)
#' defaults to 1e-6}
#' \item{wave_groups: }{Two column data.frame having start and end of
#' frequency groups (data.frame). This data.frame must have two columns
#' with the names 'start', and 'end' signifying the start and end of the
#' wave groupings. An optional third column 'multiplier' can be provided
#' to scale the particular wave group. If column names do no match, the
#' inferred column positions are start, end, multiplier.}
#' \item{catalog: }{Use the "hw95s" catalog or "ksm04" catalog (character).}
#' \item{eop: }{User defined Earth Orientation Parameter (EOP) data.frame
#' with the following columns: datetime, ddt, ut1_utc, lod, x, y, dx, dy}
#' \item{...: }{Currently not used.}
#' }
#'
#' \code{Earthtide$predict, Earthtide$analyze}
#' \describe{
#' \item{method: }{For \code{predict} and \code{analyze}. One of "gravity",
#' "tidal_potential", "tidal_tilt", "vertical_displacement",
#' "horizontal_displacement", "n_s_displacement", "e_w_displacement",
#' "vertical_strain", "areal_strain", "volume_strain", "horizontal_strain"
#' or "ocean_tides".}
#' \item{return_matrix: }{For \code{predict} and \code{analyze}. Return a
#' matrix of tidal values instead of data.frame. The datetime column will
#' not be present in this case (logical).}
#' \item{n_thread: }{For \code{predict} and \code{analyze}. Number of threads
#' to use for parallel processing.}
#' }
#'
#' @section Details:
#'
#' \code{$new(utc, latitude, longitude, elevation, azimuth, gravity,} \cr
#' \code{earth_radius, earth_eccen, cutoff, wave_groups, catalog, ...)} \cr
#' create a new \code{Earthtide} object and initialize catalog, station
#' and times.
#'
#' \code{$predict(method, astro_argument, return_matrix)} generate a combined
#' synthetic Earth tide.
#'
#' \code{$analyze(method, astro_argument, return_matrix, scale)} generate
#' components of the Earth tide for analysis.
#'
#' \code{$lod_tide()} generate components of the LOD (Length Of Day) tide.
#'
#' \code{$pole_tide()} generate components of the pole tide.
#'
#' \code{$tide()} get the tide \code{data.frame}.
#'
#' \code{$print()} print the \code{Earthtide} object.
#'
#' @docType class
#' @aliases Earthtide-class
#' @importFrom R6 R6Class
#' @importFrom RcppThread detectCores
#' @importFrom stats approx
#' @importFrom utils read.table
#' @importFrom utils read.fwf
#' @importFrom utils download.file
#' @importFrom utils data
#' @format An \code{\link{R6Class}} generator object
#'
#'
#' @references Hartmann, T., Wenzel, H.-G., 1995. The HW95 tidal potential catalogue. Geophys. Res. Lett. 22, 3553-3556. \doi{10.1029/95GL03324}
#' @references Kudryavtsev, S.M., 2004. Improved harmonic development of the Earth tide-generating potential. J. Geod. 77, 829-838. \doi{10.1007/s00190-003-0361-2}
#' @references Wenzel, H.G., 1996. The nanogal software: Earth tide data processing package ETERNA 3.30. Bull. Inf. Marées Terrestres, 124, pp.9425-9439. \url{https://www.eas.slu.edu/GGP/ETERNA34/MANUAL/ETERNA33.HTM}
#'
#' @examples
#'
#' et <- Earthtide$new(
#' utc = as.POSIXct("2017-01-01", tz = "UTC") + 0:(24 * 7) * 3600,
#' latitude = 52.3868,
#' longitude = 9.7144,
#' catalog = "ksm04",
#' wave_groups = data.frame(start = 0.0, end = 6.0)
#' )
#'
#' et$predict(method = "gravity")
#'
#' plot(gravity ~ datetime, et$tide(), type = "l")
#'
#' @name Earthtide
NULL
#' @export
Earthtide <- R6Class(
"et",
public = list(
# initialization
initialize = function(utc,
latitude = 0,
longitude = 0,
elevation = 0,
azimuth = 0,
gravity = 0,
earth_radius = 6378136.3,
earth_eccen = 6.69439795140e-3,
cutoff = 1e-6,
wave_groups = NULL,
catalog = "ksm04",
eop = NULL,
...) {
# Initialize class using input values
self$prepare_datetime(utc, eop)
self$prepare_station(
latitude, longitude, elevation, azimuth, gravity,
earth_radius, earth_eccen
)
self$prepare_astro()
self$prepare_catalog(cutoff, wave_groups, catalog)
self$love_params <- love(latitude, elevation)
invisible(self)
},
# Initialize class using input values
prepare_datetime = function(utc, eop) {
self$datetime <- .prepare_datetime(utc, eop)
self$tides <- data.frame(datetime = utc)
},
# Calculate the astronical arguments and derivative
prepare_astro = function() {
self$astro <- .prepare_astro(self)
},
# Subset values based using a cutoff amplitude and
# cutoff frequencies
prepare_catalog = function(cutoff, wave_groups = NULL, catalog = "ksm04") {
self$catalog <- .prepare_catalog(cutoff, wave_groups, catalog = catalog)
},
# Calculate properties based on the location of station
prepare_station = function(latitude, longitude, elevation, azimuth,
gravity, earth_radius, earth_eccen) {
self$station <- .prepare_station(
self, latitude, longitude, elevation,
azimuth, gravity, earth_radius,
earth_eccen
)
self$pole_quotient <- 1.16 # 1.1788, # Ducarme et al. 2006
self$pk <- rep(0, 25)
self$delta <- rep(1.0, 25)
self$deltar <- 0.0
},
gravity = function() {
self$station$dgk <- self$station$dgz
self$pk[] <- 180
self$delta[1:12] <- self$love_params$dglat
self$deltar <- self$love_params$dgr
},
tidal_potential = function() {
self$delta[1:12] <- self$love_params$dklat
self$deltar <- self$love_params$dkr
},
tidal_tilt = function() {
to_radians <- pi / 180
cos_azimuth <- cos(self$station$azimuth * to_radians)
sin_azimuth <- sin(self$station$azimuth * to_radians)
x_comp <- self$station$dgx[1:12] * cos_azimuth
y_comp <- self$station$dgy[1:12] * sin_azimuth
self$station$dgk[1:12] <- sqrt((x_comp)^2 + (y_comp)^2) * self$station$df
wh <- which(x_comp != 0 | y_comp != 0)
self$pk[] <- 0.0
self$pk[wh] <- 180 / pi * atan2(y_comp[wh], x_comp[wh])
# from etpots
self$delta[1:12] <- self$love_params$dtlat
self$deltar <- self$love_params$dkr - self$love_params$dhr
},
vertical_displacement = function() {
dfak <- 1e3 / self$station$gravity
self$station$dgk[1:12] <- self$station$dgk[1:12] *
self$love_params$dhlat[1:12] * dfak
self$pk[] <- 0.0
},
horizontal_displacement = function() {
to_radians <- pi / 180
dfak <- 1e-6 * self$station$geo_radius / self$station$gravity
cos_azimuth <- cos(self$station$azimuth * to_radians)
sin_azimuth <- sin(self$station$azimuth * to_radians)
x_comp <- self$station$dgx[1:12] * cos_azimuth
y_comp <- self$station$dgy[1:12] * sin_azimuth
# # shida numbers
# dl0 <- c(
# 0.08400, 0.08410, 0.08520, 0.01490, 0.01490, 0.01490,
# 0.01490, 0.01000, 0.01000, 0.01000, 0.01000, 0.01000
# )
#
# if(any(abs(self$love_params$dllat[1:12]-dl0) > 0.04)) {
# self$love_params$dllat[1:12] <- dl0
# }
self$station$dgk[1:12] <- sqrt((x_comp)^2 + (y_comp)^2) *
self$love_params$dllat[1:12] * dfak
wh <- which(x_comp != 0 | y_comp != 0)
self$pk[] <- 0.0
self$pk[wh] <- 180 / pi * atan2(y_comp[wh], x_comp[wh])
# from etpots (this is a modification as I think they are 0 and 1 in ETERNA)
# look into this more closely
self$delta[1:12] <- 1.0 #self$love_params$dtlat
self$deltar <- 0.0 #self$love_params$dkr - self$love_params$dhr
},
n_s_displacement = function() {
to_radians <- pi / 180
dfak <- 1e-6 * self$station$geo_radius / self$station$gravity
cos_azimuth <- cos(0 * to_radians)
sin_azimuth <- sin(0 * to_radians)
x_comp <- self$station$dgx[1:12] * cos_azimuth
y_comp <- self$station$dgy[1:12] * sin_azimuth
self$station$dgk[1:12] <- sqrt((x_comp)^2 + (y_comp)^2) *
self$love_params$dllat[1:12] * dfak
self$pk[] <- 0.0
self$pk[1:12] <- 180 / pi * atan2(y_comp, x_comp)
wh <- which(x_comp == 0 & y_comp == 0)
self$pk[wh] <- 0.0
},
e_w_displacement = function() {
to_radians <- pi / 180
dfak <- 1e-6 * self$station$geo_radius / self$station$gravity
cos_azimuth <- cos(90 * to_radians)
sin_azimuth <- sin(90 * to_radians)
x_comp <- self$station$dgx[1:12] * cos_azimuth
y_comp <- self$station$dgy[1:12] * sin_azimuth
self$station$dgk[1:12] <- sqrt((x_comp)^2 + (y_comp)^2) *
self$love_params$dllat[1:12] * dfak
self$pk[] <- 0.0
self$pk[1:12] <- 180 / pi * atan2(y_comp, x_comp)
wh <- which(x_comp == 0 & y_comp == 0)
self$pk[wh] <- 0.0
},
vertical_strain = function(poisson = 0.25) {
dfak <- 1.e9 * poisson / (poisson - 1.0)
self$strain(dfak)
},
areal_strain = function() {
dfak <- 1.e9
self$strain(dfak)
},
volume_strain = function(poisson = 0.25) {
dfak <- 1.e9 * (1.0 - 2.0 * poisson) / (1.0 - poisson)
self$strain(dfak)
},
strain = function(dfak) {
scale <- c(rep(-6, 3), rep(-12, 4), rep(-20, 5))
self$station$dgk[1:12] <- self$station$dgk[1:12] * dfak *
(2.0 * self$love_params$dhlat[1:12] +
scale * self$love_params$dllat[1:12]) /
(self$station$gravity * self$station$geo_radius)
},
horizontal_strain = function() {
to_radians <- pi / 180
dgx <- self$station$dgx
dgy <- self$station$dgy
dhlat <- self$love_params$dhlat
dllat <- self$love_params$dllat
theta <- self$station$theta * to_radians
azr <- (self$station$azimuth + 180) * to_radians
caz <- cos(azr)
saz <- sin(azr)
saz2 <- sin(2.0 * azr)
csts <- -0.5 * sin(2.0 * azr)
ct <- sin(self$station$geo_latitude * to_radians)
st <- cos(self$station$geo_latitude * to_radians)
ct2 <- ct * ct
st2 <- st * st
cc2 <- cos(2.0 * self$station$geo_latitude * to_radians)
c2t <- -cc2
cott <- 1.0 / tan(theta)
cott2 <- 1.0 / tan(2.0 * theta)
dfak <- 1e9 / (self$station$geo_radius * self$station$gravity)
dgx[1] <-
(dhlat[1] - (6.0 * dllat[1] * c2t) / (3.0 * ct2 - 1.0)) * caz^2 +
(dhlat[1] - (6.0 * dllat[1] * ct2) / (3.0 * ct2 - 1.0)) * saz^2
dgy[1] <- 0.0
dgx[2] <- (dhlat[2] - 4.0 * dllat[2]) * caz^2 +
(dhlat[2] - dllat[2] / st2 + 2.0 * dllat[2] * cott * cott2) * saz^2
dgy[2] <- 2.0 * dllat[2] * (2.0 * cott2 - cott) * csts / st
dgx[3] <- (dhlat[3] + 2.0 * dllat[3] * (cott^2 - 1.0)) * caz^2 +
(dhlat[3] - 4.0 * dllat[3] / st2 + 2.0 * dllat[3] * cott^2) * saz^2
dgy[3] <- 4.0 * dllat[3] * cott * csts / st
dgx[4] <-
(dhlat[4] + dllat[4] * (33.0 - 45.0 * ct2) / (5.0 * ct2 - 3.0)) * caz^2 +
(dhlat[4] - dllat[4] * (1.00 + 10.0 * ct2 / (5.0 * ct2 - 3.0))) * saz^2
dgy[4] <- 0.0
dgx[5] <-
(dhlat[5] - dllat[5] * (1.0 + 10.0 * (1.0 - 4.0 * ct2) / (1.0 - 5.0 * ct2))) * caz^2 +
(dhlat[5] + dllat[5] * (cott^2 - 1.0 / st2 - 10.0 * ct2 / (5.0 * ct2 - 1.0))) * saz^2
dgy[5] <- -20.0 * dllat[5] * ct * csts / (5.0 * ct2 - 1.0)
dgx[6] <- (dhlat[6] + dllat[6] * (2.0 * cott^2 - 7.0)) * caz^2 +
(dhlat[6] + dllat[6] * (2.0 * cott^2 - 1.0 - 4.0 / st2)) * saz^2
dgy[6] <- -4.0 * dllat[6] * (cott - 1.0 / cott) * csts / st
dgx[7] <- (dhlat[7] + dllat[7] * (6.0 * cott^2 - 3.0)) * caz^2 +
(dhlat[7] + dllat[7] * (3.0 * cott^2 - 9.0 / st2)) * saz^2
dgy[7] <- 12.0 * dllat[7] * cott * csts / st
dgx[8] <-
(dhlat[8] - 4.0 * dllat[8] * (4.0 - 3.0 * (5.0 * ct2 - 1.0) / (35.0 * ct2 * ct2 - 30.0 * ct2 + 3.0))) * caz^2 +
(dhlat[8] - 4.0 * dllat[8] * (1.0 + 3.0 * (5.0 * ct2 - 1.0) / (35.0 * ct2 * ct2 - 30.0 * ct2 + 3.0))) * saz^2
dgy[8] <- 0.0
dgx[9] <-
(dhlat[9] - 2.0 * dllat[9] * (8.0 - 3.0 / (7.0 * ct2 - 3.0))) * caz^2 +
(dhlat[9] - 2.0 * dllat[9] * (2.0 + 3.0 / (7.0 * ct2 - 3.0))) * saz^2
dgy[9] <- dllat[9] * 3.0 / ct * (1.0 + 2.0 / (7.0 * ct2 - 3.0)) * saz2
dgx[10] <-
(dhlat[10] - 4.0 * dllat[10] * (4.0 + 3.0 * ct2 / (7.0 * ct2^2 - 8.0 * ct2 + 1.0))) * caz^2 +
(dhlat[10] - 4.0 * dllat[10] * (1.0 - 3.0 * ct2 / (7.0 * ct2^2 - 8.0 * ct2 + 1.0))) * saz^2
dgy[10] <-
-dllat[10] * 6.0 * ct / st^2 * (1.0 - 4.0 / (7.0 * ct2 - 1.0)) * saz2
dgx[11] <- (dhlat[11] - 2.0 * dllat[11] * (8.0 - 3.0 / st2)) * caz^2 +
(dhlat[11] - 2.0 * dllat[11] * (2.0 + 3.0 / st2)) * saz^2
dgy[11] <- dllat[11] * 3.0 / ct * (3.0 - 2.0 / st2) * saz2
dgx[12] <- (dhlat[12] - 4.0 * dllat[12] * (4.0 - 3.0 / st2)) * caz^2 +
(dhlat[12] - 4.0 * dllat[12] * (1.0 + 3.0 / st2)) * saz^2
dgy[12] <- dllat[12] * 12.0 * ct / st2 * saz2
self$station$dgx <- dgx
self$station$dgy <- dgy
self$station$dgk[1:12] <-
self$station$dgk[1:12] * sqrt(dgx[1:12]^2 + dgy[1:12]^2) * dfak
self$pk[1:12] <- self$pk[1:12] + atan2(dgy[1:12], dgx[1:12]) * to_radians
},
ocean_tides = function() {
dfak <- 1e3 / self$station$gravity
self$station$dgk <- self$station$dgk * dfak
self$pk[] <- 0.0
},
predict = function(method = "gravity",
return_matrix = FALSE,
n_thread = 1) {
self$apply_method(method)
if (return_matrix) {
mat <- self$calculate(predict = TRUE, n_thread = n_thread)
colnames(mat) <- method
return(mat)
} else {
self$tides[[method]] <-
as.numeric(self$calculate(
predict = TRUE, n_thread = n_thread
))
}
# reset parameters after calculation
self$prepare_station(
self$station$latitude,
self$station$longitude,
self$station$elevation,
self$station$azimuth,
self$station$gravity,
self$station$earth_radius,
self$station$earth_eccen
)
invisible(self)
},
analyze = function(method = "gravity",
return_matrix = FALSE,
scale = TRUE,
n_thread = 1) {
self$apply_method(method)
mat <- self$calculate(predict = FALSE, scale = scale, n_thread = n_thread)
# reset parameters after calculation
self$prepare_station(
self$station$latitude,
self$station$longitude,
self$station$elevation,
self$station$azimuth,
self$station$gravity,
self$station$earth_radius,
self$station$earth_eccen
)
if (return_matrix) {
colnames(mat) <- self$catalog$col_names
return(mat)
}
self$tides[self$catalog$col_names] <- mat
invisible(self)
},
apply_method = function(method) {
if (method == "gravity") {
self$gravity()
} else if (method == "ocean_tides") {
self$ocean_tides()
} else if (method == "tidal_potential") {
self$tidal_potential()
} else if (method == "tidal_tilt") {
self$tidal_tilt()
} else if (method == "vertical_displacement") {
self$vertical_displacement()
} else if (method == "horizontal_displacement") {
self$horizontal_displacement()
} else if (method == "n_s_displacement") {
self$n_s_displacement()
} else if (method == "e_w_displacement") {
self$e_w_displacement()
} else if (method == "vertical_strain") {
self$vertical_strain()
} else if (method == "areal_strain") {
self$areal_strain()
} else if (method == "horizontal_strain") {
self$horizontal_strain()
} else if (method == "volume_strain") {
self$volume_strain()
}
},
calculate = function(predict = TRUE, scale = TRUE, n_thread = 1) {
n_cores <- detectCores()
if(n_thread > n_cores) {
warning("number of threads is more than the number of cores, setting
n_thread equal to the maximum number of cores.")
n_thread <- n_cores
}
et_calculate(
self$astro$astro,
self$astro$astro_der,
self$catalog$k,
self$pk,
self$delta,
self$deltar,
self$catalog$cc,
self$catalog$ss,
self$station$dgk,
self$catalog$jcof - 1,
self$datetime$j2000,
self$love_params$dom0,
self$love_params$domr,
self$catalog$id,
self$catalog$wave_groups$multiplier,
predict,
scale,
n_thread
)
},
pole_tide = function() {
self$tides$pole_tide <- self$pole_t
invisible(self)
},
lod_tide = function() {
self$tides$lod_tide <- self$lod_t
invisible(self)
},
tide = function() {
invisible(self$tides)
},
print = function(...) {
cat("Earthtide: \n")
cat(" Times: ", length(self$datetime$utc), "\n", sep = "")
cat(" range (utc): ", as.character(head(self$datetime$utc, 1)), " to ",
as.character(tail(self$datetime$utc, 1)), "\n",
sep = ""
)
cat(" dt[1] (sec): ", as.numeric(self$datetime$utc[2]) -
as.numeric(self$datetime$utc[1]), "\n", sep = "")
cat(" Station: \n")
cat(" latitude: ", head(self$station$latitude), "\n", sep = "")
cat(" longitude: ", head(self$station$longitude), "\n", sep = "")
cat(" elevation: ", head(self$station$elevation), "\n", sep = "")
cat(" gravity: ", head(self$station$gravity), "\n", sep = "")
cat(" azimuth: ", head(self$station$azimuth), "\n", sep = "")
cat(" Wave groups: \n")
cat(" catalog: ", head(self$catalog$catalog), "\n", sep = "")
cat(" cutoff: ", head(self$catalog$cutoff), "\n", sep = "")
cat(" n waves: ", head(self$catalog$n_constituents), "\n", sep = "")
cat(" n groups: ", nrow(self$catalog$wave_groups), "\n", sep = "")
invisible(self)
},
# time variables
datetime = list(),
# astro arguments
astro = list(),
# wavegroup catalog
catalog = list(),
# geodetic coeficients
station = list(),
# love and shida numbers
love_params = list(),
pole_quotient = 1.16, # 1.1788, # Ducarme et al. 2006
pk = rep(0, 25),
delta = rep(1.0, 25),
deltar = 0.0,
# outputs
lod_t = NA_real_,
pole_t = NA_real_,
tides = list()
)
)
Try the earthtide package in your browser
Any scripts or data that you put into this service are public.
earthtide documentation built on Nov. 16, 2023, 5:07 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.
#' @title Earthtide class
#'
#' @description Class to generate synthetic earthtide signals.
#'
#' @rdname Earthtide_class
#'
#' @section Usage:
#' \preformatted{
#' et <- Earthtide$new(
#' utc = as.POSIXct("2017-01-01", tz = "UTC") + 0:(24 * 7) * 3600,
#' latitude = 52.3868,
#' longitude = 9.7144,
#' catalog = "ksm04",
#' wave_groups = data.frame(start = 0.0, end = 6.0))
#'
#' et$predict(method = "gravity", n_thread = 1)
#' et$analyze(method = "gravity", n_thread = 1)
#' et$lod_tide()
#' et$pole_tide()
#' et$tide()
#' et$print()
#' }
#'
#' @section Arguments:
#' \code{Earthtide$new}
#' \describe{
#' \item{et: }{An \code{Earthtide} object.}
#' \item{utc: }{The date-time in UTC (POSIXct vector).}
#' \item{latitude: }{The station latitude (WGS84) (degree) (numeric)
#' defaults to 0.0}
#' \item{longitude: }{The station longitude (WGS84) (degree) (numeric)
#' defaults to 0.0}
#' \item{elevation: }{The station ellipsoidal height (WGS84) (m) (numeric)
#' defaults to 0.0}
#' \item{azimuth: }{Earth azimuth (numeric) defaults to 0 (degrees)}
#' \item{gravity: }{Gravity at the station (m/s^2) (numeric) 0 to
#' estimate gravity from elevation and latitude.}
#' \item{earth_radius: }{Radius of earth (m) (numeric) defaults to 6378136.3 }
#' \item{earth_eccen: }{Eccentricity of earth (numeric)
#' defaults to 6.69439795140e-3}
#' \item{cutoff: }{Cutoff amplitude for constituents (numeric)
#' defaults to 1e-6}
#' \item{wave_groups: }{Two column data.frame having start and end of
#' frequency groups (data.frame). This data.frame must have two columns
#' with the names 'start', and 'end' signifying the start and end of the
#' wave groupings. An optional third column 'multiplier' can be provided
#' to scale the particular wave group. If column names do no match, the
#' inferred column positions are start, end, multiplier.}
#' \item{catalog: }{Use the "hw95s" catalog or "ksm04" catalog (character).}
#' \item{eop: }{User defined Earth Orientation Parameter (EOP) data.frame
#' with the following columns: datetime, ddt, ut1_utc, lod, x, y, dx, dy}
#' \item{...: }{Currently not used.}
#' }
#'
#' \code{Earthtide$predict, Earthtide$analyze}
#' \describe{
#' \item{method: }{For \code{predict} and \code{analyze}. One of "gravity",
#' "tidal_potential", "tidal_tilt", "vertical_displacement",
#' "horizontal_displacement", "n_s_displacement", "e_w_displacement",
#' "vertical_strain", "areal_strain", "volume_strain", "horizontal_strain"
#' or "ocean_tides".}
#' \item{return_matrix: }{For \code{predict} and \code{analyze}. Return a
#' matrix of tidal values instead of data.frame. The datetime column will
#' not be present in this case (logical).}
#' \item{n_thread: }{For \code{predict} and \code{analyze}. Number of threads
#' to use for parallel processing.}
#' }
#'
#' @section Details:
#'
#' \code{$new(utc, latitude, longitude, elevation, azimuth, gravity,} \cr
#' \code{earth_radius, earth_eccen, cutoff, wave_groups, catalog, ...)} \cr
#' create a new \code{Earthtide} object and initialize catalog, station
#' and times.
#'
#' \code{$predict(method, astro_argument, return_matrix)} generate a combined
#' synthetic Earth tide.
#'
#' \code{$analyze(method, astro_argument, return_matrix, scale)} generate
#' components of the Earth tide for analysis.
#'
#' \code{$lod_tide()} generate components of the LOD (Length Of Day) tide.
#'
#' \code{$pole_tide()} generate components of the pole tide.
#'
#' \code{$tide()} get the tide \code{data.frame}.
#'
#' \code{$print()} print the \code{Earthtide} object.
#'
#' @docType class
#' @aliases Earthtide-class
#' @importFrom R6 R6Class
#' @importFrom RcppThread detectCores
#' @importFrom stats approx
#' @importFrom utils read.table
#' @importFrom utils read.fwf
#' @importFrom utils download.file
#' @importFrom utils data
#' @format An \code{\link{R6Class}} generator object
#'
#'
#' @references Hartmann, T., Wenzel, H.-G., 1995. The HW95 tidal potential catalogue. Geophys. Res. Lett. 22, 3553-3556. \doi{10.1029/95GL03324}
#' @references Kudryavtsev, S.M., 2004. Improved harmonic development of the Earth tide-generating potential. J. Geod. 77, 829-838. \doi{10.1007/s00190-003-0361-2}
#' @references Wenzel, H.G., 1996. The nanogal software: Earth tide data processing package ETERNA 3.30. Bull. Inf. Marées Terrestres, 124, pp.9425-9439. \url{https://www.eas.slu.edu/GGP/ETERNA34/MANUAL/ETERNA33.HTM}
#'
#' @examples
#'
#' et <- Earthtide$new(
#' utc = as.POSIXct("2017-01-01", tz = "UTC") + 0:(24 * 7) * 3600,
#' latitude = 52.3868,
#' longitude = 9.7144,
#' catalog = "ksm04",
#' wave_groups = data.frame(start = 0.0, end = 6.0)
#' )
#'
#' et$predict(method = "gravity")
#'
#' plot(gravity ~ datetime, et$tide(), type = "l")
#'
#' @name Earthtide
NULL
#' @export
Earthtide <- R6Class(
"et",
public = list(
# initialization
initialize = function(utc,
latitude = 0,
longitude = 0,
elevation = 0,
azimuth = 0,
gravity = 0,
earth_radius = 6378136.3,
earth_eccen = 6.69439795140e-3,
cutoff = 1e-6,
wave_groups = NULL,
catalog = "ksm04",
eop = NULL,
...) {
# Initialize class using input values
self$prepare_datetime(utc, eop)
self$prepare_station(
latitude, longitude, elevation, azimuth, gravity,
earth_radius, earth_eccen
)
self$prepare_astro()
self$prepare_catalog(cutoff, wave_groups, catalog)
self$love_params <- love(latitude, elevation)
invisible(self)
},
# Initialize class using input values
prepare_datetime = function(utc, eop) {
self$datetime <- .prepare_datetime(utc, eop)
self$tides <- data.frame(datetime = utc)
},
# Calculate the astronical arguments and derivative
prepare_astro = function() {
self$astro <- .prepare_astro(self)
},
# Subset values based using a cutoff amplitude and
# cutoff frequencies
prepare_catalog = function(cutoff, wave_groups = NULL, catalog = "ksm04") {
self$catalog <- .prepare_catalog(cutoff, wave_groups, catalog = catalog)
},
# Calculate properties based on the location of station
prepare_station = function(latitude, longitude, elevation, azimuth,
gravity, earth_radius, earth_eccen) {
self$station <- .prepare_station(
self, latitude, longitude, elevation,
azimuth, gravity, earth_radius,
earth_eccen
)
self$pole_quotient <- 1.16 # 1.1788, # Ducarme et al. 2006
self$pk <- rep(0, 25)
self$delta <- rep(1.0, 25)
self$deltar <- 0.0
},
gravity = function() {
self$station$dgk <- self$station$dgz
self$pk[] <- 180
self$delta[1:12] <- self$love_params$dglat
self$deltar <- self$love_params$dgr
},
tidal_potential = function() {
self$delta[1:12] <- self$love_params$dklat
self$deltar <- self$love_params$dkr
},
tidal_tilt = function() {
to_radians <- pi / 180
cos_azimuth <- cos(self$station$azimuth * to_radians)
sin_azimuth <- sin(self$station$azimuth * to_radians)
x_comp <- self$station$dgx[1:12] * cos_azimuth
y_comp <- self$station$dgy[1:12] * sin_azimuth
self$station$dgk[1:12] <- sqrt((x_comp)^2 + (y_comp)^2) * self$station$df
wh <- which(x_comp != 0 | y_comp != 0)
self$pk[] <- 0.0
self$pk[wh] <- 180 / pi * atan2(y_comp[wh], x_comp[wh])
# from etpots
self$delta[1:12] <- self$love_params$dtlat
self$deltar <- self$love_params$dkr - self$love_params$dhr
},
vertical_displacement = function() {
dfak <- 1e3 / self$station$gravity
self$station$dgk[1:12] <- self$station$dgk[1:12] *
self$love_params$dhlat[1:12] * dfak
self$pk[] <- 0.0
},
horizontal_displacement = function() {
to_radians <- pi / 180
dfak <- 1e-6 * self$station$geo_radius / self$station$gravity
cos_azimuth <- cos(self$station$azimuth * to_radians)
sin_azimuth <- sin(self$station$azimuth * to_radians)
x_comp <- self$station$dgx[1:12] * cos_azimuth
y_comp <- self$station$dgy[1:12] * sin_azimuth
# # shida numbers
# dl0 <- c(
# 0.08400, 0.08410, 0.08520, 0.01490, 0.01490, 0.01490,
# 0.01490, 0.01000, 0.01000, 0.01000, 0.01000, 0.01000
# )
#
# if(any(abs(self$love_params$dllat[1:12]-dl0) > 0.04)) {
# self$love_params$dllat[1:12] <- dl0
# }
self$station$dgk[1:12] <- sqrt((x_comp)^2 + (y_comp)^2) *
self$love_params$dllat[1:12] * dfak
wh <- which(x_comp != 0 | y_comp != 0)
self$pk[] <- 0.0
self$pk[wh] <- 180 / pi * atan2(y_comp[wh], x_comp[wh])
# from etpots (this is a modification as I think they are 0 and 1 in ETERNA)
# look into this more closely
self$delta[1:12] <- 1.0 #self$love_params$dtlat
self$deltar <- 0.0 #self$love_params$dkr - self$love_params$dhr
},
n_s_displacement = function() {
to_radians <- pi / 180
dfak <- 1e-6 * self$station$geo_radius / self$station$gravity
cos_azimuth <- cos(0 * to_radians)
sin_azimuth <- sin(0 * to_radians)
x_comp <- self$station$dgx[1:12] * cos_azimuth
y_comp <- self$station$dgy[1:12] * sin_azimuth
self$station$dgk[1:12] <- sqrt((x_comp)^2 + (y_comp)^2) *
self$love_params$dllat[1:12] * dfak
self$pk[] <- 0.0
self$pk[1:12] <- 180 / pi * atan2(y_comp, x_comp)
wh <- which(x_comp == 0 & y_comp == 0)
self$pk[wh] <- 0.0
},
e_w_displacement = function() {
to_radians <- pi / 180
dfak <- 1e-6 * self$station$geo_radius / self$station$gravity
cos_azimuth <- cos(90 * to_radians)
sin_azimuth <- sin(90 * to_radians)
x_comp <- self$station$dgx[1:12] * cos_azimuth
y_comp <- self$station$dgy[1:12] * sin_azimuth
self$station$dgk[1:12] <- sqrt((x_comp)^2 + (y_comp)^2) *
self$love_params$dllat[1:12] * dfak
self$pk[] <- 0.0
self$pk[1:12] <- 180 / pi * atan2(y_comp, x_comp)
wh <- which(x_comp == 0 & y_comp == 0)
self$pk[wh] <- 0.0
},
vertical_strain = function(poisson = 0.25) {
dfak <- 1.e9 * poisson / (poisson - 1.0)
self$strain(dfak)
},
areal_strain = function() {
dfak <- 1.e9
self$strain(dfak)
},
volume_strain = function(poisson = 0.25) {
dfak <- 1.e9 * (1.0 - 2.0 * poisson) / (1.0 - poisson)
self$strain(dfak)
},
strain = function(dfak) {
scale <- c(rep(-6, 3), rep(-12, 4), rep(-20, 5))
self$station$dgk[1:12] <- self$station$dgk[1:12] * dfak *
(2.0 * self$love_params$dhlat[1:12] +
scale * self$love_params$dllat[1:12]) /
(self$station$gravity * self$station$geo_radius)
},
horizontal_strain = function() {
to_radians <- pi / 180
dgx <- self$station$dgx
dgy <- self$station$dgy
dhlat <- self$love_params$dhlat
dllat <- self$love_params$dllat
theta <- self$station$theta * to_radians
azr <- (self$station$azimuth + 180) * to_radians
caz <- cos(azr)
saz <- sin(azr)
saz2 <- sin(2.0 * azr)
csts <- -0.5 * sin(2.0 * azr)
ct <- sin(self$station$geo_latitude * to_radians)
st <- cos(self$station$geo_latitude * to_radians)
ct2 <- ct * ct
st2 <- st * st
cc2 <- cos(2.0 * self$station$geo_latitude * to_radians)
c2t <- -cc2
cott <- 1.0 / tan(theta)
cott2 <- 1.0 / tan(2.0 * theta)
dfak <- 1e9 / (self$station$geo_radius * self$station$gravity)
dgx[1] <-
(dhlat[1] - (6.0 * dllat[1] * c2t) / (3.0 * ct2 - 1.0)) * caz^2 +
(dhlat[1] - (6.0 * dllat[1] * ct2) / (3.0 * ct2 - 1.0)) * saz^2
dgy[1] <- 0.0
dgx[2] <- (dhlat[2] - 4.0 * dllat[2]) * caz^2 +
(dhlat[2] - dllat[2] / st2 + 2.0 * dllat[2] * cott * cott2) * saz^2
dgy[2] <- 2.0 * dllat[2] * (2.0 * cott2 - cott) * csts / st
dgx[3] <- (dhlat[3] + 2.0 * dllat[3] * (cott^2 - 1.0)) * caz^2 +
(dhlat[3] - 4.0 * dllat[3] / st2 + 2.0 * dllat[3] * cott^2) * saz^2
dgy[3] <- 4.0 * dllat[3] * cott * csts / st
dgx[4] <-
(dhlat[4] + dllat[4] * (33.0 - 45.0 * ct2) / (5.0 * ct2 - 3.0)) * caz^2 +
(dhlat[4] - dllat[4] * (1.00 + 10.0 * ct2 / (5.0 * ct2 - 3.0))) * saz^2
dgy[4] <- 0.0
dgx[5] <-
(dhlat[5] - dllat[5] * (1.0 + 10.0 * (1.0 - 4.0 * ct2) / (1.0 - 5.0 * ct2))) * caz^2 +
(dhlat[5] + dllat[5] * (cott^2 - 1.0 / st2 - 10.0 * ct2 / (5.0 * ct2 - 1.0))) * saz^2
dgy[5] <- -20.0 * dllat[5] * ct * csts / (5.0 * ct2 - 1.0)
dgx[6] <- (dhlat[6] + dllat[6] * (2.0 * cott^2 - 7.0)) * caz^2 +
(dhlat[6] + dllat[6] * (2.0 * cott^2 - 1.0 - 4.0 / st2)) * saz^2
dgy[6] <- -4.0 * dllat[6] * (cott - 1.0 / cott) * csts / st
dgx[7] <- (dhlat[7] + dllat[7] * (6.0 * cott^2 - 3.0)) * caz^2 +
(dhlat[7] + dllat[7] * (3.0 * cott^2 - 9.0 / st2)) * saz^2
dgy[7] <- 12.0 * dllat[7] * cott * csts / st
dgx[8] <-
(dhlat[8] - 4.0 * dllat[8] * (4.0 - 3.0 * (5.0 * ct2 - 1.0) / (35.0 * ct2 * ct2 - 30.0 * ct2 + 3.0))) * caz^2 +
(dhlat[8] - 4.0 * dllat[8] * (1.0 + 3.0 * (5.0 * ct2 - 1.0) / (35.0 * ct2 * ct2 - 30.0 * ct2 + 3.0))) * saz^2
dgy[8] <- 0.0
dgx[9] <-
(dhlat[9] - 2.0 * dllat[9] * (8.0 - 3.0 / (7.0 * ct2 - 3.0))) * caz^2 +
(dhlat[9] - 2.0 * dllat[9] * (2.0 + 3.0 / (7.0 * ct2 - 3.0))) * saz^2
dgy[9] <- dllat[9] * 3.0 / ct * (1.0 + 2.0 / (7.0 * ct2 - 3.0)) * saz2
dgx[10] <-
(dhlat[10] - 4.0 * dllat[10] * (4.0 + 3.0 * ct2 / (7.0 * ct2^2 - 8.0 * ct2 + 1.0))) * caz^2 +
(dhlat[10] - 4.0 * dllat[10] * (1.0 - 3.0 * ct2 / (7.0 * ct2^2 - 8.0 * ct2 + 1.0))) * saz^2
dgy[10] <-
-dllat[10] * 6.0 * ct / st^2 * (1.0 - 4.0 / (7.0 * ct2 - 1.0)) * saz2
dgx[11] <- (dhlat[11] - 2.0 * dllat[11] * (8.0 - 3.0 / st2)) * caz^2 +
(dhlat[11] - 2.0 * dllat[11] * (2.0 + 3.0 / st2)) * saz^2
dgy[11] <- dllat[11] * 3.0 / ct * (3.0 - 2.0 / st2) * saz2
dgx[12] <- (dhlat[12] - 4.0 * dllat[12] * (4.0 - 3.0 / st2)) * caz^2 +
(dhlat[12] - 4.0 * dllat[12] * (1.0 + 3.0 / st2)) * saz^2
dgy[12] <- dllat[12] * 12.0 * ct / st2 * saz2
self$station$dgx <- dgx
self$station$dgy <- dgy
self$station$dgk[1:12] <-
self$station$dgk[1:12] * sqrt(dgx[1:12]^2 + dgy[1:12]^2) * dfak
self$pk[1:12] <- self$pk[1:12] + atan2(dgy[1:12], dgx[1:12]) * to_radians
},
ocean_tides = function() {
dfak <- 1e3 / self$station$gravity
self$station$dgk <- self$station$dgk * dfak
self$pk[] <- 0.0
},
predict = function(method = "gravity",
return_matrix = FALSE,
n_thread = 1) {
self$apply_method(method)
if (return_matrix) {
mat <- self$calculate(predict = TRUE, n_thread = n_thread)
colnames(mat) <- method
return(mat)
} else {
self$tides[[method]] <-
as.numeric(self$calculate(
predict = TRUE, n_thread = n_thread
))
}
# reset parameters after calculation
self$prepare_station(
self$station$latitude,
self$station$longitude,
self$station$elevation,
self$station$azimuth,
self$station$gravity,
self$station$earth_radius,
self$station$earth_eccen
)
invisible(self)
},
analyze = function(method = "gravity",
return_matrix = FALSE,
scale = TRUE,
n_thread = 1) {
self$apply_method(method)
mat <- self$calculate(predict = FALSE, scale = scale, n_thread = n_thread)
# reset parameters after calculation
self$prepare_station(
self$station$latitude,
self$station$longitude,
self$station$elevation,
self$station$azimuth,
self$station$gravity,
self$station$earth_radius,
self$station$earth_eccen
)
if (return_matrix) {
colnames(mat) <- self$catalog$col_names
return(mat)
}
self$tides[self$catalog$col_names] <- mat
invisible(self)
},
apply_method = function(method) {
if (method == "gravity") {
self$gravity()
} else if (method == "ocean_tides") {
self$ocean_tides()
} else if (method == "tidal_potential") {
self$tidal_potential()
} else if (method == "tidal_tilt") {
self$tidal_tilt()
} else if (method == "vertical_displacement") {
self$vertical_displacement()
} else if (method == "horizontal_displacement") {
self$horizontal_displacement()
} else if (method == "n_s_displacement") {
self$n_s_displacement()
} else if (method == "e_w_displacement") {
self$e_w_displacement()
} else if (method == "vertical_strain") {
self$vertical_strain()
} else if (method == "areal_strain") {
self$areal_strain()
} else if (method == "horizontal_strain") {
self$horizontal_strain()
} else if (method == "volume_strain") {
self$volume_strain()
}
},
calculate = function(predict = TRUE, scale = TRUE, n_thread = 1) {
n_cores <- detectCores()
if(n_thread > n_cores) {
warning("number of threads is more than the number of cores, setting
n_thread equal to the maximum number of cores.")
n_thread <- n_cores
}
et_calculate(
self$astro$astro,
self$astro$astro_der,
self$catalog$k,
self$pk,
self$delta,
self$deltar,
self$catalog$cc,
self$catalog$ss,
self$station$dgk,
self$catalog$jcof - 1,
self$datetime$j2000,
self$love_params$dom0,
self$love_params$domr,
self$catalog$id,
self$catalog$wave_groups$multiplier,
predict,
scale,
n_thread
)
},
pole_tide = function() {
self$tides$pole_tide <- self$pole_t
invisible(self)
},
lod_tide = function() {
self$tides$lod_tide <- self$lod_t
invisible(self)
},
tide = function() {
invisible(self$tides)
},
print = function(...) {
cat("Earthtide: \n")
cat(" Times: ", length(self$datetime$utc), "\n", sep = "")
cat(" range (utc): ", as.character(head(self$datetime$utc, 1)), " to ",
as.character(tail(self$datetime$utc, 1)), "\n",
sep = ""
)
cat(" dt[1] (sec): ", as.numeric(self$datetime$utc[2]) -
as.numeric(self$datetime$utc[1]), "\n", sep = "")
cat(" Station: \n")
cat(" latitude: ", head(self$station$latitude), "\n", sep = "")
cat(" longitude: ", head(self$station$longitude), "\n", sep = "")
cat(" elevation: ", head(self$station$elevation), "\n", sep = "")
cat(" gravity: ", head(self$station$gravity), "\n", sep = "")
cat(" azimuth: ", head(self$station$azimuth), "\n", sep = "")
cat(" Wave groups: \n")
cat(" catalog: ", head(self$catalog$catalog), "\n", sep = "")
cat(" cutoff: ", head(self$catalog$cutoff), "\n", sep = "")
cat(" n waves: ", head(self$catalog$n_constituents), "\n", sep = "")
cat(" n groups: ", nrow(self$catalog$wave_groups), "\n", sep = "")
invisible(self)
},
# time variables
datetime = list(),
# astro arguments
astro = list(),
# wavegroup catalog
catalog = list(),
# geodetic coeficients
station = list(),
# love and shida numbers
love_params = list(),
pole_quotient = 1.16, # 1.1788, # Ducarme et al. 2006
pk = rep(0, 25),
delta = rep(1.0, 25),
deltar = 0.0,
# outputs
lod_t = NA_real_,
pole_t = NA_real_,
tides = list()
)
)
Try the earthtide package in your browser
Any scripts or data that you put into this service are public.
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.