A SAR image looks stretched or smeared in mountains

Problem statement

In hill country a radar image stops looking like a map. Slopes facing the satellite are compressed into bright bands, ridges appear to lean towards the sensor, valleys behind them are empty, and a terrain-corrected product still shows one side of every ridge brighter than the other. Overlay it on a topographic map and features are displaced; difference two passes and every slope has changed.

None of this is a processing fault to be fixed after the fact. It is the geometry of a side-looking sensor meeting relief, and the remedies are choosing the right product, masking the pixels that carry no information, and using the pass that sees each slope best.

In the Sierra Nevada, about 1% of pixels fall in layover or shadow for a given pass at the middle of the swath — but they are concentrated exactly on the steep slopes that a mountain study is usually about, and which slopes they are changes with the orbit.

Quick answer

Compute the layover and shadow mask from a DEM and the look direction, and mask before analysis:

import numpy as np, rasterio, pystac_client, planetary_computer
from rasterio.enums import Resampling
from rasterio.vrt import WarpedVRT
from sar_io import STAC, search, passes, grid, read_pass, to_db

bbox = [-3.45, 37.00, -3.25, 37.12]
g = grid(bbox, "EPSG:25830", res=20)
cat = pystac_client.Client.open(STAC, modifier=planetary_computer.sign_inplace)
dem = np.full((g["height"], g["width"]), np.nan, "float32")
for it in cat.search(collections=["cop-dem-glo-30"], bbox=bbox).items():
    with rasterio.open(it.assets["data"].href) as src, WarpedVRT(
            src, crs=g["crs"], transform=g["transform"], width=g["width"],
            height=g["height"], resampling=Resampling.bilinear, nodata=np.nan) as v:
        a = v.read(1)
    dem[np.isnan(dem) & np.isfinite(a)] = a[np.isnan(dem) & np.isfinite(a)]
gy, gx = np.gradient(dem, 20)

def masks(look_azimuth_deg, incidence_deg=39):
    az = np.radians(look_azimuth_deg)
    rise = gx * np.sin(az) + (-gy) * np.cos(az)          # slope along the look direction
    layover = rise > np.tan(np.radians(incidence_deg))
    shadow = -rise > np.tan(np.radians(90 - incidence_deg))
    return layover, shadow

lay_a, sha_a = masks(80)                                  # ascending, looking east
lay_d, sha_d = masks(280)                                 # descending, looking west
print(f"ascending  unusable {np.mean(lay_a | sha_a):.2%}")
print(f"descending unusable {np.mean(lay_d | sha_d):.2%}")
print(f"unusable in both    {np.mean((lay_a | sha_a) & (lay_d | sha_d)):.2%}")

ascending  unusable 0.92%
descending unusable 1.31%
unusable in both    0.20%

The pixels that fail both passes are the ones no Sentinel-1 geometry can see. Everywhere else, one of the two passes is usable.

Terrain-corrected data, DEM masks, both pass directions

1. Use a terrain-corrected product

RTC places pixels correctly and removes most of the slope brightness. Geocoded-only or GRD data in mountains will not line up with anything.

2. Build the masks from a DEM

Compute the slope along the look direction and compare it with the incidence angle. Steeper towards the sensor is layover; steeper away than 90° minus incidence is shadow.

3. Mask before every statistic

Layover mixes several surfaces into one pixel and shadow contains only noise. Neither should enter a mean, a threshold or a training set.

4. Use both pass directions

A slope hidden from one direction usually faces the other. Compute per pixel which pass has the local incidence angle closest to the nominal, and take that one.

5. Expect residual brightness differences

Even in valid pixels, RTC leaves slopes facing the sensor 1–3 dB brighter. Do not compare slopes with flats, or ascending with descending, without allowing for it.

6. Check the DEM

Every correction and every mask inherits the DEM's errors. Copernicus DEM at 30 m is the usual choice; in very steep terrain its own errors matter.

Which pass sees each slope better

The data used on this page

Sentinel-1 radiometrically terrain-corrected backscatter (sentinel-1-rtc) from the Microsoft Planetary Computer, read through sar_io.py — the small module from the reading guide that searches, groups frames into passes and mosaics each pass onto a fixed grid.

%%writefile sar_io.py
"""Find, read and mosaic Sentinel-1 RTC scenes from Microsoft Planetary Computer."""
from collections import defaultdict

import numpy as np
import planetary_computer
import pystac_client
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.vrt import WarpedVRT
from rasterio.warp import transform_bounds

STAC = "https://planetarycomputer.microsoft.com/api/stac/v1"


def search(bbox, dates, orbit=None, direction=None, collection="sentinel-1-rtc"):
    """Items over bbox (lon/lat) in a date range, filtered by relative orbit and pass direction."""
    cat = pystac_client.Client.open(STAC, modifier=planetary_computer.sign_inplace)
    items = cat.search(collections=[collection], bbox=bbox, datetime=dates).item_collection()
    keep = [i for i in items
            if (orbit is None or i.properties["sat:relative_orbit"] == orbit)
            and (direction is None or i.properties["sat:orbit_state"] == direction)]
    return sorted(keep, key=lambda i: i.datetime)


def passes(items):
    """Group frames into passes: same day, same relative orbit."""
    groups = defaultdict(list)
    for i in items:
        groups[(i.datetime.date(), i.properties["sat:relative_orbit"])].append(i)
    return dict(sorted(groups.items()))


def grid(bbox, crs, res=10):
    """A fixed target grid: every pass is read onto this, so arrays line up."""
    left, bottom, right, top = transform_bounds(4326, crs, *bbox)
    width, height = int((right - left) // res), int((top - bottom) // res)
    return {"crs": crs, "transform": from_origin(left, top, res, res),
            "width": width, "height": height}


def read_pass(frames, band, g):
    """Mosaic one pass's frames for one band onto grid g. Linear power; NaN where no data."""
    out = np.full((g["height"], g["width"]), np.nan, dtype="float32")
    for item in frames:
        with rasterio.open(item.assets[band].href) as src, WarpedVRT(
                src, crs=g["crs"], transform=g["transform"], width=g["width"],
                height=g["height"], resampling=Resampling.bilinear, nodata=np.nan) as vrt:
            a = vrt.read(1)
        a[a <= 0] = np.nan                       # -32768 nodata and zero-padding
        fill = np.isnan(out) & np.isfinite(a)
        out[fill] = a[fill]
    return out


def to_db(linear):
    return 10 * np.log10(linear)

Overwriting sar_io.py

Example 1 — which pass sees each slope better

slope = np.degrees(np.arctan(np.hypot(gx, gy)))
aspect = (np.degrees(np.arctan2(-gx, gy)) + 360) % 360

def local_incidence(look_azimuth_deg, incidence_deg=39):
    """Angle between the look direction and the surface normal, in degrees."""
    az, inc = np.radians(look_azimuth_deg), np.radians(incidence_deg)
    look = np.array([-np.sin(inc) * np.sin(az), -np.sin(inc) * np.cos(az), np.cos(inc)])
    n = np.stack([gx, -gy, np.ones_like(gx)])
    n = n / np.linalg.norm(n, axis=0)
    return np.degrees(np.arccos(np.clip(np.einsum("i,ijk->jk", look, n), -1, 1)))

li_a, li_d = local_incidence(80), local_incidence(280)
better_asc = np.abs(li_a - 39) < np.abs(li_d - 39)
print(f"ascending is the better geometry for {np.nanmean(better_asc):.0%} of pixels")
steep = slope > 25
print(f"on slopes over 25°: ascending better for {np.nanmean(better_asc[steep]):.0%}, "
      f"descending for {np.nanmean(~better_asc[steep]):.0%}")

ascending is the better geometry for 54% of pixels
on slopes over 25°: ascending better for 51%, descending for 49%

Example 2 — what masking changes

P = passes(search(bbox, "2024-08-01/2024-08-07"))
db = {orbit: to_db(read_pass(frames, "vv", g)) for (day, orbit), frames in P.items() if orbit in (1, 81)}
for orbit, img in db.items():
    bad = (lay_a | sha_a) if orbit == 1 else (lay_d | sha_d)
    print(f"orbit {orbit:2d}: all pixels median {np.nanmedian(img):6.2f} dB   "
          f"masked median {np.nanmedian(img[~bad]):6.2f} dB   "
          f"layover/shadow pixels median {np.nanmedian(img[bad]):6.2f} dB")

orbit  1: all pixels median  -8.10 dB   masked median  -8.11 dB   layover/shadow pixels median  -7.50 dB
orbit 81: all pixels median  -8.11 dB   masked median  -8.14 dB   layover/shadow pixels median  -4.48 dB

Example 3 — the best-geometry mosaic

asc, dsc = db[1], db[81]
best = np.where(better_asc, asc, dsc)
bad_both = (lay_a | sha_a) & (lay_d | sha_d)
best = np.where(bad_both, np.nan, best)
print(f"combined: {np.isfinite(best).mean():.1%} of pixels usable, "
      f"median {np.nanmedian(best):.2f} dB")
print(f"difference between the two passes on steep slopes: "
      f"median {np.nanmedian(np.abs(asc - dsc)[steep]):.2f} dB")

combined: 99.8% of pixels usable, median -8.34 dB
difference between the two passes on steep slopes: median 2.54 dB

A mosaic like this is for interpretation, not for a time series: it mixes geometries by construction, and the seam follows the terrain.

Why layover cannot be undone

Why the distortion cannot be undone

Layover is a many-to-one mapping: several parts of the ground return at the same range and are summed into one pixel. Terrain correction can place that pixel correctly, but the sum cannot be separated into its parts. Shadow is the opposite — nothing was measured — and no interpolation can invent it. Both are information that the acquisition never recorded.

Why the mask changes with the pass

A slope facing east is compressed under an ascending pass looking east and stretched under a descending pass looking west. Because Sentinel-1 is right-looking on both, the two geometries are roughly opposite, which is exactly why combining them works: a slope is rarely too steep for both.

Why residual brightness differences remain

Terrain flattening normalises by the illuminated area, which handles the geometric part of the effect. It does not model how a rough surface scatters differently when the wave arrives at a shallow local angle. That residual is what leaves sensor-facing slopes a decibel or two brighter in even the best products.

Swath incidence, snow and DEM voids

• Incidence varies across the swath, so the same slope can be usable at far range and in layover at near range.
• Cities have their own layover: tall buildings lean towards the sensor over the streets in front.
• Snow and ice add their own backscatter changes on top of the geometry.
• DEM voids in steep terrain propagate into both the correction and the mask.
• Sentinel-1 provides no layover mask on Planetary Computer's RTC items; compute one.
• For InSAR, the same geometry limits apply and coherence is lost first on the steepest slopes.

Internal links

• SAR geometry explained — layover, foreshortening and shadow
• Sentinel-1 products explained — what RTC corrects
• SAR scenes from different dates do not line up — mixing geometries
• Digital elevation models explained — the DEM behind every correction
• How to calculate slope and aspect in Python — the terrain derivatives used here
• How to create a hillshade in Python — the optical analogue of illumination

FAQ

Why does my SAR image look stretched in mountains?

Foreshortening: slopes facing the sensor span less range than ground distance, so they compress into fewer, brighter pixels.

Can terrain correction remove layover?

No. It places pixels correctly but cannot separate signals that were summed into one.

How do I find layover and shadow pixels?

Compute the slope along the look direction from a DEM and compare it with the incidence angle.

Should I combine ascending and descending passes?

For interpretation and coverage, yes. For a time series, no.

How much of a mountain scene is unusable?

Around one percent at mid-swath in the Sierra Nevada, concentrated on the steepest slopes.

Does the DEM matter?

Yes. Both the correction and the masks inherit its errors.