---
title: "fastfocal: Fast Multi-scale Raster Extraction and Moving Window Analysis with FFT"
subtitle: "Version 0.1.3 (2025)"
author: "Ho Yi Wan"
output: rmarkdown::html_vignette
bibliography: references.bib
vignette: >
  %\VignetteIndexEntry{fastfocal: Fast Multi-scale Raster Extraction and Moving Window Analysis with FFT}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = TRUE,
  fig.width = 6,
  fig.height = 4,
  out.width = "100%"
)
library(terra)
library(fastfocal)
```

## What is fastfocal?

I created the `fastfocal` R package (Wan 2025) for multi-scale ecological modeling that requires focal statistics at various spatial scales. As rasters and kernels grow larger, traditional methods can become slow. `fastfocal` is designed for high-performance focal and buffer-based raster processing. It provides fast moving-window operations and point-based extraction with a variety of window (kernel) shapes. The package can automatically switch to a Fast Fourier Transform (FFT) backend for large-kernel operations, and otherwise use the C++ backend via `terra` for smaller cases. It is built on `terra` and works with `SpatRaster` and `SpatVector` objects.

This vignette walks through the core functionality using simulated raster layers and step-by-step examples. We highlight how `fastfocal()` and `fastextract()` support multi-scale analysis with flexible kernel definitions.

---

## Key features

- Multi-layer raster support: apply focal operations across stacked rasters.
- Flexible kernels: choose from many window types (e.g., circle, gaussian, logistic, quartic).
- Map-unit radius: specify smoothing or extraction radius directly in meters (map units).
- Point-based extraction: extract raster summaries at points with optional buffer windows.
- FFT backend: automatically selected for large kernels; C++ backend otherwise.
- Identical geometry output: results maintain resolution, extent, and CRS.

---

## Installation

```{r, eval=FALSE}
# install.packages("devtools")
# devtools::install_github("hoyiwan/fastfocal")
```

```{r}
library(fastfocal)
library(terra)
```

---

## Simulating rasters for this vignette

We generate three example rasters with values resembling ecological variables. Users can replace these with their own data.

```{r simulate-coarse-fine}
set.seed(888)
ndvi_coarse  <- rast(matrix(runif(100, 0.1, 0.9), 10, 10))            # NDVI (0.1 - 0.9)
slope_coarse <- rast(matrix(rbeta(100, 2, 5) * 40, 10, 10))            # slope (degrees)
temp_coarse  <- rast(matrix(seq(12, 24, length.out = 100) + rnorm(100, 0, 1), 10, 10))  # temperature with gradient

# Disaggregate to finer resolution
ndvi  <- disagg(ndvi_coarse, fact = 10, method = "bilinear")
slope <- disagg(slope_coarse, fact = 10, method = "bilinear")
temp  <- disagg(temp_coarse, fact = 10, method = "bilinear")

# Set a 3 km by 3 km extent for all layers (projected units, e.g., meters)
ext(ndvi) <- ext(0, 3000, 0, 3000)
ext(slope) <- ext(ndvi)
ext(temp) <- ext(ndvi)
crs(ndvi) <- "EPSG:32611"  # UTM Zone 11N
crs(slope) <- crs(ndvi)
crs(temp) <- crs(ndvi)

# Combine into a multi-layer SpatRaster
r <- c(ndvi, slope, temp)
names(r) <- c("ndvi", "slope", "temp")

plot(r)
```

---

## Point and buffer-based extraction

Use `fastextract()` to extract values at points and, optionally, within buffers of different sizes.

```{r simulate-points}
# Simulate some points
pts <- vect(matrix(c(500,500, 1500,1500, 2500,2500), ncol=2, byrow=TRUE), type="points", crs=crs(r))

oldpar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
plot(r[[1]], main = "NDVI");  plot(pts, add = TRUE, pch = 20, cex = 2)
plot(r[[2]], main = "Slope"); plot(pts, add = TRUE, pch = 20, cex = 2)
plot(r[[3]], main = "Temp");  plot(pts, add = TRUE, pch = 20, cex = 2)
par(oldpar)
```

```{r extract}
fastextract(r, pts)                               # default d = 0 (extract at the point)
fastextract(r, pts, d = c(0, 100, 500))           # multiple distances
fastextract(r, pts, d = c(0, 100, 500, 1000), w = "gaussian", fun = "sd")
```

---

## Moving-window analysis

Apply `fastfocal()` to compute moving-window statistics over the stack using different kernels and summary functions.

```{r focal}
# Focal means using a circular window (300 m radius)
f_circ <- fastfocal(r, d = 300, w = "circle",   fun = "mean")

# Focal means using a gaussian window (300 m radius)
f_gaus <- fastfocal(r, d = 300, w = "gaussian", fun = "mean")

oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 3), mar = c(2, 2, 2, 2))
plot(r[[1]], main = "NDVI");  plot(f_circ[[1]], main = "NDVI - Circle");  plot(f_gaus[[1]], main = "NDVI - Gaussian")
plot(r[[2]], main = "Slope"); plot(f_circ[[2]], main = "Slope - Circle"); plot(f_gaus[[2]], main = "Slope - Gaussian")
plot(r[[3]], main = "Temp");  plot(f_circ[[3]], main = "Temp - Circle");  plot(f_gaus[[3]], main = "Temp - Gaussian")
par(oldpar)
```

---

## Comparing statistics on one layer

```{r multi-stat}
f_mean   <- fastfocal(r, d = 300, fun = "mean",   na.rm = TRUE)
f_sd     <- fastfocal(r, d = 300, fun = "sd",     na.rm = TRUE)
f_median <- fastfocal(r, d = 300, fun = "median", na.rm = TRUE)

oldpar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
plot(r[[1]],       main = "NDVI")
plot(f_mean[[1]],  main = "NDVI - Mean (300 m)")
plot(f_sd[[1]],    main = "NDVI - SD (300 m)")
plot(f_median[[1]],main = "NDVI - Median (300 m)")
par(oldpar)
```

---

## Visualizing supported kernels

This section visualizes the built-in window types supported by `fastfocal_weights()`. We keep the radius modest so the code runs quickly in the vignette.

```{r visualize-kernels}
center_r <- rast(ext(0, 90, 0, 90), resolution = 30, crs = "EPSG:32611")
kernel_types <- c("circle", "rectangle", "gaussian", "pareto", "idw", "exponential",
                  "triangular", "cosine", "logistic", "cauchy", "quartic", "epanechnikov")

pad_kernel <- function(k, pad = 2) {
  nr <- nrow(k); nc <- ncol(k)
  out <- matrix(0, nrow = nr + 2 * pad, ncol = nc + 2 * pad)
  out[(pad + 1):(pad + nr), (pad + 1):(pad + nc)] <- k
  out
}

oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 4), mar = c(2, 2, 2, 2))
for (w in kernel_types) {
  k_raw <- fastfocal_weights(center_r, d = 150, w = w, normalize = FALSE)  # ~11x11 at 30 m
  k_pad <- pad_kernel(k_raw, pad = 2)
  k <- k_pad / max(k_pad, na.rm = TRUE)
  image(k, main = w, col = hcl.colors(20, "YlOrRd", rev = TRUE), zlim = c(0, 1), asp = 1, cex.main = 1.2)
}
par(oldpar)
```

---

## Notes

- Multi-layer support is built into `fastfocal()`.
- Input can be a single-layer raster or a stack.
- Output preserves layer names and returns a `SpatRaster` with the same geometry.
- Kernel size `d` is in map units (e.g., meters).

---

## Explore more

See the benchmark comparison: `vignettes/benchmark.Rmd` (built as `inst/doc/benchmark.html`).

---

## Session info

```{r}
sessionInfo()
```

## Citation

To cite the package:

Wan, H. Y. (2025). *fastfocal: Fast Multi-scale Raster Extraction and Moving Window Analysis with FFT*. R package version 0.1.3. Zenodo. https://doi.org/10.5281/zenodo.17074691


```{r}
citation("fastfocal")
```