## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4,
  dev = "svglite",
  fig.ext = "svg"
)

## ----libraries----------------------------------------------------------------
library(topocast)
library(terra)

## ----terrain------------------------------------------------------------------
set.seed(1)
fine <- rast(nrows = 120, ncols = 120,
             xmin = 0, xmax = 12000, ymin = 0, ymax = 12000,
             crs = "EPSG:32632")
xy <- crds(fine)
elev_fine <- setValues(fine, 1200 + 600 * sin(xy[, 1] / 2500) +
                               400 * cos(xy[, 2] / 2000))
names(elev_fine) <- "elev"
slope_fine <- terrain(elev_fine, v = "slope", unit = "degrees")
names(slope_fine) <- "slope"

## ----aggregate----------------------------------------------------------------
elev_coarse  <- aggregate(elev_fine,  fact = 6, fun = "mean")
slope_coarse <- aggregate(slope_fine, fact = 6, fun = "mean")

## ----climate------------------------------------------------------------------
prec_coarse <- 900 - 0.18 * elev_coarse + 4 * slope_coarse +
  setValues(elev_coarse, rnorm(ncell(elev_coarse), 0, 20))
names(prec_coarse) <- "prec"

## ----stacks-------------------------------------------------------------------
coarse  <- c(prec_coarse, elev_coarse, slope_coarse)
terrain <- c(elev_fine, slope_fine)
coarse

## ----single-------------------------------------------------------------------
fine_prec <- topocast(prec ~ elev, data = coarse, onto = terrain, radius = 4)
fine_prec

## ----plot-single, fig.height = 3.2--------------------------------------------
op <- par(mfrow = c(1, 2), mar = c(2, 2, 2, 4))
plot(prec_coarse, main = "coarse precipitation (600 m)")
plot(fine_prec,   main = "downscaled (100 m)")
par(op)

## ----multi--------------------------------------------------------------------
fine_prec2 <- topocast(prec ~ elev + slope, data = coarse, onto = terrain,
                       radius = 4)

## ----missing, error = TRUE----------------------------------------------------
try({
topocast(prec ~ elev + aspect, data = coarse, onto = terrain, radius = 4)
})

## ----several-responses--------------------------------------------------------
tmin_coarse <- 8 - 0.006 * elev_coarse +
  setValues(elev_coarse, rnorm(ncell(elev_coarse), 0, 0.5))
names(tmin_coarse) <- "tmin"
climate <- c(prec_coarse, tmin_coarse, elev_coarse, slope_coarse)

both <- topocast(cbind(prec, tmin) ~ elev + slope, data = climate, onto = terrain,
                 radius = 4)
both

## ----shortcut-----------------------------------------------------------------
shortcut <- topocast(prec ~ elev + slope, data = prec_coarse, onto = terrain,
                     radius = 4)

## ----shortcut-check-----------------------------------------------------------
round(max(abs(values(shortcut) - values(fine_prec2)), na.rm = TRUE), 4)

## ----coefficients, fig.height = 3.2-------------------------------------------
grids <- topocast(prec ~ elev, data = coarse, onto = terrain, radius = 4,
                  coefficients = TRUE)
names(grids)

op <- par(mfrow = c(1, 2), mar = c(2, 2, 2, 4))
plot(grids[["prec"]], main = "downscaled precipitation")
plot(grids[["elev"]], main = "local slope (mm per m)")
par(op)

## ----diagnostics, fig.height = 3.6--------------------------------------------
fit_quality <- topocast(prec ~ elev + slope, data = coarse, onto = terrain,
                        radius = 4, diagnostics = TRUE)
names(fit_quality)
plot(fit_quality[["r.squared"]], main = "local R-squared", range = c(0, 1))

## ----series-------------------------------------------------------------------
jan <- prec_coarse * 0.4
jul <- prec_coarse * 1.6
months <- c(jan, jul)
names(months) <- c("jan", "jul")

series <- topocast(prec ~ elev + slope, data = coarse, onto = terrain,
                   radius = 4, anomaly = months, type = "ratio")
nlyr(series)

## ----plot-series, fig.height = 3.2--------------------------------------------
op <- par(mfrow = c(1, 2), mar = c(2, 2, 2, 4))
plot(series[["jan"]], main = "January (100 m)")
plot(series[["jul"]], main = "July (100 m)")
par(op)

## ----sf-available, include = FALSE--------------------------------------------
has_sf <- requireNamespace("sf", quietly = TRUE)

## ----points, eval = has_sf----------------------------------------------------
library(sf)

plots <- st_as_sf(
  data.frame(x = c(2000, 4000, 6000, 8000, 10000, 5000),
             y = c(3000, 8000, 5000, 9000,  2000, 6000)),
  coords = c("x", "y"), crs = "EPSG:32632")

at_terrain <- terra::extract(terrain, vect(plots))
plots$elev  <- at_terrain$elev
plots$slope <- at_terrain$slope

at_plots <- topocast(prec ~ elev + slope, data = coarse, onto = plots, radius = 4)
at_plots

## ----plot-points, fig.height = 3.6, eval = has_sf-----------------------------
pv <- vect(plots)
plot(fine_prec2, main = "downscaled precipitation, with plot predictions")
plot(pv, add = TRUE, pch = 21, cex = 1.3)
text(crds(pv), labels = round(at_plots$prec), pos = 3, cex = 0.8)

## ----points-df, eval = has_sf-------------------------------------------------
topocast(prec ~ elev + slope, data = coarse, onto = plots, radius = 4,
         output = "data.frame")

## ----engine-------------------------------------------------------------------
y <- as.matrix(prec_coarse, wide = TRUE)
x <- list(as.matrix(elev_coarse, wide = TRUE),
          as.matrix(slope_coarse, wide = TRUE))
fit <- window_regression(y, x, radius = 4)
str(fit, max.level = 2)

## ----engine-recover-----------------------------------------------------------
big <- window_regression(y, x, radius = 19)
round(median(big$slope[[1]], na.rm = TRUE), 3)   # simulated -0.18