SAR values come out negative, huge or averaged wrongly

Problem statement

A Sentinel-1 array should hold small positive numbers, and yet the minimum is −32768, the maximum is in the thousands, a histogram is a single spike at zero, or a zonal mean comes out several decibels away from what the same zone looks like on screen. These are four faces of one problem: backscatter arrives in linear power with a nodata convention, and every one of those symptoms comes from treating it as something else.

On a Sentinel-1 RTC frame over Valencia, a window at the frame edge is 32.7% nodata; left in, it drags the window's mean to −10,715. Fifty-one of its valid pixels sit above 0 dB, which is a corner reflector and not an error. And averaging the open-water class of the same scene in decibels rather than in power moves its mean by 8.6 dB.

Quick answer

Mask nodata, keep linear power for arithmetic, and convert to decibels only for display and thresholds:

import numpy as np
import rasterio
from sar_io import search

item = search([-0.46, 39.25, -0.33, 39.45], "2024-10-20/2024-10-20", orbit=110)[0]
with rasterio.open(item.assets["vv"].href) as src:
    raw = src.read(1, window=((192, 592), (4800, 5200)))    # straddles the frame edge
    nodata = src.nodata

print(f"raw: min {raw.min():.0f}, max {raw.max():.1f}, nodata {nodata}, "
      f"pixels at nodata {np.mean(raw == nodata):.1%}, at zero {np.mean(raw == 0):.1%}")

vv = np.where(raw > 0, raw, np.nan)                 # nodata and zeros out
db = 10 * np.log10(vv)
print(f"valid: {np.nanmin(vv):.5f} to {np.nanmax(vv):.1f} linear "
      f"= {np.nanmin(db):.1f} to {np.nanmax(db):.1f} dB")

raw: min -32768, max 2.9, nodata -32768.0, pixels at nodata 32.7%, at zero 0.0%
valid: 0.00734 to 2.9 linear = -21.3 to 4.7 dB

Positive decibels are not an error. A corner formed by two flat surfaces can return more power than the diffuse reference the coefficient is defined against.

Read nodata, drop zeros, compute in linear power

1. Read the nodata value from the file

Do not hard-code it. src.nodata is −32768.0 for these products; other providers differ.

2. Drop zeros as well

Frame padding and some processing steps leave exact zeros. They are not measurements, and their logarithm is minus infinity.

3. Expect a wide range

Valid linear values span roughly 0.0001 to several thousand. Water sits near 0.005 (−23 dB), fields near 0.1 (−10 dB), and corner reflectors above 1 (positive dB).

4. Do arithmetic in linear power

Means, composites, resampling and filters operate on power. Convert the result once.

5. Check the scale convention

Power values use 10·log₁₀; amplitude values use 20·log₁₀. If a provider's decibels are twice what you expect, the product holds amplitude.

6. Check for pre-scaled integers

Some archives store decibels as scaled integers with a scale_factor. Reading them as power gives values of thousands that are not corner reflectors.

What surviving nodata does

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 — what nodata does if it survives

bad_mean = raw.mean()
good_mean = np.nanmean(vv)
print(f"mean with nodata included: {bad_mean:,.1f} (in linear units!)   "
      f"mean of valid pixels: {good_mean:.4f} = {10 * np.log10(good_mean):.2f} dB")
with np.errstate(divide="ignore", invalid="ignore"):
    print("log of nodata:", np.log10(np.array([-32768.0, 0.0])))

mean with nodata included: -10,714.9 (in linear units!)   mean of valid pixels: 0.1164 = -9.34 dB
log of nodata: [ nan -inf]

Example 2 — the bright pixels are real

bright = vv[vv > 1]
print(f"{bright.size} pixels above 1 (0 dB) out of {np.isfinite(vv).sum():,}: "
      f"up to {np.nanmax(vv):.0f} linear = {10 * np.log10(np.nanmax(vv)):.1f} dB")

51 pixels above 1 (0 dB) out of 107,681: up to 3 linear = 4.7 dB

A handful of pixels far above the rest is what a metal roof, a ship or a pylon looks like. Clipping them changes every mean that includes them.

Example 3 — the decibel mean, once more

zone = vv[np.isfinite(vv)]
print(f"mean of power -> {10 * np.log10(zone.mean()):.2f} dB")
print(f"mean of dB    -> {np.mean(10 * np.log10(zone)):.2f} dB")
print(f"median (same either way) -> {10 * np.log10(np.median(zone)):.2f} dB")

mean of power -> -9.34 dB
mean of dB    -> -9.94 dB
median (same either way) -> -9.92 dB

Why nodata is a number and bright values are real

Why nodata is −32768 and not NaN

The value comes from the integer conventions of earlier products and has been kept for compatibility. It is far outside the valid range on purpose, so it is obvious when it survives — as long as someone looks. In a mean it is not obvious at all: it drags the result towards a large negative number that still prints as a number.

Why values above one are legitimate

Backscatter coefficients compare the power returned by a surface with that of an ideal diffuse reflector of the same area. A corner reflector directs almost all the energy back to the sensor rather than scattering it, so the ratio exceeds one. Positive decibels over cities, ships and metal structures are expected; over a field they would not be.

Why the median is safe in either scale

The logarithm is a monotonic transformation, so it preserves order and therefore preserves the median and every percentile. Only operations that add values — means, sums, convolutions — care which scale they are in. That makes percentile statistics a convenient way to report backscatter without worrying about it.

Provider conventions and amplitude products

• Different providers, different nodata: 0, −9999 and NaN all appear; always read it from the file.
• Float32 precision is ample for backscatter; the range is not the problem, the conventions are.
• Clipping for display is fine; clipping before statistics is not.
• Negative decibels are normal — most of the land is between −25 and 0 dB.
• Sigma nought and gamma nought differ; mixing them shifts values by the incidence angle's cosine.
• Amplitude products need 20·log₁₀ and have a different noise floor.

Internal links

• How to convert SAR backscatter to decibels — the full treatment
• How to read Sentinel-1 RTC data in Python — masking on read
• A SAR image looks noisy after processing — the other symptom of bad units
• Sentinel-1 products explained — sigma nought against gamma nought
• Rasterio returns wrong values because of nodata and scaling — nodata in raster reads generally
• Speckle explained — why zonal spread is large

FAQ

What does -32768 mean in Sentinel-1 data?

It is the nodata value. Mask it before any arithmetic.

Why are some SAR values greater than one?

Corner reflectors — buildings, ships, metal — return more power than the diffuse reference. Positive decibels are legitimate.

Why is my zonal mean wrong?

Almost always because it was computed on decibels. Average power and convert once.

Should I clip bright pixels?

For display, yes. For statistics, no — they are real measurements.

Is 10·log10 or 20·log10 right?

10·log₁₀ for power, 20·log₁₀ for amplitude. Sentinel-1 RTC holds power.

Are medians affected by the scale?

No. Percentiles are preserved by the logarithm.