import numpy as np
import intake
import json
import matplotlib.pyplot as plt
import holoviews as hv
from holoviews.streams import RangeXY
from scipy.ndimage import uniform_filter
from functools import partial
"bokeh") hv.extension(
This notebook will provide an empirical demonstration of speckle - how it originates, how it visually and statistically looks like, and some of the most common approaches to filter it.
Speckle is defined as a kind of noise that affects all radar images. Given the multiple scattering contributions originating from the various elementary objects present within a resolution cell, the resulting backscatter signal can be described as a random constructive and destructive interference of wavelets. As a consequence, speckle is the reason why a granular pattern normally affects SAR images, making it more challenging to interpret and analyze them.
Credits: ESRI
Let’s make an example of a cornfield (with a typical backscattering value of about -10 dB). According to the following equation:
\[ \sigma^0 = \frac{1}{\text{area}} \sum_{n \in \text{area}} \sigma_n \]
We should ideally have a uniform discrete sigma naught \(\sigma^0\) value, given that the cornfield pixel is the only individual contributor.
However, since we already learned from the previous notebooks that a pixel’s ground size can be in the order of tens of meters (i.e., 10 meters for Sentinel-1), we can imagine that different distributed targets in the scene contribute to the global backscattered information.
Let´s replicate this behavior with an ideal uniform area constituted by 100 pixels and then by adding 30% of speckle.
= -10 # in dB, a typical value for cornfields
ideal_backscatter = 12
width = (width, width)
size = np.full(size, ideal_backscatter)
ideal_data = 10 ** (
ideal_data_linear / 10
ideal_data # Convert dB to linear scale for speckle addition
)
= 0.3
speckle_fraction = int(
num_speckled_pixels 0] * size[1] * speckle_fraction
size[# Rayleigh speckle noise
) = np.random.choice(
speckled_indices * width, num_speckled_pixels, replace=False
width # random indices for speckle
)
# Initialize speckled data as the same as the ideal data
= ideal_data_linear.copy()
speckled_data_linear
= np.random.gumbel(scale=1.0, size=num_speckled_pixels)
speckle_noise *= (
speckled_data_linear.ravel()[speckled_indices] # Add speckle to the selected pixels
speckle_noise
)
= 10 * np.log10(ideal_data_linear)
ideal_data_dB = 10 * np.log10(speckled_data_linear)
speckled_data_dB =(16, 10))
plt.figure(figsize
# Ideal data
2, 2, 1)
plt.subplot(="gray", vmin=-20, vmax=0)
plt.imshow(ideal_data_dB, cmap"Ideal Backscatter (Cornfield)")
plt.title(="Backscatter (dB)")
plt.colorbar(label
# Speckled data
2, 2, 2)
plt.subplot(="gray", vmin=-20, vmax=0)
plt.imshow(speckled_data_dB, cmapf"Speckled Backscatter ({int(speckle_fraction * 100)}% of Pixels)")
plt.title(="Backscatter (dB)")
plt.colorbar(label
= 25
bins = np.histogram(ideal_data_dB.ravel(), bins=bins,
hist_ideal, bins_ideal range=(-20, 0))
= np.histogram(
hist_speckled, bins_speckled =bins, range=(-20, 0)
speckled_data_dB.ravel(), bins
)= max(
max_freq max(), hist_speckled.max()
hist_ideal.# maximum frequency for normalization
)
# Histogram for ideal data
2, 2, 3)
plt.subplot(=bins, range=(-20, 0), color="gray",
plt.hist(ideal_data_dB.ravel(), bins=0.7)
alpha0, max_freq)
plt.ylim("Histogram of Ideal Backscatter")
plt.title("Backscatter (dB)")
plt.xlabel("Frequency")
plt.ylabel(
# Histogram for speckled data
2, 2, 4)
plt.subplot(=bins, range=(-20, 0), color="gray",
plt.hist(speckled_data_dB.ravel(), bins=0.7)
alpha0, max_freq)
plt.ylim(
plt.title(f"Histogram of Speckled Backscatter ({int(speckle_fraction * 100)}%)"
)"Backscatter (dB)")
plt.xlabel("Frequency")
plt.ylabel(
plt.tight_layout()
/tmp/ipykernel_3697/3970694764.py:26: RuntimeWarning: invalid value encountered in log10
speckled_data_dB = 10 * np.log10(speckled_data_linear)
Figure 1: Synthetic data that emulates speckles in microwave backscattering
We can imagine that the second plot represents a real SAR acquisition over a cornfield, while the first plot represents an ideal uniform SAR image over a cornfield land (no speckle). The introduction of a simulated 30% speckle noise could be related to the presence of distributed scatterers of any sort present in the scene, which would cause a pixel-to-pixel variation in terms of intensity.
All the random contributions (such as the wind) would result in a different speckle pattern each time a SAR scene is acquired over the same area. Many subpixel contributors build up a complex scattered pattern in any SAR image, making it erroneous to rely on a single pixel intensity for making reliable image analysis. In order to enhance the degree of usability of a SAR image, several techniques have been put in place to mitigate speckle. We will now show two of the most common approaches: the temporal and the spatial filter.
7.1 Lake Neusiedl data
We load a dataset that contains the CORINE land cover and Sentinel-1 \(\sigma^0_E\) at a 20 meter resolution. This is the same data presented in notebook 6.
= "https://huggingface.co/datasets/martinschobben/microwave-remote-sensing/resolve/main/microwave-remote-sensing.yml"
url = intake.open_catalog(url)
cat = cat.speckle.read().compute()
fused_ds fused_ds
<xarray.Dataset> Size: 108MB Dimensions: (y: 1221, x: 1230, time: 8) Coordinates: spatial_ref int64 8B 0 * time (time) datetime64[ns] 64B 2023-08-05T16:51:22 ... 2023-10-28... * x (x) float64 10kB 5.282e+06 5.282e+06 ... 5.294e+06 5.294e+06 * y (y) float64 10kB 1.571e+06 1.571e+06 ... 1.559e+06 1.559e+06 Data variables: land_cover (y, x) float64 12MB 12.0 12.0 12.0 12.0 ... 41.0 41.0 41.0 41.0 sig0 (time, y, x) float64 96MB -6.99 -7.32 -8.78 ... -14.32 -14.22
We also create the same dashboard for backscatter of different landcover types over time. In order to make this code reusable and adaptable we will define the following function plot_variability
, which allows the injection of a spatial and/or temporal filter. It is not important to understand all the code of the following cell!
# Load encoding
with cat.corine_cmap.read()[0] as f:
= json.load(f)
color_mapping_data
# Get mapping
= {item["value"]: item for item in
color_mapping "land_cover"]}
color_mapping_data[
# Get landcover codes present in the image
= np.unique(
present_landcover_codes ~np.isnan(fused_ds.land_cover.values)].
fused_ds.land_cover.values[int)
astype(
)
def load_image(var_ds, time, land_cover, x_range, y_range,
=None):
filter_fun_spatial"""
Callback Function Landcover.
Parameters
----------
time: panda.datetime
time slice
landcover: int
land cover type
x_range: array_like
longitude range
y_range: array_like
latitude range
Returns
-------
holoviews.Image
"""
if time is not None:
= var_ds.sel(time=time)
var_ds
if land_cover == "\xa0\xa0\xa0 Complete Land Cover":
= var_ds.sig0
sig0_selected_ds else:
= int(land_cover.split()[0])
land_cover_value = var_ds.land_cover == land_cover_value
mask_ds = var_ds.sig0.where(mask_ds)
sig0_selected_ds
if filter_fun_spatial is not None:
= filter_fun_spatial(sig0_selected_ds.values)
sig0_np else:
= sig0_selected_ds.values
sig0_np
# Convert unfiltered data into Holoviews Image
= hv.Dataset(
img "x"], sig0_selected_ds["y"], sig0_np),
(sig0_selected_ds["x", "y"],
["sig0"
)
if x_range and y_range:
= img.select(x=x_range, y=y_range)
img
return hv.Image(img)
def plot_variability(var_ds, filter_fun_spatial=None,
=None):
filter_fun_temporal
= var_ds.sig0.quantile(0.02).item()
robust_min = var_ds.sig0.quantile(0.98).item()
robust_max
= [
bin_edges + j * 0.5
i for i in range(int(robust_min) - 2, int(robust_max) + 2)
for j in range(2)
]
= {"\xa0\xa0\xa0 Complete Land Cover": 1}
land_cover
land_cover.update(
{f"{int(value): 02} {color_mapping[value]['label']}": int(value)
for value in present_landcover_codes
}
)= var_ds.sig0["time"].values
time
= RangeXY()
rangexy
if filter_fun_temporal is not None:
= filter_fun_temporal(var_ds)
var_ds = partial(
load_image_ =var_ds, filter_fun_spatial=filter_fun_spatial,
load_image, var_ds=None
time
)= (
dmap =["Landcover"], streams=[rangexy])
hv.DynamicMap(load_image_, kdims=land_cover)
.redim.values(Landcover=True, bins=bin_edges)
.hist(normed
)
else:
= partial(
load_image_ =var_ds, filter_fun_spatial=filter_fun_spatial
load_image, var_ds
)= (
dmap =["Time", "Landcover"],
hv.DynamicMap(load_image_, kdims=[rangexy])
streams=time, Landcover=land_cover)
.redim.values(Time=True, bins=bin_edges)
.hist(normed
)
= hv.opts.Image(
image_opts ="Greys_r",
cmap=True,
colorbar=["hover"],
tools=(robust_min, robust_max),
clim="equal",
aspect=False,
framewise=500,
frame_height=500,
frame_width
)
= hv.opts.Histogram(width=350, height=555)
hist_opts
return dmap.opts(image_opts, hist_opts)
Now, lets work on the real-life dataset to see how speckle actually looks like.
plot_variability(fused_ds)