Custom Sentinel-1 Pipeline: Change Detection¶

In this example we show how to use EO-Tools to insert custom processing steps in the InSAR pipeline to add incoherent amplitude change detection maps.

The process_insar and process_slc function apply fixed processing chains. They internally call the prepare_insar and prepare_slc functions which preprocess the data in the SAR geometry (coregisration and phase corrections for InSAR, burst stitching and lookup table computation for both InSAR pairs and single SLC). The lookup tables allow the projection of each subswath from the SAR geometry to the geographic coordinate system of the DEM. These functions apply further processing like phase, coherence and amplitude computation and call geocode_and_merge_iw to reproject, merge and crop the images to the geographic coordinates.

In this notebook we show a more advanced way to set up processing pipelines which allows the user to define their own processing functions.

Instead of calling the process_insar function, we use the prepare_insar and apply some custom processing steps to the files in the SAR geometry: we compute the coherence and amplitudes of the SLC products and apply a custom change_detection function to obtain a change map. Finally, the files are projected and merged in the geographic coordinate system thanks to geocode_and_merge_iw.

In [ ]:

%load_ext autoreload
%autoreload 2
import logging
logging.basicConfig(level=logging.INFO)
import matplotlib.pyplot as plt

import geopandas as gpd
from eodag import EODataAccessGateway
import rioxarray as riox
import numpy as np

# credentials need to be stored in the following file (see EODAG docs)
confpath = "/data/eodag_config.yml"
dag = EODataAccessGateway(user_conf_file_path=confpath)
# make sure cop_dataspace will be used 
dag.set_preferred_provider("cop_dataspace")
log = logging.getLogger(__name__)

%load_ext autoreload %autoreload 2 import logging logging.basicConfig(level=logging.INFO) import matplotlib.pyplot as plt import geopandas as gpd from eodag import EODataAccessGateway import rioxarray as riox import numpy as np # credentials need to be stored in the following file (see EODAG docs) confpath = "/data/eodag_config.yml" dag = EODataAccessGateway(user_conf_file_path=confpath) # make sure cop_dataspace will be used dag.set_preferred_provider("cop_dataspace") log = logging.getLogger(__name__)

Set up parameters and output dir¶

In [2]:

# change to your custom locations
data_dir = "/data/S1"

ids = [
 "S1A_IW_SLC__1SDV_20230904T063730_20230904T063757_050174_0609E3_DAA1", 
 "S1A_IW_SLC__1SDV_20230916T063730_20230916T063757_050349_060FCD_6814"
]
primary_path = f"{data_dir}/{ids[0]}.zip"
secondary_path = f"{data_dir}/{ids[1]}.zip"
output_dir="/data/res/test-change-detection-pipeline"

# change to your custom locations data_dir = "/data/S1" ids = [ "S1A_IW_SLC__1SDV_20230904T063730_20230904T063757_050174_0609E3_DAA1", "S1A_IW_SLC__1SDV_20230916T063730_20230916T063757_050349_060FCD_6814" ] primary_path = f"{data_dir}/{ids[0]}.zip" secondary_path = f"{data_dir}/{ids[1]}.zip" output_dir="/data/res/test-change-detection-pipeline"

Download S-1 products¶

In [ ]:

# load a geometry
aoi_file = "https://raw.githubusercontent.com/odhondt/eo_tools/refs/heads/main/data/Morocco_AOI.geojson"
shp = gpd.read_file(aoi_file).geometry[0]

search_criteria = {
    "productType": "S1_SAR_SLC",
    "start": "2023-09-03",
    "end": "2023-09-17",
    "geom": shp
}

results = dag.search(**search_criteria)
to_dl = [it for it in results if it.properties["id"] in ids]
print(f"{len(to_dl)} products to download")
dag.download_all(to_dl, output_dir="/data/S1/", extract=False)

# load a geometry aoi_file = "https://raw.githubusercontent.com/odhondt/eo_tools/refs/heads/main/data/Morocco_AOI.geojson" shp = gpd.read_file(aoi_file).geometry[0] search_criteria = { "productType": "S1_SAR_SLC", "start": "2023-09-03", "end": "2023-09-17", "geom": shp } results = dag.search(**search_criteria) to_dl = [it for it in results if it.properties["id"] in ids] print(f"{len(to_dl)} products to download") dag.download_all(to_dl, output_dir="/data/S1/", extract=False)

Pre-process InSAR pair¶

This function outputs two SLC (primary and secondary) image per subswath and polarization in the SAR geometry which will be further processed and projected in a geographic coodrinate system.

In this example we extract both polarimetric channels for the sake of demonstration.

In [ ]:

from eo_tools.S1.process import prepare_insar

out_dir = prepare_insar(
    prm_path=primary_path,
    sec_path=secondary_path,
    output_dir=output_dir,
    aoi_name=None,
    shp=shp,
    pol="full",
    subswaths=["IW1", "IW2", "IW3"],
    cal_type="sigma",
    apply_fast_esd=False,
    dem_upsampling=1.8,
    dem_force_download=False,
    dem_buffer_arc_sec=40,
    warp_kernel="bicubic",
)

from eo_tools.S1.process import prepare_insar out_dir = prepare_insar( prm_path=primary_path, sec_path=secondary_path, output_dir=output_dir, aoi_name=None, shp=shp, pol="full", subswaths=["IW1", "IW2", "IW3"], cal_type="sigma", apply_fast_esd=False, dem_upsampling=1.8, dem_force_download=False, dem_buffer_arc_sec=40, warp_kernel="bicubic", )

Define a simple amplitude change detection function¶

In [ ]:

def change_detection(amp_prm_file, amp_sec_file, out_file):
    log.info("Smoothing amplitudes")
    amp_prm = riox.open_rasterio(amp_prm_file)[0].rolling(x=7, y=7, center=True).mean()
    amp_sec = riox.open_rasterio(amp_sec_file)[0].rolling(x=7, y=7, center=True).mean()
    log.info("Incoherent changes")
    ch = np.log(amp_prm+1e-10) - np.log(amp_sec+1e-10)
    ch.rio.to_raster(out_file)

def change_detection(amp_prm_file, amp_sec_file, out_file): log.info("Smoothing amplitudes") amp_prm = riox.open_rasterio(amp_prm_file)[0].rolling(x=7, y=7, center=True).mean() amp_sec = riox.open_rasterio(amp_sec_file)[0].rolling(x=7, y=7, center=True).mean() log.info("Incoherent changes") ch = np.log(amp_prm+1e-10) - np.log(amp_sec+1e-10) ch.rio.to_raster(out_file)

Apply processing chains: coherence and change detection¶

Here we use two helper functions apply_to_patterns_for_pair and apply_to_patterns_for_single to avoid writing loops over polarizations and subswaths.

In [ ]:

from eo_tools.S1.process import coherence, amplitude
from eo_tools.S1.process import apply_to_patterns_for_pair, apply_to_patterns_for_single
from pathlib import Path

out_dir = f"{output_dir}/S1_InSAR_2023-09-04-063730__2023-09-16-063730/sar"
geo_dir = Path(out_dir).parent

# compute interferometric coherence
apply_to_patterns_for_pair(
    coherence,
    out_dir=out_dir,
    prm_file_prefix="slc_prm",
    sec_file_prefix="slc_sec",
    out_file_prefix="coh",
    box_size=[3, 3],
    multilook=[1, 4],
)

# compute primary amplitude
apply_to_patterns_for_single(
    amplitude,
    out_dir=out_dir,
    in_file_prefix="slc_prm",
    out_file_prefix="amp_prm",
    multilook=[2, 8],
)

# compute secondary amplitude
apply_to_patterns_for_single(
    amplitude,
    out_dir=out_dir,
    in_file_prefix="slc_sec",
    out_file_prefix="amp_sec",
    multilook=[2, 8],
)

# compute incoherent changes
apply_to_patterns_for_pair(
    change_detection,
    out_dir=out_dir,
    prm_file_prefix="amp_prm",
    sec_file_prefix="amp_sec",
    out_file_prefix="change",
)

from eo_tools.S1.process import coherence, amplitude from eo_tools.S1.process import apply_to_patterns_for_pair, apply_to_patterns_for_single from pathlib import Path out_dir = f"{output_dir}/S1_InSAR_2023-09-04-063730__2023-09-16-063730/sar" geo_dir = Path(out_dir).parent # compute interferometric coherence apply_to_patterns_for_pair( coherence, out_dir=out_dir, prm_file_prefix="slc_prm", sec_file_prefix="slc_sec", out_file_prefix="coh", box_size=[3, 3], multilook=[1, 4], ) # compute primary amplitude apply_to_patterns_for_single( amplitude, out_dir=out_dir, in_file_prefix="slc_prm", out_file_prefix="amp_prm", multilook=[2, 8], ) # compute secondary amplitude apply_to_patterns_for_single( amplitude, out_dir=out_dir, in_file_prefix="slc_sec", out_file_prefix="amp_sec", multilook=[2, 8], ) # compute incoherent changes apply_to_patterns_for_pair( change_detection, out_dir=out_dir, prm_file_prefix="amp_prm", sec_file_prefix="amp_sec", out_file_prefix="change", )

Apply geocoding, merge and crop subswaths¶

The processed images are still in the SAR geometry, we need to project them, merge the subswaths and crop to the user given geometry shp. This can be done in a single step with the geocode_and_merge_iw function.

In [ ]:

from eo_tools.S1.process import geocode_and_merge_iw
geo_dir = Path(out_dir).parent
geocode_and_merge_iw(geo_dir, shp=shp, var_names=["coh","change"])

from eo_tools.S1.process import geocode_and_merge_iw geo_dir = Path(out_dir).parent geocode_and_merge_iw(geo_dir, shp=shp, var_names=["coh","change"])

Visualize¶

In [7]:

arr_amp = riox.open_rasterio(f"{geo_dir}/change_vv.tif", masked=True)[0]
arr_amp.plot.imshow(vmin=-0.25, vmax=0.25, cmap="RdBu_r")

arr_amp = riox.open_rasterio(f"{geo_dir}/change_vv.tif", masked=True)[0] arr_amp.plot.imshow(vmin=-0.25, vmax=0.25, cmap="RdBu_r")

Out[7]:

<matplotlib.image.AxesImage at 0x7fd9f801e720>

In [8]:

arr_amp = riox.open_rasterio(f"{geo_dir}/change_vh.tif", masked=True)[0]
arr_amp.plot.imshow(vmin=-0.25, vmax=0.25, cmap="RdBu_r")

arr_amp = riox.open_rasterio(f"{geo_dir}/change_vh.tif", masked=True)[0] arr_amp.plot.imshow(vmin=-0.25, vmax=0.25, cmap="RdBu_r")

Out[8]:

<matplotlib.image.AxesImage at 0x7fd9f7856cf0>

In [9]:

arr_coh = riox.open_rasterio(f"{geo_dir}/coh_vv.tif", masked=True)[0]
arr_coh.plot.imshow(vmin=0,vmax=1, cmap="gray")

arr_coh = riox.open_rasterio(f"{geo_dir}/coh_vv.tif", masked=True)[0] arr_coh.plot.imshow(vmin=0,vmax=1, cmap="gray")

Out[9]:

<matplotlib.image.AxesImage at 0x7fd9f5f69730>

In [10]:

arr_coh = riox.open_rasterio(f"{geo_dir}/coh_vh.tif", masked=True)[0]
arr_coh.plot.imshow(vmin=0,vmax=1, cmap="gray")

arr_coh = riox.open_rasterio(f"{geo_dir}/coh_vh.tif", masked=True)[0] arr_coh.plot.imshow(vmin=0,vmax=1, cmap="gray")

Out[10]:

<matplotlib.image.AxesImage at 0x7fd9f803a0f0>