Skip to contents

Model T_i(t) = theta(t) + b_i + W(x_i, t) (Rohde et al. 2013). Starting from station means, each iteration fits the weather correlation R(d) = c exp(-d / L) by weighted least squares to pairwise residual correlations (weights the overlap lengths; grid then golden-section search on L), estimates theta(t) by ordinary block kriging of the domain mean over a grid, re-estimates the baselines and centres theta. Identical to the Python arm morie.fn.berkeley.

Usage

BerkeleyEarth(
  stations,
  series,
  grid = NULL,
  spherical = FALSE,
  n_iter = 20L,
  tol = 1e-10
)

Arguments

stations

Matrix of coordinates (x, y), or (lat, lon) with spherical.

series

Matrix of station temperatures (stations by times, NA missing).

grid

Optional matrix of domain points (default 10 x 10 on the bounding box).

spherical

Great-circle kilometres from latitude and longitude.

n_iter, tol

Iteration cap and convergence tolerance.

Value

List: theta, theta_se, baselines, correlation_c, correlation_length, nugget, correlation_sse, iterations.

References

Rohde, R. et al. (2013). Berkeley Earth temperature averaging process. Geoinformatics and Geostatistics: An Overview 1(2).

Examples

st <- rbind(c(0, 0), c(1, 0), c(0, 1), c(1, 1), c(0.5, 0.5))
th <- c(0, 0.3, -0.2, 0.5, 0.1, 0.4)
T <- outer(10 + 2 * (0:4), th, "+") + 0.05 * sin(outer(1:5, 1:6))
BerkeleyEarth(st, T)$theta
#> [1] -0.17768920  0.12373335 -0.37557544  0.32375943 -0.08017146  0.18594332