--- 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.