I am downscaling (i.e., increasing the spatial resolution) VIIRS nighttime lights imagery. The transfer function between the coarse and fine resolution is approximated by a Gaussian filter whose width σ varies with the per‑column viewing angle θ. This is an analytical, geometric relationship, not an estimate.
In the cross‑track direction (left→right), the sigma scales as
where θ is the per‑pixel viewing angle. The filter should act only horizontally (along x), pooling all y‑values within the column uniformly. My earlier attempt simply multiplied the raster values by
but I now realise the operation is a convolution with a 1‑D Gaussian kernel whose width depends on the column’s θ.
Given a value raster x and an angle raster theta with identical dimensions, what is the correct and efficient way in R (using terra) to apply a 1‑D Gaussian blur on each row, where the kernel’s standard deviation is determined by theta at the central pixel’s column?
Current (incomplete) attempt:
library(terra)
# ---- 1. Create example rasters ----
set.seed(42)
nrows <- 5
ncols <- 200
x <- rast(nrows = nrows, ncols = ncols, vals = runif(nrows * ncols))
theta_vals <- rep(seq(20, 26, length.out = ncols), each = nrows)
theta <- rast(nrows = nrows, ncols = ncols, vals = theta_vals)
# ---- 2. Spatially varying Gaussian blur function ----
# sigma0 = standard deviation of Gaussian at nadir (in pixel units)
# theta: raster of angles in degrees (same dimensions as x)
anisotropic_blur <- function(x, theta, sigma0 = 10) {
# Convert to matrix (rows = y, cols = x)
mat_x <- as.matrix(x, wide = TRUE)
mat_angle <- as.matrix(theta, wide = TRUE) # same dimensions
nr <- nrow(mat_x)
nc <- ncol(mat_x)
mat_out <- matrix(NA_real_, nr, nc)
# For each row, apply convolution with column-dependent sigma
for (r in 1:nr) {
row_vals <- mat_x[r, ]
for (c in 1:nc) {
theta_c <- mat_angle[r, c] # angle at this pixel (degrees)
sigma <- sigma0 / (cos(theta_c * pi/180)^2)
halfwin <- ceiling(3 * sigma) # kernel half‑width
col_idx <- max(1, c - halfwin) : min(nc, c + halfwin)
dist <- abs(col_idx - c)
w <- exp(-0.5 * (dist / sigma)^2)
w <- w / sum(w)
mat_out[r, c] <- sum(w * row_vals[col_idx])
}
}
# Return as raster with same properties
rast(mat_out, crs = crs(x), ext = ext(x))
}
# ---- 3. Apply the filter ----
x_blurred <- anisotropic_blur(x, theta, sigma0 = 10)
# ---- 4. Plot side‑by‑side ----
par(mfrow = c(1, 2))
plot(x, main = "Original predictor")
plot(x_blurred, main = expression("Filtered: Gaussian "*sigma*" = "*sigma[0]/cos^2*theta))
I need to actually perform a 1‑D Gaussian convolution along rows with a column‑varying (\sigma). What is an idiomatic way to do this in terra?