This guide has 2 routes. Pick one, or read both. Each route is complete on its own.
Optical satellites such as Sentinel-2 and Landsat see sunlight reflected from the ground, so clouds block them and night stops them. Radar sends its own microwaves, and measures what comes back. Microwaves go through cloud and rain, and they work by night. That is why radar is the first data used in a flood, when the sky is still full of cloud. Sentinel-1 is a pair of European radar satellites, with a cell of 10 m and a visit about every 6 to 12 days over most of the world.
This guide gets a Sentinel-1 image for a place in Punjab, and shows how to read it.
How to read a radar image
A radar image records backscatter: the share of the signal that returns to the satellite. It does not look like a photograph. Three things set the value of a cell.
- Smooth surfaces send the signal away from the satellite, so little returns. Calm water is smooth, so it is dark. This is why radar finds water so well.
- Rough surfaces scatter in all directions, and some returns. Crops and bare ploughed soil are in between.
- Corners return strongly. A wall and the ground in front of it make a right angle that sends the signal back, so buildings are very bright.
Wet soil returns more than dry soil, and a tall crop returns differently from a short one, so a radar image also carries information about moisture and crop growth.
Sentinel-1 sends and receives in two polarizations, which are the directions in which the wave vibrates. VV means vertical send, vertical receive. VH means vertical send, horizontal receive. VV is strong over smooth and flat surfaces and over water, and VH responds more to vegetation. Use VV to find water and VH to look at crops and trees.
The product
Sentinel-1 scenes are sold at several processing levels. Ground Range Detected, or GRD, is the basic image. It needs more steps before the values from different dates can be compared, because hills and slopes change the backscatter. Radiometrically Terrain Corrected, or RTC, has already been corrected for terrain, and it is the one to use. The Microsoft Planetary Computer serves a ready-made RTC product for the world, so you do not need to process anything.
Route A: a portal
The official archive is the Copernicus Data Space Ecosystem. The Copernicus Browser, the same site as for Sentinel-2, has a Sentinel-1 data source. You draw your area, choose the dates, and download the GRD scene, which you then process yourself, for example with the Sentinel Application Platform, SNAP. The Alaska Satellite Facility also serves Sentinel-1 scenes. Both need a free account, and their menus change, so the steps are in their own documentation. A GRD scene is large, 1 GB or more. For most uses, route B is faster.
Route B: Python, no account
The catalogue is a STAC catalogue, as in the Sentinel-2 and Landsat guides. The files are protected by a short-lived access key, which the script requests for you. You need Python with requests, numpy, rasterio, scipy and shapely.
import requests
from shapely.geometry import shape, box
# Step 1: your area, as west, south, east, north in longitude and latitude
study_area = [75.40, 30.74, 75.55, 30.84] # a box of about 14 km by 11 km near Jagraon, Punjab
# Step 2: ask the catalogue for Sentinel-1 RTC scenes of February 2023 over the box
url = "https://planetarycomputer.microsoft.com/api/stac/v1/search"
search = {"collections": ["sentinel-1-rtc"], "bbox": study_area, "datetime": "2023-02-01T00:00:00Z/2023-02-28T23:59:59Z", "limit": 50}
scenes = requests.post(url, json=search, timeout=60).json()["features"]
# Step 3: keep only the scenes that cover the whole box, and print what matters for comparing them
covering = [s for s in scenes if shape(s["geometry"]).contains(box(*study_area))]
print(len(scenes), "scenes touch the box,", len(covering), "cover all of it")
for scene in covering:
p = scene["properties"]
print(p["datetime"][:10], p["sat:orbit_state"], "relative orbit", p["sat:relative_orbit"], p["sar:polarizations"], p["platform"])
7 scenes touch the box, 7 cover all of it
2023-02-27 ascending relative orbit 100 ['VV', 'VH'] SENTINEL-1A
2023-02-23 descending relative orbit 34 ['VV', 'VH'] SENTINEL-1A
2023-02-22 ascending relative orbit 27 ['VV', 'VH'] SENTINEL-1A
2023-02-15 ascending relative orbit 100 ['VV', 'VH'] SENTINEL-1A
2023-02-11 descending relative orbit 34 ['VV', 'VH'] SENTINEL-1A
2023-02-10 ascending relative orbit 27 ['VV', 'VH'] SENTINEL-1A
2023-02-03 ascending relative orbit 100 ['VV', 'VH'] SENTINEL-1A
The lines show three things that you must read before you pick a scene. The date. The orbit state, which is ascending when the satellite moves north and descending when it moves south. And the relative orbit, the number of the track on the ground. The satellite looks at the ground from a different angle on each track and in each direction, and the backscatter changes with the angle. To compare two dates, take scenes of the same orbit state and the same relative orbit. The same track at a different date gives the same viewing angle.
Pick two scenes that can be compared
Take a scene of February and a scene of August from the same track.
# Step 1: a helper that finds the scenes of one month that cover the box
def find(start, end):
query = {"collections": ["sentinel-1-rtc"], "bbox": study_area, "datetime": f"{start}T00:00:00Z/{end}T23:59:59Z", "limit": 50}
found = requests.post(url, json=query, timeout=60).json()["features"]
return [s for s in found if shape(s["geometry"]).contains(box(*study_area))]
# Step 2: choose ascending scenes of relative orbit 100 in both months, as in the list above
def pick(scenes, orbit=100, state="ascending"):
chosen = [s for s in scenes if s["properties"]["sat:relative_orbit"] == orbit and s["properties"]["sat:orbit_state"] == state]
return sorted(chosen, key=lambda s: s["properties"]["datetime"])[-1] # the latest one of the month
winter = pick(find("2023-02-01", "2023-02-28"))
monsoon = pick(find("2023-08-01", "2023-08-31"))
print("winter scene:", winter["properties"]["datetime"][:10], "| monsoon scene:", monsoon["properties"]["datetime"][:10])
print("assets:", sorted(winter["assets"].keys()))
winter scene: 2023-02-27 | monsoon scene: 2023-08-26
assets: ['rendered_preview', 'tilejson', 'vh', 'vv']
The scene has several assets. The two that matter are vv and vh, the two polarizations.
Read a window
import numpy as np
import rasterio
from rasterio.windows import from_bounds
from rasterio.warp import transform_bounds
# Step 1: the access key. It is a short piece of text that unlocks the files for a limited time.
token = requests.get("https://planetarycomputer.microsoft.com/api/sas/v1/token/sentinel-1-rtc", timeout=60).json()["token"]
def read_window(scene, band):
asset = scene["assets"][band]
with rasterio.open(asset["href"] + "?" + token) as src:
# Step 2: the box is in longitude and latitude and the file is in metres (a UTM zone), so convert the box
left, bottom, right, top = transform_bounds("EPSG:4326", src.crs, *study_area)
window = from_bounds(left, bottom, right, top, src.transform).round_offsets().round_lengths()
values = src.read(1, window=window).astype("float32")
# Step 3: cells with no data become NaN
values[values == src.nodata] = np.nan
return values, src.window_transform(window), src.crs
vv_winter, transform, crs = read_window(winter, "vv")
vh_winter, _, _ = read_window(winter, "vh")
vv_monsoon, _, _ = read_window(monsoon, "vv")
vh_monsoon, _, _ = read_window(monsoon, "vh")
print("window:", vv_winter.shape[1], "x", vv_winter.shape[0], "cells | CRS: EPSG", crs.to_epsg(), "| cells without data:", round(float(np.isnan(vv_winter).mean()) * 100, 1), "percent")
print("linear VV: lowest", round(float(np.nanmin(vv_winter)), 4), "| median", round(float(np.nanmedian(vv_winter)), 3), "| highest", round(float(np.nanmax(vv_winter)), 1))
window: 1440 x 1114 cells | CRS: EPSG 32643 | cells without data: 0.0 percent
linear VV: lowest 0.0024 | median 0.064 | highest 92.3
The numbers are linear backscatter. Most are small, below 0.2, and a few are above 1 and some far above. A few very bright cells, such as the corner of a metal roof or a tower, dominate the range. A scale like that is hard to use, which is why radar is shown in decibels.
Convert to decibels
Decibels are ten times the base 10 logarithm of the value. They turn the huge range into a handful of numbers, and the typical land cover falls between about -25 and 0.
# Step 1: linear to decibels. Cells at zero or below cannot have a logarithm, so they become NaN.
def to_db(linear):
return 10 * np.log10(np.where(linear > 0, linear, np.nan))
vv_db, vh_db = to_db(vv_winter), to_db(vh_winter)
# Step 2: the typical values. 5 percent of the cells are below the first number, and 5 percent are above the second.
print("VV in dB: median", round(float(np.nanmedian(vv_db)), 1), "| 5th and 95th percentile", np.nanpercentile(vv_db, [5, 95]).round(1).tolist())
print("VH in dB: median", round(float(np.nanmedian(vh_db)), 1), "| 5th and 95th percentile", np.nanpercentile(vh_db, [5, 95]).round(1).tolist())
VV in dB: median -12.0 | 5th and 95th percentile [-16.2, -6.0]
VH in dB: median -18.2 | 5th and 95th percentile [-22.5, -12.8]
VH is weaker than VV, because the cross-polarized return is weak, so its decibel values are lower, by about 6 dB here. A value of -12 dB in VV and -18 dB in VH is typical of farmland. Water is far lower, below about -18 dB in VV, and buildings are far higher, often above -5 dB.
Reduce the speckle
A radar image has speckle, a grain that comes from the way the waves add up in each cell. Neighbouring cells of one field have very different values, which hides the pattern. The simplest cure is to average each cell with its neighbours, on the linear values and not on the decibels.
from scipy.ndimage import uniform_filter
# Step 1: average each cell with its 5 x 5 neighbours, on the linear values. NaN cells are filled first with the median.
def smooth(linear, size=5):
filled = np.where(np.isnan(linear), np.nanmedian(linear), linear)
return uniform_filter(filled, size=size)
# Step 2: the spread of the values in decibels over the whole window, before and after smoothing
vv_smooth_db = to_db(smooth(vv_winter))
print("standard deviation of VV in dB, all cells: raw", round(float(np.nanstd(vv_db)), 2), "| smoothed", round(float(np.nanstd(vv_smooth_db)), 2))
standard deviation of VV in dB, all cells: raw 3.18 | smoothed 2.64
Smoothing lowers the spread, here from 3.18 to 2.64 dB, and lets the shapes of fields show. Part of what it removes is grain and part is real detail, and it costs resolution: a 5 by 5 average blurs the 10 m cell to 50 m. Choose the size to suit the job: a small one for narrow features such as roads, a large one for big fields. More careful filters exist in SNAP and in other packages.
Winter against monsoon
With scenes from the same track, the difference between the two dates is meaningful.
# Step 1: the change in VV and in VH from February to August, in dB, on the smoothed images
vv_change = to_db(smooth(vv_monsoon)) - to_db(smooth(vv_winter))
vh_change = to_db(smooth(vh_monsoon)) - to_db(smooth(vh_winter))
print("median change, VV:", round(float(np.nanmedian(vv_change)), 1), "dB | VH:", round(float(np.nanmedian(vh_change)), 1), "dB")
# Step 2: how many cells got much brighter, and how many much darker, in VH
print("VH more than 3 dB brighter in August:", round(float((vh_change > 3).mean()) * 100, 1), "percent of cells")
print("VH more than 3 dB darker in August:", round(float((vh_change < -3).mean()) * 100, 1), "percent of cells")
median change, VV: -0.2 dB | VH: 0.5 dB
VH more than 3 dB brighter in August: 7.9 percent of cells
VH more than 3 dB darker in August: 4.3 percent of cells
The median change over the window is small, about a fifth of a decibel in VV and half a decibel in VH. The change is in the tails: about 8 percent of the cells are more than 3 dB brighter in August and about 4 percent are more than 3 dB darker. Crops, soil moisture and ploughing all change the backscatter, so a change map shows where something changed and does not say what. Maps of this kind are the basis of crop monitoring and flood mapping: a cell that was land and is now very dark in VV has probably gone under water.
Save the result
# Step 1: write the smoothed VV and VH of February, in dB, as a two-band GeoTIFF with the position of the window
profile = {"driver": "GTiff", "height": vv_db.shape[0], "width": vv_db.shape[1], "count": 2, "dtype": "float32",
"crs": crs, "transform": transform, "nodata": np.nan, "compress": "deflate"}
with rasterio.open("s1_rtc_2023-02_vv_vh_db.tif", "w", **profile) as dst:
dst.write(to_db(smooth(vv_winter)).astype("float32"), 1)
dst.write(to_db(smooth(vh_winter)).astype("float32"), 2)
dst.set_band_description(1, "VV dB, 5 x 5 mean")
dst.set_band_description(2, "VH dB, 5 x 5 mean")
# Step 2: read it back to check
with rasterio.open("s1_rtc_2023-02_vv_vh_db.tif") as check:
print("saved:", check.count, "bands |", check.width, "x", check.height, "cells | CRS EPSG", check.crs.to_epsg())
saved: 2 bands | 1440 x 1114 cells | CRS EPSG 32643
Open the file in QGIS and use Singleband gray with a stretch from -25 to 0 dB for VV. Calm water is black, buildings are white, and the fields have the grey of their crops.
What radar cannot tell you
- A dark cell is not always water. Smooth bare ground, a tarmac road and radar shadow behind a hill are also dark. Use a threshold together with other data, such as a DEM and a slope map.
- The values depend on the angle and the direction. Compare only scenes of the same track.
- Wet soil and rain on the leaves change the backscatter. A change between two dates is not always a change in the crop.
- A decibel is a logarithm. The average of decibels is not the decibel of an average, so average on the linear values, as the smoothing step does.
Common mistakes
- Comparing scenes from different tracks, or from ascending and descending passes.
- Averaging in decibels instead of linear values.
- Reading VH and VV as if they meant the same.
- Taking the logarithm of zero or negative values, which gives errors.
- Calling every dark cell water.
- Skipping the speckle filter, or using a filter so strong that narrow features vanish.
- Using the plain GRD product for a comparison between dates, when RTC is available.