p <- expand.grid(143:149, -44:-40)
plot(p, asp = 1/cos(-40 * pi/180))
maps::map(add = TRUE)
text(p, lab = geographiclib::mgrs_fwd(as.matrix(p) * 1.0, precision = 0), pos = 1, cex = .7)
Contrast stretching for a Sentinel-2 scene
Finding a scene
First we write a quick plot to find the MGRS code at (bottom left corner of) the tile we want.
Then look up a scene for that tile on a date we already know has one.
## remotes::install_github("hypertidy/sds")
qu <- sds::stacit(c("55GDN"), "2026-02-27")
## get the STAC monstrosity and find the href:
js <- jsonlite::fromJSON(qu)
sprintf("/vsicurl/%s", js$features$assets$red$href)[1] "/vsicurl/https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/55/G/DN/2026/2/S2C_T55GDN_20260227T000650_L2A/B04.tif"Stretch functions
All functions return values in [0, 1]. Named with the s2_ prefix to avoid clashing with terra::stretch().
library(terra)terra 1.8.97# 1. Linear percent clip
# low/high are the natural UX knobs; 2%/98% is a solid default,
# tighter for low-contrast scenes.
s2_stretch_linear <- function(x, low = 0.02, high = 0.98) {
v <- terra::values(x, na.rm = TRUE)
lo <- quantile(v, low, na.rm = TRUE)
hi <- quantile(v, high, na.rm = TRUE)
terra::clamp(x, lo, hi, values = TRUE) |>
(\(r) (r - lo) / (hi - lo))()
}
# 2. Log stretch (per-band minimum offset, guards log(0))
# Use only on integer overview values, not native-res reflectance — see gotcha below.
s2_stretch_log <- function(x, offset = NULL) {
if (is.null(offset)) {
mn <- min(terra::global(x, "min", na.rm = TRUE)[[1]])
offset <- if (mn <= 0) abs(mn) + 1 else 0
}
log1p(x + offset) |> s2_stretch_linear()
}
# 3. Log-relative: subtract global minimum before log.
# + 1 guards against log(0) when a pixel equals the minimum exactly.
s2_stretch_log_relative <- function(x, low = 0.02, high = 0.98) {
mn <- min(terra::global(x, "min", na.rm = TRUE)[[1]])
log(x - mn + 1) |> s2_stretch_linear(low, high)
}
# Joint variant: single global minimum AND a single quantile clip across all
# bands combined — preserves inter-band luminance relationships for true-colour RGB.
# In practice the visual difference from per-band is subtle for well-exposed scenes;
# the real fix for colour balance is the shared minimum subtraction.
s2_stretch_log_relative_joint <- function(x, low = 0.02, high = 0.98) {
mn <- min(terra::global(x, "min", na.rm = TRUE)[[1]])
logged <- log(x - mn + 1)
v <- as.vector(terra::values(logged, na.rm = TRUE)) # flatten all bands
lo <- quantile(v, low, na.rm = TRUE)
hi <- quantile(v, high, na.rm = TRUE)
terra::clamp(logged, lo, hi, values = TRUE) |>
(\(r) (r - lo) / (hi - lo))()
}
# 4. Square root — gentler than log, good middle ground.
s2_stretch_sqrt <- function(x) {
mn <- min(terra::global(x, "min", na.rm = TRUE)[[1]])
offset <- if (mn < 0) abs(mn) else 0
sqrt(x + offset) |> s2_stretch_linear()
}
# 5. Histogram equalisation.
# Good when interesting detail spans simultaneously deep shadow and bright areas.
s2_stretch_histeq <- function(x, n = 256) {
v <- terra::values(x, na.rm = TRUE)
breaks <- quantile(v, probs = seq(0, 1, length.out = n + 1), na.rm = TRUE)
out <- x
terra::values(out) <- findInterval(terra::values(x), unique(breaks)) /
length(unique(breaks))
out
}Dispatcher
joint defaults to TRUE for multi-band input so that RGB composites share a common stretch across all three channels.
s2_stretch <- function(x,
method = c("linear", "log", "sqrt", "histeq", "log_relative"),
low = 0.02, high = 0.98,
joint = terra::nlyr(x) > 1, ...) {
method <- match.arg(method)
switch(method,
linear = if (joint) {
v <- terra::values(x, na.rm = TRUE)
lo <- quantile(v, low, na.rm = TRUE)
hi <- quantile(v, high, na.rm = TRUE)
terra::clamp(x, lo, hi, values = TRUE) |>
(\(r) (r - lo) / (hi - lo))()
} else s2_stretch_linear(x, low, high),
log_relative = if (joint) s2_stretch_log_relative_joint(x, low, high)
else s2_stretch_log_relative(x, low, high),
log = s2_stretch_log(x, ...),
sqrt = s2_stretch_sqrt(x),
histeq = s2_stretch_histeq(x, ...)
)
}
# terra::plotRGB() requires [0, 255]; terra::stretch() is used solely for that
# byte-scaling side effect, applied to already-stretched [0, 1] data.
pltrgb <- function(x, ...) {
terra::plotRGB(terra::stretch(x), ...)
}Comparing stretch methods
library(terra)
band <- c("B04", "B03", "B02")
dsn <- glue::glue("vrt:///vsicurl/https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/55/G/DN/2026/2/S2C_T55GDN_20260227T000650_L2A/{band}.tif?ovr=2")
r <- terra::rast(dsn)
mar <- par("mar")
op <- par(mfcol = c(2, 3), mar = mar + c(1, 0, 1, 0) * 2)
pltrgb(s2_stretch(r, "linear"), main = "linear")
pltrgb(s2_stretch(r, "log"), main = "log")
pltrgb(s2_stretch(r, "sqrt"), main = "sqrt")
pltrgb(s2_stretch(r, "histeq"), main = "histeq")
pltrgb(s2_stretch(r, "log_relative", joint = FALSE), main = "log_relative")
pltrgb(s2_stretch(r, "log_relative"), main = "log_relative_joint")