Skip to content

How to "translate" math to R code to apply blurring to spatraster

1 message · Nikolaos Tziokas

Original
Nikolaos Tziokas
# 
I manage to create a code the looks about right, but I am not 100% sure
about it.

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)

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)
values(va_viirs) <- rep(seq(22.5, 24.5, length.out=ncol(va_viirs)),
times=nrow(va_viirs))
va_high_res <- resample(va_viirs, high_res_covariate, method="near")

R_E <- 6371         # km
h_VIIRS <- 824      # km (Note: image shows 833)
sigma_nadir <- 0.742 # km

# 1. Calculate sigma_cross per pixel (Your exact function)
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 km -> m -> pixels
}
sigma_cross_pixels <- app(va_high_res, calc_sigma_cross)



# 2. SEPARABLE FILTERING: Y-Direction (Constant sigma = 2 pixels)
sigma_along_pixels <- 2
radius_y <- ceiling(3 * sigma_along_pixels) # 3 standard deviations
y_seq <- seq(-radius_y, radius_y)
kernel_y <- dnorm(y_seq, mean = 0, sd = sigma_along_pixels)
kernel_y <- matrix(kernel_y / sum(kernel_y), ncol = 1) # Vertical kernel

# Apply vertical blur
y_blurred <- focal(high_res_covariate, w = kernel_y, na.rm = TRUE)



# 3. SEPARABLE FILTERING: X-Direction (Varying sigma)
# Extract to matrices for fast row-by-row processing
val_mat <- as.matrix(y_blurred, wide = TRUE)
sig_mat <- as.matrix(sigma_cross_pixels, wide = TRUE)

out_mat <- matrix(NA, nrow = nrow(val_mat), ncol = ncol(val_mat))
cols <- 1:ncol(val_mat)

# Pre-calculate a meshgrid of column indices
j_mat <- matrix(cols, nrow = length(cols), ncol = length(cols), byrow =
TRUE)  # Source columns
i_mat <- matrix(cols, nrow = length(cols), ncol = length(cols), byrow =
FALSE) # Target columns

for (r in 1:nrow(val_mat)) {
  row_vals <- val_mat[r, ]
  row_sigs <- sig_mat[r, ]

  # Broadcast the target pixel's standard deviation across the weight matrix
  sig_i <- matrix(row_sigs, nrow = length(cols), ncol = length(cols), byrow
= FALSE)

  # Calculate Gaussian weights (W[i,j] is the weight of source j on target
i)
  W <- exp(- ((j_mat - i_mat)^2) / (2 * sig_i^2))

  # Normalize rows so weights sum to 1 (this naturally handles image edges)
  W <- W / rowSums(W)

  # Apply weights to the row values using matrix multiplication
  out_mat[r, ] <- as.numeric(W %*% row_vals)
}

# Assign back to a SpatRaster
final_filtered <- y_blurred
values(final_filtered) <- as.vector(t(out_mat)) # t() ensures correct
row-major assignment

plot(final_filtered, main="Anisotropic Filtered Output")

# Ensure consistent color scale
range_vals <- range(
  values(high_res_covariate),
  values(final_filtered),
  na.rm = TRUE
)

# Side-by-side layout
par(mfrow = c(1, 2), mar = c(4, 4, 3, 5))

plot(
  high_res_covariate,
  main = "Original Predictor (100m)",
  col = hcl.colors(100, "YlOrRd"),
  zlim = range_vals
)

plot(
  final_filtered,
  main = "Anisotropic Gaussian Blurred",
  col = hcl.colors(100, "YlOrRd"),
  zlim = range_vals
)

par(mfrow = c(1, 1))  # reset layout

???? ??? 27 ??? 2026 ???? 5:56??.?., ?/? Nikolaos Tziokas <
nikos.tziokas at gmail.com> ??????: