I need to apply an anisotropic filter to a raster where:
?_along (along-track, y-direction) is constant = 2 pixels
?_cross (cross-track, x-direction) varies per pixel based on viewing angle
I have the mathematical formulas to calculate ?_cross for each pixel:
Given:
R_E = 6371 km (Earth's radius)
h_VIIRS = 824 km (satellite altitude)
? = viewing angle from nadir (varies per pixel)
?_nadir = 742 m
D_view = (R_E + h_VIIRS) ? cos(?) - ?[R_E? - (R_E + h_VIIRS)? ? sin?(?)] ?
= arcsin[(R_E + h_VIIRS)/R_E ? sin(?)] ?_cross = ?_nadir ? (D_view /
h_VIIRS) ? [1 / cos(?)]
Please see this link (https://www.desmos.com/calculator/zwvoyxzpoz) for the
actual maths I need to translate into R code to blur the spatraster.
My current pixel-by-pixel nested loop implementation produces diagonal
striping artifacts, suggesting I'm not correctly translating the
mathematics into working code.
Reproducible example:
library(terra)
# High-resolution covariate (100m)
set.seed(123)
high_res_covariate <- rast(nrows=230, ncols=255,
xmin=17013000, xmax=17038500,
ymin=-3180000, ymax=-3157000,
crs="EPSG:3857")
res(high_res_covariate) <- c(100, 100)
values(high_res_covariate) <- runif(ncell(high_res_covariate), 0, 100)
# Viewing angle raster (500m, varies left to right)
viirs_ntl <- rast(nrows=46, ncols=51,
xmin=17013000, xmax=17038500,
ymin=-3180000, ymax=-3157000,
crs="EPSG:3857")
res(viirs_ntl) <- c(500, 500)
va_viirs <- rast(viirs_ntl)
va_values <- rep(seq(22.5, 24.5, length.out=ncol(va_viirs)),
times=nrow(va_viirs))
values(va_viirs) <- va_values
va_high_res <- resample(va_viirs, high_res_covariate, method="near")
# Parameters
R_E <- 6371 # km
h_VIIRS <- 824 # km
sigma_nadir <- 0.742 # km
# Calculate sigma_cross per pixel
calc_sigma_cross <- function(theta_deg) {
theta_rad <- theta_deg * pi / 180
D_view <- (R_E + h_VIIRS) * cos(theta_rad) -
sqrt(R_E^2 - (R_E + h_VIIRS)^2 * sin(theta_rad)^2)
beta <- asin((R_E + h_VIIRS) / R_E * sin(theta_rad))
sigma_cross <- sigma_nadir * (D_view / h_VIIRS) * (1 / cos(beta))
return(sigma_cross * 1000 / 100) # Convert to pixels
}
sigma_cross_pixels <- app(va_high_res, calc_sigma_cross)
sigma_along_pixels <- 2 # constant
How do I use terra::focal() with these spatially-varying sigma values?
SessionInfo
R version 4.5.2 (2025-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)
Matrix products: default
LAPACK version 3.12.1
locale:
[3] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United
States.utf8 LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C LC_TIME=English_United
States.utf8
time zone: Europe/Budapest
tzcode source: internal
attached base packages:
[3] stats graphics grDevices utils datasets methods base
other attached packages:
[3] terra_1.8-93
loaded via a namespace (and not attached):
[3] compiler_4.5.2 cli_3.6.5 ragg_1.5.0 tools_4.5.2
rstudioapi_0.18.0 Rcpp_1.1.1 codetools_0.2-20
[8] textshaping_1.0.4 lifecycle_1.0.5 rlang_1.1.7 systemfonts_1.3.1