⚙️ How-To Radar & SAR Raster & rasterio Radar & SAR Imagery · step 9 Notebook By · Updated

How to convert SAR backscatter to decibels, and average it correctly

On this page

Problem statement

Radar backscatter comes in linear power units: small positive numbers, mostly between 0.001 and 1, with a long tail of very bright targets. Plotted as they are, they show a black image with a few white dots. Converted to decibels they become a readable range from about −25 to +5, and every threshold, colour scale and published value you will compare against is in decibels.

The conversion is one line. What goes wrong is everything done after it: averages, resampling and zonal statistics computed on decibels instead of power. Those are biased, and by an amount that depends on how speckled the data is. On a Sentinel-1 scene south of Valencia, a zonal mean computed in decibels is between 0.9 and 8.6 dB darker than the correct one, depending on the land cover.

Quick answer

Convert to decibels for display and thresholds, but do arithmetic in linear units:

import numpy as np

def to_db(linear):
    with np.errstate(divide="ignore", invalid="ignore"):
        return np.where(linear > 0, 10 * np.log10(linear), np.nan)

def from_db(db):
    return 10 ** (db / 10)

def mean_db(linear, axis=None):
    """The right way to average backscatter: mean in linear, then convert."""
    return to_db(np.nanmean(linear, axis=axis))

print(to_db(np.array([1.0, 0.1, 0.01, 0.0, -32768.0])))
print(from_db(np.array([0.0, -10.0, -20.0])))
[  0. -10. -20.  nan  nan]
[1.   0.1  0.01]

Zero and the nodata value become NaN rather than -inf or an error. A difference of 10 dB is a factor of ten in power; 3 dB is a factor of two.

Linear backscatter from 0.001 to 10 mapped to decibels from −30 to +10, with typical surfaces marked. linear power and decibels-30 dB0.001-20 dB0.01-10 dB0.10 dB110 dB10open water -20.2cropland -9.5built-up -4.3brightest 1% +5.9Medians and the 99th percentile from the 20 October 2024 scene south of Valencia spatialworkflow.io
Decibels turn four orders of magnitude of power into a readable range.

Mask, keep linear values, average and resample in power

1. Mask before you convert

Nodata (-32768 in Sentinel-1 RTC) and zeros must become NaN first. The logarithm of zero is minus infinity and of a negative number is undefined, and both silently spread through later statistics.

2. Keep a linear copy

Keep the linear array alongside the decibel one, or convert back with 10 ** (db / 10) before any arithmetic.

3. Average in linear units

Means over zones, over time and over neighbourhoods are computed on power and converted to decibels at the end. The mean of decibels is the geometric mean, which is always lower.

4. Resample in linear units

Downsampling with average resampling is an average; do it on linear values. Nearest-neighbour resampling does not average and can be done in either.

5. Subtract decibels for ratios

A ratio of two powers — post over pre, VH over VV — is a difference in decibels. That is the one arithmetic operation decibels make easier.

6. Label the units

Name files and variables with the unit: vv_db, vv_lin. The two look nothing alike numerically, which is some protection, but not enough when an average of the wrong one is 2 dB off.

How much darker a zonal mean computed in decibels is than the decibels of the linear mean, by land cover. error from averaging dB instead of powercropland−0.91 dBtree cover−1.36 dBgrassland−2.07 dBbuilt-up−4.94 dBopen water−8.57 dB05 dB10 dBThe more variable the class, the larger the error; it is never zero spatialworkflow.io
Averaging decibels darkens every class, by more where the backscatter is more variable.

Zonal means and downsampling done both ways

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 — zonal means done both ways

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

bbox = [-0.46, 39.25, -0.33, 39.45]
g = grid(bbox, "EPSG:25830")
vv = read_pass(next(iter(passes(search(bbox, "2024-10-20/2024-10-20", orbit=110)).values())), "vv", g)
cat = pystac_client.Client.open(STAC, modifier=planetary_computer.sign_inplace)
wc_item = next(cat.search(collections=["esa-worldcover"], bbox=bbox).items())
with rasterio.open(wc_item.assets["map"].href) as src, WarpedVRT(
        src, crs=g["crs"], transform=g["transform"], width=g["width"], height=g["height"],
        resampling=Resampling.nearest) as vrt:
    wc = vrt.read(1)

for code, name in {80: "open water", 40: "cropland", 30: "grassland", 10: "tree cover", 50: "built-up"}.items():
    m = (wc == code) & np.isfinite(vv)
    right, wrong = mean_db(vv[m]), np.nanmean(to_db(vv[m]))
    print(f"{name:10}  linear mean {right:6.2f} dB   mean of dB {wrong:6.2f} dB   error {wrong - right:+.2f}")
open water  linear mean -10.80 dB   mean of dB -19.37 dB   error -8.57
cropland    linear mean  -8.71 dB   mean of dB  -9.62 dB   error -0.91
grassland   linear mean  -7.69 dB   mean of dB  -9.76 dB   error -2.07
tree cover  linear mean  -7.67 dB   mean of dB  -9.03 dB   error -1.36
built-up    linear mean   1.02 dB   mean of dB  -3.92 dB   error -4.94

Every class comes out darker, from 0.9 dB for cropland to 8.6 dB for open water. The water figure is the reminder that a correct mean can still be the wrong statistic: a few very bright pixels in the class — boats, reed beds, the lagoon's edges, one of them at +33 dB — dominate its linear mean. For skewed zones report the median, which is the same whichever scale it is taken in.

Example 2 — downsampling to 40 m both ways

def block_mean(a, k):
    h, w = (a.shape[0] // k) * k, (a.shape[1] // k) * k
    return np.nanmean(a[:h, :w].reshape(h // k, k, w // k, k), axis=(1, 3))

good = to_db(block_mean(vv, 4))
bad = block_mean(to_db(vv), 4)
d = bad - good
print(f"40 m cells: mean error {np.nanmean(d):+.2f} dB, worst {np.nanmin(d):+.2f} dB")
40 m cells: mean error -0.37 dB, worst -9.37 dB

A 0.37 dB average error looks small until the worst cell: 9.4 dB, in a block mixing very bright and very dark pixels — the only way the two averages can drift that far apart.

Example 3 — the ratio is a subtraction

vh = read_pass(next(iter(passes(search(bbox, "2024-10-20/2024-10-20", orbit=110)).values())), "vh", g)
a = to_db(vh) - to_db(vv)
b = to_db(vh / vv)
print(f"largest difference between the two: {np.nanmax(np.abs(a - b)):.1e} dB (float32 rounding)")
largest difference between the two: 7.6e-06 dB (float32 rounding)

Why the mean of decibels is too low

Why decibels at all

Backscatter spans several orders of magnitude within one scene: calm water returns a thousandth of what a building corner does. A logarithmic scale compresses that into a range that fits a colour ramp and turns multiplicative effects — speckle, calibration factors, ratios — into additive ones. Thresholds are therefore quoted in decibels almost universally.

Why the mean of decibels is too low

The logarithm is concave, so by Jensen's inequality the average of logs is at most the log of the average, with equality only if every value is the same. Speckle guarantees they are not. The more variable the backscatter within the zone, the larger the gap — which is why the built-up and water classes, the most variable, show the largest errors.

Why the error is not a constant you can subtract

The bias depends on the distribution within each zone: on speckle, on how many looks the data has, and on the real texture of the surface. It differs by land cover, by date and by how much filtering came before. Correcting it after the fact is harder than averaging correctly in the first place.

Common operations on backscatter and the units to do each in. which units for which operationoperationdo it inmean, zonal mean, compositelinearaverage resampling, filteringlinearratio (post/pre, VH/VV)dB subtractionthreshold, colour scaledBmedian, percentileeither — they agreeArithmetic in power; comparison, ratio and display in decibels spatialworkflow.io
Arithmetic in linear units; comparisons, ratios and display in decibels.

Amplitude versus power, medians and integer storage

  • Amplitude vs power. Some products store amplitude; decibels are then 20·log₁₀, not 10·log₁₀.
  • Medians commute with the logarithm, so a median in decibels equals the decibels of the median.
  • Positive decibels are legitimate for corner reflectors, ships and some buildings.
  • Integer storage. Some archives store decibels scaled to integers; check the scale factor.
  • Colour ramps usually span −25 to 0 dB for VV and −30 to −5 dB for VH.
  • Standard deviations in dB and in linear units describe different things; say which.

FAQ

How do I convert Sentinel-1 backscatter to dB?

10 * log10(value) for power values, after setting nodata and zeros to NaN.

Should I average in dB?

No. Average in linear units and convert the result to decibels.

Why is my dB mean lower than expected?

The mean of decibels is the geometric mean, which is always below the arithmetic mean on speckled data.

Is the median affected?

No. The median is the same whichever scale you take it in.

How do I compute a ratio of two SAR images?

Subtract their decibel values, or divide their linear values.

What does a positive dB value mean?

More power returned than a perfect diffuse reflector would give, typical of corner reflectors like buildings and ships.