An Introduction to SAR Change Detection Methods with NISAR

This notebook introduces several commonly used SAR change detection methods that can be applied to radar image time series data. We provide both time-series and pair-wise SAR change detection workflows using NISAR GCOV data products and xarray-based raster analysis.

More specifically, this notebook covers the following SAR change detection concepts:

Time series percentile difference and variance thresholding
Time series coefficient of variation thresholding
Log ratio-based pair-wise change detection from image pairs
Threshold-based SAR change mask generation
Exporting change detection products as GeoTIFF files

Overview¶

Prerequisites
Load NISAR Data Stack
Open and Organize GCOV data
Define Helper Python Functions
Visualize the GCOV Time Series
Computation and Visualization of Time Series Metrics
Time Series Change Detection Methods
Pair-Wise Change Detection Methods
Write Change Detection Results and Metrics to GeoTIFF Files
Summary
Resources and References

1. Prerequisites¶

Prerequisite	Importance	Notes
The `isce3` software environment for this cookbook	Necessary
How to search and access NISAR data in your AOI	Necessary
Familiarity with xarray	Helpful

Rough Notebook Time Estimate: 20 minutes

2. Load NISAR Data Stack¶

If you do not have NISAR GCOV data prepared, you can prepare data via the Search and Download NISAR GCOV Data with asf_search notebook from the NISAR Cookbook.

To learn more about using NISAR data, please visit the NISAR Data User Guide. These services are provided by the Alaska Satellite Facility.

Note: This notebook requires dual-pol data¶

Select the directory holding your NISAR GCOV data

Click the Select button
Navigate to your data directory
Click the Select button
Confirm that the desired path appears in green text
Click the Change button to alter your selection

If you receive a “widget not found” message, try refreshing your browser tab

from pathlib import Path
from ipyfilechooser import FileChooser
fc = FileChooser(Path.cwd())
display(fc)

Prepare output paths and GCOV .h5 file list

data_dir = Path(fc.selected)

# create a folder for saving output figures
plot_dir = data_dir / "plots"
plot_dir.mkdir(exist_ok=True)

# write path to GCOV files
gcov_files = sorted(list(data_dir.glob("*.h5")))
gcov_files

[PosixPath('/home/jovyan/NISAR_change_detection_example/NISAR_L2_PR_GCOV_004_042_A_019_2005_DHDH_A_20251101T022728_20251101T022803_X05010_N_F_J_001.h5'),
 PosixPath('/home/jovyan/NISAR_change_detection_example/NISAR_L2_PR_GCOV_005_042_A_019_2005_DHDH_A_20251113T022729_20251113T022803_X05009_N_F_J_001.h5'),
 PosixPath('/home/jovyan/NISAR_change_detection_example/NISAR_L2_PR_GCOV_006_042_A_019_2005_DHDH_A_20251125T022729_20251125T022804_X05009_N_F_J_001.h5'),
 PosixPath('/home/jovyan/NISAR_change_detection_example/NISAR_L2_PR_GCOV_007_042_A_019_2005_DHDH_A_20251207T022730_20251207T022804_X05009_N_F_J_001.h5'),
 PosixPath('/home/jovyan/NISAR_change_detection_example/NISAR_L2_PR_GCOV_008_042_A_019_2005_DHDH_A_20251219T022730_20251219T022805_X05009_N_F_J_001.h5'),
 PosixPath('/home/jovyan/NISAR_change_detection_example/NISAR_L2_PR_GCOV_009_042_A_019_2005_DHDH_A_20251231T022731_20251231T022805_X05009_N_F_J_001.h5')]

3. Open and Organize GCOV files¶

This section opens the NISAR GCOV .h5 files and organizes data of interest into a data cube.

%%time
import xarray as xr
import rioxarray

group_path = "/science/LSAR/GCOV/grids/frequencyA" # change this to any GCOV HDF5 group you wish

datatrees = [
    xr.open_datatree(
        f,
        engine="h5netcdf",
        decode_timedelta=False,
        phony_dims="access",
        chunks="auto",
        group=group_path,
    )
    for f in gcov_files
]

CPU times: user 1.17 s, sys: 136 ms, total: 1.31 s
Wall time: 1.03 s

import re
from urllib.parse import urlparse
from datetime import datetime
from pathlib import PurePosixPath

NISAR_TS_RE = re.compile(r"_(\d{8}T\d{6})_")

def nisar_start_time_from_url(pth: str) -> datetime:
    name = pth.name
    
    m = NISAR_TS_RE.search(name)
    if not m:
        raise ValueError(f"No NISAR timestamp found in: {pth}")
    
    return datetime.strptime(m.group(1), "%Y%m%dT%H%M%S")

dts = [nisar_start_time_from_url(pth) for pth in gcov_files]

dataarrays = [
    tree.ds.assign_coords(time=dt).expand_dims(time=1)
    for dt, tree in zip(dts, datatrees)
]

dataarrays[0].time

datacube = xr.concat(dataarrays, dim="time")
projection = datacube.isel(time=0).projection.item()
datacube

Select a Spatial Subset. This notebook allows you to subset via pixel index (option 1) or by coordinate values (option 2)¶

The full GCOV data cube is large, so this example selects a smaller subset for analysis.

Subset Option 1: Subset by Pixel Index¶

%%time
polarization = "HVHV"

datacube = datacube[polarization].isel(
    xCoordinates=slice(8000, 11000),
    yCoordinates=slice(8000, 11000)
)

datacube

CPU times: user 596 μs, sys: 82 μs, total: 678 μs
Wall time: 682 μs

Subset Option 2: Subset by Coordinate Values¶

Uncomment to use this option

# # print the coordinate range from which to select values
# print(datacube.xCoordinates.min().values, datacube.xCoordinates.max().values)
# print(datacube.yCoordinates.min().values, datacube.yCoordinates.max().values)

# %%time
# polarization = "HVHV"

# datacube = datacube[polarization].sel(
#     xCoordinates=slice(600000, 700000),
#     yCoordinates=slice(3900000, 3800000)
# )

# datacube

4. Define Some Python Helper Functions for this Notebook¶

We are defining two helper functions for this notebook:

CreateGeoTiff() to write out images
timeseries_metrics() to compute various metrics from a data stack

def create_geotiff(name, array, data_type, ndv, bandnames=None, 
                   ref_image=None, geo_t=None, projection=None):
    # If it's a 2D image we fake a third dimension:
    if len(array.shape) == 2:
        array = np.array([array])
    if ref_image == None and (geo_t == None or projection == None):
        raise RuntimeWarning('ref_image or settings required.')
    if bandnames != None:
        if len(bandnames) != array.shape[0]:
            raise RuntimeError(f'Need {Array.shape[0]} bandnames. {len(bandnames)} given')
    else:
        bandnames = [f'Band {i+1}' for i in range(array.shape[0])]
    if ref_image != None:
        refimg = gdal.Open(ref_image)
        geo_t = refimg.GetGeoTransform()
        Projection = refimg.GetProjection()
    driver = gdal.GetDriverByName('GTIFF')
    array[np.isnan(array)] = ndv
    dataset = driver.Create(name, array.shape[2], array.shape[1], 
                            array.shape[0], data_type)
    dataset.SetGeoTransform(geo_t)
    dataset.SetProjection(projection)
    for i, image in enumerate(array, 1):
        dataset.GetRasterBand(i).WriteArray(image)
        dataset.GetRasterBand(i).SetNoDataValue(ndv)
        dataset.SetDescription(bandnames[i-1])
    dataset.FlushCache()
    return name

def timeseries_metrics(datacube):     
    q = datacube.quantile([0.05, 0.5, 0.95], dim="time")
    p05 = q.sel(quantile=0.05, drop=True)
    p50 = q.sel(quantile=0.50, drop=True)
    p95 = q.sel(quantile=0.95, drop=True)
    prange = p95 - p05
    
    max = datacube.max(dim="time")
    min = datacube.min(dim="time")
    range = max - min

    mean = datacube.mean(dim="time")
    var = datacube.var(dim="time")
    std = datacube.std(dim="time")
    cv = std / xr.where(mean != 0, mean, np.nan)
    
    return xr.Dataset(
        dict(
            min=min, max=max, range=range,
            mean=mean, var=var, std=std, CV=cv,
            p05=p05, p50=p50, p95=p95,
            prange=p95 - p05
        )
    )

5. Visualize the GCOV Time Series¶

This section prepares the selected GCOV backscatter stack for visualization and creates an animation to show how radar backscatter changes through time.

5.1 Convert Backscatter to dB¶

GCOV backscatter values are stored in linear power units. For visualization and interpretation, SAR backscatter is commonly converted to decibels (dB), which compresses the range of values and makes temporal differences easier to see.

%%time
import numpy as np

datacube_db = 10*np.log10(datacube)
datacube_db

CPU times: user 1.39 ms, sys: 197 μs, total: 1.59 ms
Wall time: 1.6 ms

5.2 Create a Time Series Animation¶

The animation displays each acquisition in sequence using the selected polarization, frequency, and spatial subset. This provides a quick visual check of temporal changes and helps identify areas with strong backscatter variability before applying change detection methods.

%%capture
from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.animation import FFMpegWriter, PillowWriter

vmin, vmax = (datacube_db.quantile([0.05, 0.95], skipna=True).values.astype(float))

fig, ax = plt.subplots(figsize=(6, 6))
frame_0 = datacube_db.isel(time=0).values
im = ax.imshow(frame_0, vmin=vmin, vmax=vmax, origin="upper", cmap="grey")
cb = fig.colorbar(im, ax=ax, label="γ⁰ (dB)")
ax.set_title(np.datetime_as_string(datacube_db.time.values[0], unit='s'))
ax.set_axis_off()

def animate(i):
    frame = datacube_db.isel(time=i).values
    im.set_data(frame)
    ax.set_title(np.datetime_as_string(datacube_db.time.values[i], unit='s'))
    return (im,)

ani = FuncAnimation(fig, animate, frames=datacube_db.sizes["time"], interval=300, blit=False)

plt.show()

from IPython.display import HTML

HTML(ani.to_jshtml())

Save the animation (animation.gif):

ani.save(plot_dir/f'animation.gif', writer='pillow', fps=2)

6. Computation and Visualization of Time Series Metrics¶

Once the GCOV time series stack has been constructed, we can compute a set of temporal statistics for each pixel in the stack. These metrics summarize how radar backscatter changes through time and can help identify areas with persistent variability or anomalous temporal behavior.

The metrics computed in this section include:

Mean
Median
Maximum
Minimum
Range (Maximum − Minimum)
5th Percentile
95th Percentile
Percentile Range (95th − 5th Percentile)
Variance
Standard Deviation
Coefficient of Variation (Standard Deviation / Mean)

These statistics are computed independently for every pixel across the time dimension of the GCOV stack.

%%time

# compute temporal metrics from the GCOV stack
metrics = timeseries_metrics(datacube)

# load results into memory
metrics.load()

metrics

CPU times: user 17.2 s, sys: 3.53 s, total: 20.7 s
Wall time: 15.6 s

# access the coefficient of variation raster
cv = metrics["CV"]

The coefficient of variation and variance layers can now be visualized and thresholded to identify areas exhibiting stronger temporal variability within the SAR time series.

np.min(cv)

np.max(cv)

len(cv)

3000

from matplotlib import pyplot as plt

fig, ax = plt.subplots(1,2,figsize=(9, 3))

var = metrics["var"].values.ravel()
ax[0].hist(var, bins=100, range=np.nanpercentile(var, [1,99]))
ax[0].set_title('Variance')

cv = metrics["CV"].values.ravel()
ax[1].hist(cv, bins=100, range=np.nanpercentile(cv, [1,99]))
ax[1].set_title('Coefficient of Variation')

plt.show()

# List the metrics keys you want to plot
metric_keys=['mean', 'p50', 'max', 'min', 
             'p95', 'p05', 'prange', 'var', 'std', 'CV']
fig= plt.figure(figsize=(16,40))
idx=1
for i in metric_keys:
    ax = fig.add_subplot(5,2,idx)
    vmin, vmax = np.nanpercentile(metrics[i], [1, 99])
    ax.imshow(metrics[i],vmin=vmin,vmax=vmax,cmap='inferno')
    ax.set_title(i.upper())
    ax.axis('off')
    idx+=1

7. Time-Series SAR Change Detection Methods¶

Time-series methods analyze variability across the full GCOV acquisition stack to identify pixels with persistent or statistically significant temporal changes through time. These approaches are useful for studying seasonal behavior, long-term variability, recurring disturbances, or areas that exhibit unusually dynamic radar backscatter across the observation period.

7.1 Change Detection with the Percentile Difference and the Variance Threshold Method¶

Time-series methods identify areas with strong temporal variability by calculating statistics across the SAR acquisition stack.

In this section, we use metrics including percentile range (p95 - p05), variance and standard deviation to identify pixels with significant temporal changes.

def plot_histogram_cdf(metric='std'):
    plt.rcParams.update({'font.size': 12})
    fig = plt.figure(figsize=(14, 4)) # Initialize figure with a size
    ax1 = fig.add_subplot(121)  # 121 determines: 2 rows, 2 plots, first plot
    ax2 = fig.add_subplot(122)

    h = ax1.hist(
        metrics[metric].values.ravel(), range=np.nanpercentile(metrics[metric], [1, 99]))
    ax1.xaxis.set_label_text(f'{metric}')
    ax1.set_title('Histogram')

    n, bins, patches = ax2.hist(
        metrics[metric].values.ravel(), range=np.nanpercentile(metrics[metric], [1, 99]),
        cumulative='True', density='True', histtype='step', label='Empirical')

    ax2.xaxis.set_label_text(f'{metric}')
    ax2.set_title('CDF')

    outind = np.where(n > 0.95)   
    threshind = np.min(outind)
    thresh = bins[threshind]
    ax1.axvline(thresh,color='red')
    _ = ax2.axvline(thresh,color='red')
    plt.savefig(plot_dir/f'{metric}_histogram.png',
                dpi=200, transparent='true')

Now let’s look at the 95th - 5th percentile range of radar backscatter values across the designated temporal period

plot_histogram_cdf(metric='prange')

Here we visualize the 5% of all pixels with the largest (95th - 5th percentile) difference in the time series. We will refer to the pixels (x,y) that exceed this threshold $t$ as likely change pixels (cp):

${cp}_{x,y} = P_{x,y}^{95th} - P_{x,y}^{5th} > t$

If the threshold $t$ is set to the top 5% of pixels with the largest $(p95 - p05)$ values, the resulting change mask looks like:

def plot_threshold_classifier(metric='prange', percentage_cutoff=5):
    plt.figure(figsize=(8,8))
    thresh = np.nanpercentile(metrics[metric], 100 - percentage_cutoff)
    mask = metrics[metric] < thresh # For display we prepare the inverse mask
    plt.imshow(mask, cmap='gray')
    _=plt.title(f'Threshold Classifier on {metric} > %1.8f' % thresh)
    plt.savefig(plot_dir/f'changes_{metric}.png',
            dpi=200, transparent='true')
    return np.logical_not(mask)

metric = 'prange'
masks = {metric: plot_threshold_classifier(metric=metric)}

Thresholding other metrics
Instead of applying a threshold on the 95th - 5th percentile difference data, we can also attempt to threshold other metrics. Pixels with high standard deviation are areas where radar backscatter changes strongly throughout the time series.

Here, change pixels are defined as:

${cp}_{x,y} = \sigma > t$

with $t=CDF_{\sigma} > 0.95$ , where the threshold $t$ corresponds to the top 5% of pixels with the highest standard deviation values.

plot_histogram_cdf(metric='std')

metric = 'std'
masks[metric] = plot_threshold_classifier(metric=metric)

7.2 Change Detection with the Coefficient of Variation Method¶

We can apply a threshold to the coefficient of variation (CV) image to identify pixels with strong temporal variability in the SAR time series:

CV_{x,y} = \frac{\sigma_{x,y}}{\bar{X}_{x,y}} > t

(1)

The coefficient of variation (CV) measures the ratio of temporal variability to the mean backscatter value for each pixel.

Pixels with high CV values may indicate areas experiencing stronger temporal variability or surface change throughout the acquisition period. The coefficient of variation can be especially useful because it normalizes variability relative to the mean signal level, making it easier to compare regions with different average backscatter values.

Let’s now examine the histogram and cumulative distribution function (CDF) of the CV values.

plot_histogram_cdf(metric='CV')

$t=CDF_{CV} > 0.95$

With $t$ defined as the top 5% of pixels with the highest CV values, the resulting change mask looks like:

metric = 'CV'
masks[metric] = plot_threshold_classifier(metric=metric)

8. Pair-Wise SAR Change Detection¶

Pair-wise change detection compares two individual acquisitions to identify localized temporal changes between observation dates. In radar remote sensing, this is commonly performed using ratios in the linearly scaled power domain or differences in the logarithmic dB domain. Pair-wise approaches are useful when the objective is to isolate changes associated with a specific event or time interval, such as flooding, vegetation disturbance, snow accumulation, surface deformation or infrastructure changes.

This section identifies the dataset scenes closest to two user-defined acquisition dates within the GCOV time series. They can later be used for image differencing, visualization and threshold-based change detection analysis.

from datetime import datetime

# extract dataset acquisition times
tindex = datacube.time.values

# define target dates for comparison
dates = ('2025-11-15', '2025-12-01')

# convert strings into NumPy-compatible datetime objects
dates_ = [np.datetime64(datetime.strptime(date, '%Y-%m-%d')) for date in dates]

# find the nearest acquisition index for each target date
dates_ind = [np.argmin(np.abs(date - tindex)) for date in dates_]

# print the matching acquisition dates from the dataset
print(f'Comparing image {dates_ind[0]} from {tindex[dates_ind[0]]} with {dates_ind[1]} from {tindex[dates_ind[1]]}')

Comparing image 1 from 2025-11-13T02:27:29.000000000 with 2 from 2025-11-25T02:27:29.000000000

Compute the log ratio in dB, corresponding to the difference in dB

# select the dataset scenes closest to the user-defined dates
date_0 = datacube.isel(time=dates_ind[0])
date_1 = datacube.isel(time=dates_ind[1])

# compute the logarithmic ratio in dB between the two acquisitions
ratio_dB = 10 * np.log10(date_1 / date_0)

ratio_dB

Manually Select a Threshold

thresh sets the threshold value (for example, -2 dB)
thresh_type controls whether values below (lower) or above (upper) the threshold are masked

# compute the 1st and 99th percentile values of the ratio image
low, high = (ratio_dB.quantile([0.01, 0.99]))

print(low.compute().item())
print(high.compute().item())

# define a manual threshold for change detection
thresh = -2

# identify pixels meeting the threshold condition
a = ratio_dB < thresh

# check whether any pixels meet the threshold condition
a.max().compute().item()

-6.177210426330566
3.5505395698547346

True

thresh = -2
thresh_type = 'lower'  # 'lower': mask values below thresh, 'upper': mask values above thresh

# use the absolute dB ratio to set a symmetric display range
abs_ratio_dB = np.abs(ratio_dB)
dynamic_range = abs_ratio_dB.quantile(0.99, skipna=True).compute().item()

# create a 3-panel figure: histogram, ratio image, and threshold mask
fig, axs = plt.subplots(ncols=3, nrows=1)
fig.set_size_inches(20, 4)
plt.subplots_adjust(hspace=0.4, right=0.85)

# plot histogram of ratio values
axs[0].hist(
    ratio_dB.values.ravel(),
    bins=10,
    range=np.nanpercentile(ratio_dB, [0.1, 99.9])
)

axs[0].set_xlabel('difference [dB]')
axs[0].set_title('Histogram')

# display the log-ratio image
im0 = axs[1].imshow(
    ratio_dB,
    cmap='RdBu',
    vmin=-dynamic_range,
    vmax=dynamic_range
)

cbar = fig.colorbar(im0, orientation='vertical', ax=axs.ravel().tolist(), shrink=0.7)
cbar.set_label('[dB]')
axs[1].set_title('Image')

# create a binary mask from the selected threshold
mask = (ratio_dB < thresh if thresh_type == 'lower' else ratio_dB > thresh).astype(np.int8)

# display the threshold mask
axs[2].imshow(mask, cmap='gray')
axs[2].set_title('Mask')

# convert numpy datetime64 values to date strings for the title and filename
date0_str = np.datetime_as_string(tindex[dates_ind[0]], unit='D')
date1_str = np.datetime_as_string(tindex[dates_ind[1]], unit='D')

fig.suptitle(f'{date0_str} {date1_str}')

logratiolabel = f'logratio_{date0_str}_{date1_str}'

plt.savefig(
    plot_dir / f'{logratiolabel}.png',
    dpi=200,
    transparent=True
)

# store the inverse of the mask for later use
masks[logratiolabel] = np.logical_not(mask)

Do you think the threshold is appropriate? If not, try adjusting it.

9. Write Change Detection Results and Metrics to GeoTIFF Files¶

This section exports the computed change detection results as georeferenced GeoTIFF files.

Both time-series metrics and pair-wise change detection outputs can be saved as GIS-ready raster products for visualization and downstream analysis.

9.1 Determine Output Geometry¶

First, we define the geotransform and projection information for the output GeoTIFFs.

Because the workflow operates directly from GCOV data opened with xarray, the output geometry is derived from the dataset coordinate arrays.

# Calculate pixel spacing from the coordinate arrays
xres = float(datacube.xCoordinates[1] - datacube.xCoordinates[0])
yres = float(datacube.yCoordinates[1] - datacube.yCoordinates[0])

# Build the geotransform for the output rasters
geotrans = [
    float(datacube.xCoordinates.min()),   # upper-left x coordinate
    xres,                # pixel width
    0,
    float(datacube.yCoordinates.max()),   # upper-left y coordinate
    0,
    -abs(yres)           # negative pixel height for north-up rasters
]

geotrans

[690650.0, 20.0, 0, 3851110.0, 0, -20.0]

9.2 Export Time-Series Metrics¶

This section exports the time-series metric rasters, such as variance and coefficient of variation, as GeoTIFF files for GIS visualization and analysis.analysis.

# Create a directory for the exported metric images
metricDir = plot_dir / "ts_metrics"
metricDir.mkdir(exist_ok=True)

print(metricDir)

/home/jovyan/NISAR_change_detection_example/plots/ts_metrics

from osgeo import gdal, osr

srs = osr.SpatialReference()
srs.ImportFromEPSG(projection)
proj = srs.ExportToWkt()

/home/jovyan/NISAR_GCOV_Cookbook/.pixi/envs/isce3/lib/python3.13/site-packages/osgeo/osr.py:424: FutureWarning: Neither osr.UseExceptions() nor osr.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default.
  warnings.warn(

labeldb = polarization

names = []  # keep track of exported filenames

for metric in metrics:

    name_ = metricDir / f'{metric}_{labeldb}.tif'

    create_geotiff(
        str(name_),
        metrics[metric],
        gdal.GDT_Float32,
        np.nan,
        [metric],
        geo_t=geotrans,
        projection=proj
    )

    names.append(str(name_))

9.3 Build a Virtual Raster Stack (VRT)¶

A Virtual Raster (VRT) can be created to combine the exported metric GeoTIFFs into a single virtual raster stack.

The VRT references the original rasters without duplicating the underlying data.

# Build a virtual raster stack (VRT) from the exported metric GeoTIFFs
cmd = f'gdalbuildvrt -separate -overwrite -vrtnodata nan {str(metricDir)}/{labeldb}.vrt ' + ' '.join(names)

# Print the GDAL command for reference
print(cmd)

# Execute the GDAL VRT-building command
!{cmd}

gdalbuildvrt -separate -overwrite -vrtnodata nan /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/HVHV.vrt /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/min_HVHV.tif /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/max_HVHV.tif /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/range_HVHV.tif /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/mean_HVHV.tif /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/var_HVHV.tif /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/std_HVHV.tif /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/CV_HVHV.tif /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/p05_HVHV.tif /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/p50_HVHV.tif /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/p95_HVHV.tif /home/jovyan/NISAR_change_detection_example/plots/ts_metrics/prange_HVHV.tif
0...10...20...30...40...50...60...70...80...90...100 - done.

9.4 Export Pair-Wise Change Detection Masks¶

This section writes the threshold-based pair-wise change detection masks to GeoTIFF files. These rasters identify pixels that met the selected change threshold between the two acquisitions and preserve the spatial reference information needed for GIS visualization and downstream analysis.

# Write each threshold mask as a GeoTIFF
for metric in masks:

    # Define the output filename for the threshold mask
    fnmetric = plot_dir / f'{metric}_thresholds.tif'

    print(fnmetric)

    # Export the binary mask as a georeferenced GeoTIFF
    create_geotiff(
        str(fnmetric),
        masks[metric],
        gdal.GDT_Byte,
        np.nan,
        [metric],
        geo_t=geotrans,
        projection=proj
    )

/home/jovyan/NISAR_change_detection_example/plots/prange_thresholds.tif
/home/jovyan/NISAR_change_detection_example/plots/std_thresholds.tif
/home/jovyan/NISAR_change_detection_example/plots/CV_thresholds.tif
/home/jovyan/NISAR_change_detection_example/plots/logratio_2025-11-13_2025-11-25_thresholds.tif

10. Summary¶

This notebook demonstrated several common SAR change detection methods using a NISAR GCOV time series opened and analyzed with xarray.

The workflow included:

loading and subsetting NISAR data
visualizing temporal changes
calculating time-series statistics
applying threshold-based change detection methods
performing a pair-wise analysis
exporting the resulting products as georeferenced GeoTIFF files

11. Resources and references¶

Author: Alex Lewandowski, Julia White

Adapted for NISAR data from earlier work by:
- Franz J Meyer; University of Alaska Fairbanks