--- title: "Sentinel-1 SAR flood mapping in Python: the 2022 Pakistan floods" description: "Map floods with free Sentinel-1 SAR in Python: dB conversion, Lee speckle filter, Otsu threshold and flood area, using the 2022 Pakistan floods." url: https://www.technicalgis.com/tutorials/sentinel-1-sar-flood-mapping/ author: TechnicalGIS Team publisher: TechnicalGIS (https://www.technicalgis.com/about/) published: 2026-10-15 updated: 2026-10-15 language: en-GB --- # Sentinel-1 SAR flood mapping in Python: the 2022 Pakistan floods > Map floods with free Sentinel-1 SAR in Python: dB conversion, Lee speckle filter, Otsu threshold and flood area, using the 2022 Pakistan floods. - Level: Intermediate - Tested with: Python 3.12 - Source code: [Download scripts](https://www.technicalgis.com/code/sentinel-1-sar-flood-mapping/) ## Key takeaways 1. Calm water is a mirror to radar, so it looks dark in Sentinel-1 images, through cloud, day or night. A threshold on VV backscatter in decibels separates water from land. 2. Otsu's method found a threshold of −14.1 dB from the flood image itself; with a pre-flood image as a baseline, 1,933 km² of the 8,095 km² study area was newly flooded on 8 September 2022. 3. Average in linear power, not in decibels, and filter speckle (a Lee filter here) before thresholding; then clean the mask with a morphological opening and a minimum patch size. 4. A 12-day time series from one orbit shows the flood rising from late August, peaking on 8 September and still covering 332 km² of new ground on Christmas Day. ## Before you start - Comfortable with NumPy arrays and rasterio (the Sentinel-2 Lake Mead tutorial introduces STAC and WarpedVRT) - Python 3.10+ with pystac-client, planetary-computer, rasterio, numpy, scipy, scikit-image, pandas, geopandas and matplotlib - An internet connection: imagery is streamed from the Microsoft Planetary Computer, with no account or API key **Short answer:** to map a flood with Sentinel-1 in Python, load a radar image from before the flood and one from during it, convert the VV backscatter to decibels (`10 * log10`), reduce speckle with a filter such as Lee, and choose a water threshold automatically with Otsu’s method (`skimage.filters.threshold_otsu`). Pixels below the threshold are open water. Pixels that are water now but were not before are flood. Multiply the count by the pixel area to get the flooded area. In August and September 2022, record monsoon rain and glacier melt put about a third of Pakistan under water. Over 1,700 people died and more than 30 million were affected. Cloud covered Sindh for weeks, which made optical satellites such as Sentinel-2 almost useless. Radar sees through cloud. In this tutorial we use free **Sentinel-1** data to map the flood around Manchar Lake and the town of Dadu, and we follow it from June to December. You will learn how to: 1. search the Microsoft Planetary Computer’s **Sentinel-1 RTC** collection with `pystac-client`, and why you should stick to one orbit; 2. convert backscatter to **decibels**, and why you average in linear power; 3. reduce **speckle** with a Lee filter written in a few lines of NumPy and SciPy; 4. pick a water threshold automatically with **Otsu’s method**; 5. turn the result into **flood extent maps, areas, a GeoTIFF and a GeoPackage**, and a **time series** for the whole monsoon. All outputs below are real results from [`sar_flood_mapping.py`](https://www.technicalgis.com/code/sentinel-1-sar-flood-mapping/), run on 15 October 2026 with Python 3.12, pystac-client 0.9.0 and rasterio 1.5.2. ## Why radar for floods? Sentinel-1 carries a C-band **synthetic aperture radar (SAR)**. It sends microwave pulses obliquely towards the ground and records how much energy comes back: the **backscatter**. Rough surfaces such as fields, towns and trees scatter energy in all directions, so some returns to the satellite and they look bright. **Calm open water** acts like a mirror and reflects the pulse away, so it looks very dark. Because the radar makes its own illumination at a wavelength of 5.6 cm, it works at night and through cloud and rain. We use the **VV** polarisation (vertical transmit, vertical receive), which gives the clearest contrast between water and land. Sentinel-1 also records VH, which helps with vegetation. ## Step 1: search for Sentinel-1 RTC scenes with STAC The Planetary Computer’s `sentinel-1-rtc` collection holds Sentinel-1 GRD scenes that have already been **radiometrically terrain corrected (RTC)**: geocoded to a map grid and corrected for slope brightness, as Cloud Optimised GeoTIFFs. Doing that yourself from raw GRD data is a job for ESA SNAP or `sarsen`. Different orbits look at the ground from different angles, which changes the backscatter. So we keep **one relative orbit** (42, ascending) and compare like with like: `sar_flood_mapping.py` ```python STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1" BBOX = (67.45, 26.15, 68.25, 27.05) # lon/lat: Manchar Lake, Dadu and the Indus in Sindh, Pakistan SEASON = "2022-06-01/2022-12-31" ORBIT = 42 # one relative orbit = identical viewing geometry on every date BEFORE, AFTER = "2022-06-16", "2022-09-08" RES = 40 # analysis resolution in metres CRS = "EPSG:32642" # UTM zone 42N, the native CRS of these scenes catalog = pystac_client.Client.open(STAC_URL, modifier=planetary_computer.sign_inplace) search = catalog.search(collections=["sentinel-1-rtc"], bbox=BBOX, datetime=SEASON, query={"sat:relative_orbit": {"eq": ORBIT}, "sat:orbit_state": {"eq": "ascending"}}) items = sorted(search.items(), key=lambda i: i.datetime) ``` `planetary_computer.sign_inplace` adds a free, short-lived access token to every asset URL; no account is needed. `Output — STAC search` ```text 18 scenes from relative orbit 42 (ascending), 04 Jun to 25 Dec 2022 S1A_IW_GRDH_1SDV_20220908T133610_20220908T133635_044914_055D60_rtc assets: vh, vv, tilejson, rendered_preview polarisations: ['VV', 'VH'] | mode IW | pixel 10 m | ENL 4.4 | EPSG:32642 vv raster:bands: [{'nodata': -32768, 'data_type': 'float32', 'spatial_resolution': 10}] ``` That is one scene every **12 days**. Sentinel-1B failed in December 2021, so in 2022 only Sentinel-1A was flying. Since Sentinel-1C began operations in 2025, the revisit is back to 6 days. The `sar:looks_equivalent_number` property (ENL 4.4) tells us how much the speckle has already been averaged down. We use it in the filter below. ## Step 2: understand decibels and speckle The RTC values are **linear power** (γ⁰). Water is around 0.001–0.01 and a bright building can exceed 10. That range is too wide to see or threshold, so we convert to **decibels**: ```python def to_db(power: np.ndarray) -> np.ndarray: with np.errstate(divide="ignore", invalid="ignore"): return 10 * np.log10(np.where(power > 0, power, np.nan)) ``` Every SAR image also has **speckle**, a grainy texture from interference between the many scatterers in each pixel. It is not sensor noise you can calibrate away; it is part of how coherent imaging works. The **Lee filter** (Lee, 1980) smooths homogeneous areas while keeping edges. It compares the local variance with the variance expected from speckle alone (1 / ENL), and smooths only where the two are similar: ```python def lee_filter(img: np.ndarray, size: int = 7, enl: float = 4.4) -> np.ndarray: """Lee (1980) speckle filter on linear power. enl = equivalent number of looks.""" nan = np.isnan(img) img = np.where(nan, np.nanmedian(img), img) mean = ndimage.uniform_filter(img, size) var = np.maximum(ndimage.uniform_filter(img * img, size) - mean * mean, 0) ci2 = var / np.maximum(mean * mean, 1e-12) # local coefficient of variation, squared w = np.clip(1 - (1 / enl) / np.maximum(ci2, 1e-12), 0, 1) out = mean + w * (img - mean) # smooth flat areas, keep edges out[nan] = np.nan return out ``` To see what it does, we read a 6 km chip around Dadu at the full 10 m resolution, using a rasterio window so only that part of the COG is downloaded: ```python CHIP = (67.75, 26.70, 67.81, 26.76) # 6 km box around Dadu town and the floodwater surrounding it cx, cy = to_utm.transform([CHIP[0], CHIP[2]], [CHIP[1], CHIP[3]]) with rasterio.open(example.assets["vv"].href) as src: chip = src.read(1, window=from_bounds(cx[0], cy[0], cx[1], cy[1], src.transform)).astype("float32") chip[chip == -32768] = np.nan chip_lee = lee_filter(chip, size=7, enl=ENL) def local_std(db: np.ndarray, size: int = 5) -> float: """Median of the 5 x 5 moving-window standard deviation: pixel-to-pixel noise, not scene contrast.""" m = ndimage.uniform_filter(db, size) return float(np.nanmedian(np.sqrt(np.maximum(ndimage.uniform_filter(db * db, size) - m * m, 0)))) print(f"pixel-to-pixel noise (median 5x5 std): raw {local_std(to_db(chip)):.2f} dB, Lee {local_std(to_db(chip_lee)):.2f} dB") ``` `Output — speckle (10 m chip)` ```text chip 603 x 659 pixels; linear gamma0 0.0009..30.54 pixel-to-pixel noise (median 5x5 std): raw 1.84 dB, Lee 0.65 dB ``` The filter cuts the pixel-to-pixel variation from 1.84 dB to 0.65 dB, while the edges of the town, fields and floodwater stay sharp. ![Two greyscale radar images of a 6 km area around Dadu on 8 September 2022. On the left, the raw image is grainy; on the right, the Lee-filtered version is smoother but keeps the sharp edge between the bright town and the black floodwater to its east.](https://www.technicalgis.com/_astro/02-speckle-lee-filter.DP69uAxg_Z2v3qMi.webp) _Raw VV backscatter at 10 m (left) and after a 7 × 7 Lee filter (right). Dadu town is the bright area; the black area to the east is floodwater, and the sharp, straight edge between them is a man-made line such as a road or embankment._ ## Step 3: read every date onto one 40 m grid As in the [Sentinel-2 Lake Mead tutorial](https://www.technicalgis.com/tutorials/sentinel-2-stac-lake-mead/), rasterio’s `WarpedVRT` streams each COG straight onto a fixed grid, reading from the internal overviews so only a fraction of the data is downloaded. Our study area is 81 × 100 km, and 40 m pixels are plenty for a flood of this size. ```python def read_gamma0(href: str) -> np.ndarray: """Warp a remote RTC COG onto the analysis grid (averaging linear power), NaN where empty.""" with rasterio.open(href) as src, WarpedVRT(src, crs=CRS, transform=TRANSFORM, width=WIDTH, height=HEIGHT, resampling=Resampling.average, src_nodata=src.nodata, nodata=np.nan, dtype="float32") as vrt: return vrt.read(1) vv_db = {} for it in items: d = it.datetime.strftime("%Y-%m-%d") power = read_gamma0(it.assets["vv"].href) # Each 40 m pixel averages 16 pixels (~4 independent resolution cells), so ENL rises ~4x. vv_db[d] = to_db(lee_filter(power, size=5, enl=ENL * 4)) ``` ![Three panels of the 81 by 100 km study area. Left: the pre-flood radar image of 16 June 2022, mostly mid-grey, with Manchar Lake as a dark patch. Middle: 8 September 2022, with huge black areas along the Indus and west of Dadu. Right: a colour composite where the flooded areas are bright red.](https://www.technicalgis.com/_astro/01-vv-before-after.BIGSVZKP_Z1f85hb.webp) _VV backscatter before (16 June) and during (8 September) the flood, and a change composite with red = before and green and blue = after. Anything that was bright before and dark after shows as red, which here is the flood._ The composite on the right is a quick, threshold-free way to see change: areas that stayed the same are grey, and areas that went dark are red. ## Step 4: choose the water threshold with Otsu’s method The histogram of the flood image has two peaks: water near −20 dB and land near −8 dB. **Otsu’s method** finds the threshold that best separates them, by maximising the variance between the two classes: ```python before, after = vv_db[BEFORE], vv_db[AFTER] valid = np.isfinite(before) & np.isfinite(after) thr = threshold_otsu(after[valid]) ``` `Output — Otsu threshold` ```text VV threshold from 2022-09-08: -14.1 dB median VV: water-side -19.8 dB, land-side -8.3 dB ``` ![Overlaid histograms of VV backscatter in decibels. The pre-flood histogram in grey has one peak near −10 dB. The flood histogram in blue has a second peak near −20 dB. A vertical orange line marks the Otsu threshold at −14.1 dB, in the valley between the peaks.](https://www.technicalgis.com/_astro/03-otsu-histogram.Bm6G-bsC_14a9Lu.webp) _Histograms of Lee-filtered VV backscatter before and during the flood. The flood adds a second, dark peak; Otsu places the threshold at −14.1 dB in the valley between the two._ The threshold comes from the flood image, where water is plentiful. We apply it to every date, so all dates are classified the same way. Commonly published thresholds for VV water are around −15 to −18 dB, so −14.1 dB is in the expected range for a scene with this much water. ## Step 5: build a clean water mask and the flood extent A raw threshold leaves single dark pixels everywhere: radar shadow, smooth roads and tarmac. A **morphological opening** removes isolated pixels, and a minimum patch size of 7 pixels (about 1 ha) removes small speckles: ```python def water_mask(db: np.ndarray) -> np.ndarray: w = np.isfinite(db) & (db < thr) w = ndimage.binary_opening(w, iterations=1) # drop isolated dark pixels lab, n = ndimage.label(w) sizes = ndimage.sum(w, lab, index=np.arange(1, n + 1)) return np.isin(lab, 1 + np.flatnonzero(sizes >= 7)) # keep patches >= ~1 ha water_before, water_after = water_mask(before) & valid, water_mask(after) & valid flood = water_after & ~water_before ``` **Flood = water now that was not water before.** Comparing with a pre-flood image removes permanent water, such as Manchar Lake and the Indus channel. It also removes surfaces that are always dark, such as smooth sand or tarmac. `Output — flood extent` ```text valid area: 8,095 km² water on 2022-06-16: 958 km² water on 2022-09-08: 2,470 km² newly flooded: 1,933 km² (23.9% of the area) water before but not after: 421 km² pixels at least 3 dB darker: 2,432 km²; of the Otsu flood pixels, 97% dropped by 3 dB or more ``` ![Map of the study area over a greyscale radar background. Manchar Lake and the Indus are blue (permanent water). A huge orange flood area surrounds the lake and stretches north along both sides of the Indus, encircling the labelled town of Dadu.](https://www.technicalgis.com/_astro/04-flood-extent-map.BPVQCM2L_Z1pfMLY.webp) _Flood extent on 8 September 2022. Blue is water on both dates (Manchar Lake and the Indus); orange is newly flooded land, 1,933 km², nearly a quarter of the study area. Dadu town itself stayed dry._ Two checks give confidence in the result: - **A threshold-free cross-check.** 97% of the flood pixels became at least 3 dB darker, which is the classic change-detection criterion. - **Lost water is plausible.** 421 km² was water in June but not in September. Some of it is seasonal water that drained or was replaced by growing crops. Some is probably flooded vegetation, which looks bright rather than dark (see the limitations below). Finally, save the result as a classified GeoTIFF and as vector polygons for use in QGIS or a web map: ```python shapes = rasterio.features.shapes(flood.astype("uint8"), mask=flood, transform=TRANSFORM) polys = gpd.GeoDataFrame.from_features([{"geometry": g, "properties": {"class": "flood"}} for g, _ in shapes], crs=CRS) polys["area_km2"] = polys.area / 1e6 polys.to_file(OUT / "flood_extent_2022-09-08.gpkg", driver="GPKG") ``` `Output — vector output` ```text 8,320 flood polygons, total 1,933 km², largest 1,208 km² ``` ## Step 6: follow the flood through time Because every scene is on the same grid, a time series is just a loop: ```python rows = [] for d in dates: v = np.isfinite(vv_db[d]) & valid w = water_mask(vv_db[d]) & v rows.append({"date": pd.Timestamp(d), "water_km2": round(w.sum() * PIXEL_KM2), "flooded_km2": round((w & ~water_before).sum() * PIXEL_KM2), "valid_pct": round(v.mean() * 100, 1)}) ts = pd.DataFrame(rows) ``` `Output — water area through the 2022 monsoon` ```text date water_km2 flooded_km2 valid_pct 2022-06-04 925 89 100.0 2022-06-16 958 0 100.0 2022-06-28 1065 162 100.0 2022-07-10 767 150 100.0 2022-07-22 1012 374 100.0 2022-08-03 1191 562 100.0 2022-08-15 912 490 100.0 2022-08-27 2156 1608 100.0 2022-09-08 2470 1933 100.0 2022-09-20 1807 1369 100.0 2022-10-02 1428 1022 100.0 2022-10-14 1281 872 100.0 2022-10-26 1112 713 100.0 2022-11-07 985 586 100.0 2022-11-19 874 482 100.0 2022-12-01 816 419 100.0 2022-12-13 768 366 100.0 2022-12-25 758 332 100.0 peak: 08 September 2022, 1,933 km² newly flooded ``` ![Line chart of water area from June to December 2022. Total open water hovers around 800–1,200 km² in June to mid-August, jumps to about 2,150 km² on 27 August and peaks at 2,470 km² on 8 September, then declines steadily to about 760 km² by late December. The newly flooded line peaks at 1,933 km².](https://www.technicalgis.com/_astro/05-flood-area-timeseries.CRP6m7kW_RwCW9.webp) _Open water (blue area) and newly flooded land (orange line) on all 18 dates. The flood arrives in late August, peaks on 8 September and drains slowly: on 25 December, 332 km² of normally dry land was still under water._ The curve tells the story of the disaster. Rain in July and early August flooded fields locally. The main flood wave arrived in late August, as the Indus and hill torrents from the west filled the plain. In early September, the authorities breached the embankment of Manchar Lake to protect the towns of Sehwan and Bhan Saeedabad, and the water drained only slowly through the autumn. Months later, land that should have been sown with wheat was still under water. ## Limitations and how to improve them - **Dark is not always water.** Smooth tarmac, dry sand and radar shadow behind hills are dark too. Change detection removes most of these; a **HAND** (height above nearest drainage) or slope mask from the [Copernicus DEM](https://www.technicalgis.com/data-sources/#src-copernicus-dem-glo-30-glo-90) removes the rest. - **Wind roughens water.** Strong wind can make a lake look bright and be missed. - **Flooded vegetation and towns look bright, not dark.** Radar bounces between the water surface and stems or walls (**double bounce**), so flooded crops, forest and streets are under-mapped by a darkness threshold. VH and the VV/VH ratio help. - **One threshold for a huge area** can be biased. Use tiled Otsu, or validate against a Copernicus Emergency Management Service map of the same event. ## Frequently asked questions ### Can Sentinel-1 see through clouds? Yes. Sentinel-1 is a C-band radar with a 5.6 cm wavelength, which passes through cloud, haze and smoke, and it makes its own illumination, so it works at night too. Only very heavy rain cells cause minor effects. That is why it is the standard satellite for flood mapping. ### What dB threshold should I use for water in Sentinel-1? There is no universal value. Open water in VV is usually below about −15 to −18 dB, and in VH below about −22 to −25 dB, but it depends on wind, incidence angle and processing. Compute the threshold from your own image with Otsu’s method or a tiled variant, as in this tutorial, rather than copying a number. ### What is the difference between Sentinel-1 GRD and RTC? **GRD** (Ground Range Detected) is ESA’s standard product: calibrated backscatter in the radar’s own ground-range geometry, projected onto an Earth ellipsoid without terrain correction. **RTC** (radiometrically terrain corrected) is GRD that has also been geocoded with a DEM and corrected for the brightening and darkening of slopes. For flood mapping in hilly areas, use RTC; on a flat plain like Sindh, GRD is very similar. ### Why is my SAR image so grainy? That is speckle, a feature of every coherent radar image. Reduce it with a speckle filter (Lee, refined Lee, Gamma-MAP), by averaging pixels (multilooking), or by averaging several dates. Filter the linear power values, not the decibels. ### Is Sentinel-1 data free? Yes. All Sentinel data is free and open under the Copernicus licence. You can download it from the [Copernicus Data Space Ecosystem](https://www.technicalgis.com/data-sources/#src-copernicus-data-space-ecosystem), from ASF, or stream analysis-ready RTC from the [Microsoft Planetary Computer](https://www.technicalgis.com/data-sources/#src-microsoft-planetary-computer) as in this tutorial. ## Next steps - Add VH and map flooded vegetation with the VV/VH ratio. - Mask steep and high ground with HAND computed from the Copernicus DEM; [Terrain analysis in Python with rasterio and rioxarray](https://www.technicalgis.com/tutorials/raster-terrain-analysis-rasterio/) shows how to stream and process it. - Compare with optical water indices in [Measuring Lake Mead’s shrinking shoreline with Sentinel-2 and STAC](https://www.technicalgis.com/tutorials/sentinel-2-stac-lake-mead/), or combine both on cloud-free days. - Go beyond water extent to water quality in [Monitoring reservoir water quality from space with Sentinel-2](https://www.technicalgis.com/tutorials/reservoir-water-quality-sentinel-2/).- Map elevation in 3D with lidar in [LiDAR point clouds in Python with laspy](https://www.technicalgis.com/tutorials/lidar-point-clouds-laspy/). - Overlay the flood polygons on settlements and roads from [OSMnx](https://www.technicalgis.com/tutorials/street-networks-osmnx/) to estimate exposure. ## Data and licences - **Sentinel-1 RTC**: contains modified Copernicus Sentinel data 2022, processed by ESA; radiometric terrain correction by Catalyst/Microsoft. Licensed CC BY 4.0. Accessed from the [Microsoft Planetary Computer](https://planetarycomputer.microsoft.com/dataset/sentinel-1-rtc) on 15 October 2026. - Reference: Lee, J.-S. (1980), “Digital image enhancement and noise filtering by use of local statistics”, _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 2(2); Otsu, N. (1979), _IEEE Transactions on Systems, Man, and Cybernetics_ 9(1). - Code: [`sar_flood_mapping.py`](https://www.technicalgis.com/code/sentinel-1-sar-flood-mapping/), tested with Python 3.12.0, pystac-client 0.9.0, planetary-computer 1.0.0, rasterio 1.5.2, scipy 1.18.1 and scikit-image 0.26.0. ## Key terms - **SAR**: Synthetic aperture radar: an active microwave sensor that images the ground through cloud and at night by measuring the echo of its own pulses. - **Backscatter (γ⁰)**: The share of radar energy returned towards the satellite. Gamma-nought is backscatter normalised for the illuminated area; it is usually shown in decibels. - **RTC**: Radiometric terrain correction: processing that removes the brightness effects of slopes facing towards or away from the radar. - **Speckle**: The grainy, salt-and-pepper texture of every coherent radar image, caused by interference between many scatterers in each pixel. - **Otsu's method**: An automatic threshold that splits a histogram into two classes by maximising the variance between them; it works best on a bimodal histogram. ## Data for this tutorial - [Microsoft Planetary Computer](https://www.technicalgis.com/data-sources/#src-microsoft-planetary-computer) · Sentinel-1 RTC as Cloud Optimised GeoTIFFs via STAC - [Copernicus Data Space Ecosystem](https://www.technicalgis.com/data-sources/#src-copernicus-data-space-ecosystem) · Sentinel-1 GRD and SLC from ESA - [Copernicus DEM (GLO-30 / GLO-90)](https://www.technicalgis.com/data-sources/#src-copernicus-dem-glo-30-glo-90) · Useful for HAND and slope masks Licences and formats for each source are listed in the [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). --- Source: https://www.technicalgis.com/tutorials/sentinel-1-sar-flood-mapping/