---
title: "How moving-window downscaling works"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{How moving-window downscaling works}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r libraries, include = FALSE}
library(topocast)
library(terra)
```

Downscaling a climate grid with terrain rests on a long-standing idea: a climate
variable often tracks elevation through a roughly linear relationship, so a
high-resolution elevation model can carry the variable to a finer grid through a
regression. This is the regression-on-elevation step behind high-resolution
climate surfaces such as CHELSA (Karger et al. 2017). `topocast` generalizes the
predictor from elevation alone to any set of aligned terrain covariates, fits the
relationship locally rather than globally, and keeps the per-cell cost
independent of the neighborhood size.

This vignette describes the method: the local regression, the summed-area tables
that make it fast, the numerical choices, and what the output does and does not
preserve.

## A regression in every neighborhood

A global regression of the coarse response $y$ on $k$ predictors $x_1, \dots,
x_k$ assumes one relationship over the whole map. Terrain-climate relationships
are not that stable: a lapse rate on a windward slope differs from one in a rain
shadow. `topocast` instead fits a separate ordinary least squares regression in a
square window around every coarse cell. For the cell at row $r$, column $c$, the
window of radius $\rho$ holds the cells within $\rho$ rows and columns, and the
fit solves

$$
\min_{\beta} \sum_{(i,j) \in W_{rc}} \Big( y_{ij} - \beta_0 - \sum_{m=1}^{k} \beta_m \, x_{m,ij} \Big)^2 .
$$

The result is a grid of intercepts $\hat\beta_0(r,c)$ and a grid of slopes
$\hat\beta_m(r,c)$ for each predictor, all at the coarse resolution. These
coefficient grids are the locally varying relationship between the response and
the terrain, and `topocast(..., coefficients = TRUE)` returns them on the fine
grid alongside the fitted layer.

The downscaling step resamples each coefficient grid to the fine resolution and
evaluates the regression on the fine predictors:

$$
\hat y(f) = \hat\beta_0(f) + \sum_{m=1}^{k} \hat\beta_m(f)\, x_m(f),
$$

where $f$ indexes fine cells and the coefficients are the resampled coarse grids.
The fine-scale structure of the output comes entirely from the fine predictors;
the coefficients carry the relationship.

## Constant cost with summed-area tables

A naive moving window re-sums every cell in the neighborhood for every center, so
its cost grows with the window area. For a radius of 20 that is over 1,600 cells
touched per output cell. `topocast` avoids this with summed-area tables, also
called integral images.

The normal equations for the window fit are $A \beta = b$, where $A$ holds the
within-window sums $\sum 1$, $\sum x_m$, and $\sum x_m x_{m'}$, and $b$ holds
$\sum y$ and $\sum x_m y$. Each of these is a sum of a fixed quantity over the
window. For any grid $g$, the summed-area table $S$ stores the sum over the
rectangle from the top-left corner to each cell:

$$
S(r,c) = \sum_{i \le r,\, j \le c} g_{ij}.
$$

With $S$ built once in a single pass, the sum over any rectangular window is four
lookups:

$$
\sum_{(i,j)\in W} g_{ij} = S(r_1,c_1) - S(r_0-1,c_1) - S(r_1,c_0-1) + S(r_0-1,c_0-1).
$$

`topocast` builds one table for the cell count, one for $\sum y$, one per
predictor for $\sum x_m$ and $\sum x_m y$, and one per predictor pair for $\sum
x_m x_{m'}$. Every window then assembles $A$ and $b$ from four lookups each,
solves a $(k{+}1) \times (k{+}1)$ system, and stores the coefficients. The total
cost is proportional to the number of cells and the small per-cell solve, with no
dependence on the radius.

The radius below changes the smoothness of the slope grid but not the runtime.

```{r cost, fig.height = 3.2}
set.seed(1)
g <- rast(nrows = 80, ncols = 80, xmin = 0, xmax = 8000, ymin = 0, ymax = 8000,
          crs = "EPSG:32632")
elev <- setValues(g, 1000 + 500 * sin(crds(g)[, 1] / 1500))
prec <- 800 - 0.15 * elev + setValues(g, rnorm(ncell(g), 0, 15))

y <- as.matrix(prec, wide = TRUE)
x <- as.matrix(elev, wide = TRUE)
small <- window_regression(y, x, radius = 3)
large <- window_regression(y, x, radius = 15)

op <- par(mfrow = c(1, 2), mar = c(2, 2, 2, 4))
plot(setValues(elev, as.vector(t(small$slope[[1]]))), main = "slope, radius 3")
plot(setValues(elev, as.vector(t(large$slope[[1]]))), main = "slope, radius 15")
par(op)
```

## Numerical and masking choices

**Global centering.** Before accumulation the response and each predictor are
centered on their global mean over the valid cells. Centering keeps the entries of
$A$ small even when elevation runs into the thousands, which conditions the
solve. The intercept is mapped back to the raw response and predictor scale on
output, so the returned coefficients are on the original units.

**Complete-case masking.** A cell contributes to a window only where the response
and every predictor are finite. A separate count table tracks how many valid
cells fall in each window, and invalid cells are set to zero in the centered
copies so they drop out of the sums without distorting them.

**Degenerate windows.** A cell is returned as `NA` in three cases: the window
holds fewer valid cells than the model needs (`k + 1`, raised by `min_cells`); a
predictor has within-window variance below `min_variance`, so its slope is not
identified; or the assembled system is singular. These guards keep a flat or
data-starved neighborhood from producing a meaningless coefficient.

## Reporting fit quality

The same sufficient statistics give the local quality of each fit at no extra
window cost. One more summed-area table, of the squared centered response, yields
the within-window total sum of squares $\mathrm{SST} = \sum (y - \bar y_W)^2$,
and the residual sum of squares follows from the fitted coefficients as
$\mathrm{SSE} = y'y - \hat\beta' b$, where $b$ is the response side of the normal
equations. Their ratio is the coefficient of determination,

$$
R^2(r,c) = 1 - \frac{\mathrm{SSE}}{\mathrm{SST}},
$$

returned by `diagnostics = TRUE` as an `r.squared` grid on the fine resolution. It
maps where the terrain explains the response and where the downscaled field
carries little more than the coarse level.

The same statistics also give the residual standard deviation,
$\hat\sigma(r,c) = \sqrt{\mathrm{SSE} / (n - k - 1)}$, in the response's own units
rather than R-squared's unitless, variance-relative scale, and the window's valid
cell count $n$ itself. `diagnostics = TRUE` returns both, as `residual.sd` and
`n.valid`. A low R-squared is ambiguous on its own: it can mean the terrain-response
relationship is genuinely weak in that neighborhood, or that the window barely
cleared the minimum valid-cell count and the fit is just noisy from too little data.
`n.valid` tells the two apart. Because the valid-cell mask is complete-case across
every response fit in one call, `n.valid` is one shared grid rather than one per
response.

## Carrying a time series

A climatology is one field. A time series shares the terrain relationship across
many periods, so refitting the regression each period would be wasteful and would
let the spatial pattern drift between steps. `topocast` instead downscales the
baseline once and carries each period's coarse anomaly onto the fine baseline.

For a period with coarse value $v$, coarse baseline $b$, and fine baseline
$\hat y_b$, the ratio path returns

$$
\hat y = \hat y_b \cdot \frac{v}{b},
$$

and the additive path returns $\hat y = \hat y_b + (v - b)$. The ratio suits
non-negative variables such as precipitation; the additive form suits variables
such as temperature where a difference is the natural anomaly. Both carry the
coarse anomaly verbatim, so the downscaled series tracks the coarse values over
time even though the baseline itself is a regression surface. A zero or missing
baseline in the ratio path returns `NA` rather than dividing by zero.

## What the output preserves

The downscaled field is the regression evaluated on the fine predictors. Two
consequences are worth stating plainly.

The output is not mass conserving. Averaging the fine field back to the coarse
grid does not reproduce the coarse input, because the coarse residual, the part
of $y$ the regression did not explain, is discarded. For a smoothed climate
product this is usually acceptable, and the `anomaly` path restores fidelity to
the coarse values for the time-varying part.

The slopes are applied to the full fine-predictor range. Where the fine
elevation reaches values the coarse window never saw, the regression extrapolates
linearly. The fine output is only as trustworthy as a linear relationship over
that extended range, which is a reason to keep predictors that stay roughly
linear with the response and to inspect the result where the terrain is extreme.
Setting `clamp = TRUE` bounds the downscaled field to the observed range of the
coarse response, capping such extrapolation at the values the variable actually
takes.

## Reference

Karger, D. N., Conrad, O., Bohner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W.,
Zimmermann, N. E., Linder, H. P., and Kessler, M. (2017). Climatologies at high
resolution for the earth's land surface areas. *Scientific Data* 4, 170122.
<doi:10.1038/sdata.2017.122>
