Map Flood Areas Using Sentinel-1 SAR Imagery

This example uses:

This example shows how to map flooded areas using Sentinel-1 Synthetic Aperture Radar (SAR) images.

This example requires the Image Processing Toolbox™ Hyperspectral Imaging Library. You can install the Image Processing Toolbox Hyperspectral Imaging Library from Add-On Explorer. For more information about installing add-ons, see Get and Manage Add-Ons. The Image Processing Toolbox Hyperspectral Imaging Library requires desktop MATLAB®, as MATLAB® Online™ and MATLAB® Mobile™ do not support the library.

In this example, you perform these tasks to map flooded areas using Sentinel-1 SAR images.

Read SAR images of the same region taken before and during a flood.
Select a region of interest (ROI) in the images.
Perform preprocessing operations on the SAR images.
Identify the flooded areas by taking the difference of the preprocessed SAR images of the ROI.
Project the identified flooded areas on the map.

Download the Data Set

This example uses SAR images from the Copernicus Open Access Hub [1], which show the flood that occurred in Western Germany in the middle of July, 2021. The images focus on the Euskirchen district in North Rhine-Westphalia, an area on the west side of the Rhine river heavily impacted by the flood. The first image is the archive image, acquired on July 3, before the flood. The second image is the crisis image, acquired on July 15, during the flood.

Download the SAR data using the helperDownloadDataset helper function. The helper function is attached to this example as a supporting file.

dataDir = fullfile(pwd,"data");
helperDownloadDataset(dataDir)

Started downloading the datasets(can take several minutes)...
Completed.

The helper function downloads and unzips two folders, corresponding to the two instances of the acquired SAR data. The SAR data for the archive image is in the folder S1A_IW_GRDH_1SDV_20210703T055051_20210703T055116_038609_048E45_3ADC.SAFE. The SAR data for the crisis image is in the folder S1A_IW_GRDH_1SDV_20210715T055052_20210715T055117_038784_049389_2074.SAFE. In this example, you use only the vertical copolarized (VV) images from the data, because the cross-polarized bands do not provide any additional information for flood mapping.

Read Sentinel-1 SAR Data from Archive Image

Read the meta information for the SAR data of the archive image, acquired before the flood, by parsing the manifest.safe file in the data folder using the parseSafeFile helper function. The helper function is attached to this example as a supporting file.

archiveDataDir = fullfile(dataDir,"S1A_IW_GRDH_1SDV_20210703T055051_20210703T055116_038609_048E45_3ADC.SAFE");
metaDataArchive = parseSafeFile(archiveDataDir)

metaDataArchive = struct with fields:
            rootDir: "data\S1A_IW_GRDH_1SDV_20210703T055051_20210703T055116_038609_048E45_3ADC.SAFE"
           platform: 'Sentinel-1A'
          orbitPass: 'descending'
        productType: 'GRD'
            acqMode: 'IW'
           numSwath: 1
            swathID: {'IW'}
    numPolarization: 2
       polarization: {'VV'  'VH'}
       acqStartTime: '2021-07-03T05:50:51.668638'
        acqStopTime: '2021-07-03T05:51:16.667578'
        swathCoords: [4×2 double]
           dataFile: {'measurement/s1a-iw-grd-vv-20210703t055051-20210703t055116-038609-048e45-001.tiff'  'measurement/s1a-iw-grd-vh-20210703t055051-20210703t055116-038609-048e45-002.tiff'}
        prodAnnFile: {'annotation/s1a-iw-grd-vv-20210703t055051-20210703t055116-038609-048e45-001.xml'  'annotation/s1a-iw-grd-vh-20210703t055051-20210703t055116-038609-048e45-002.xml'}
         calAnnFile: {'annotation/calibration/calibration-s1a-iw-grd-vv-20210703t055051-20210703t055116-038609-048e45-001.xml'  'annotation/calibration/calibration-s1a-iw-grd-vh-20210703t055051-20210703t055116-038609-048e45-002.xml'}
       noiseAnnFile: {'annotation/calibration/noise-s1a-iw-grd-vv-20210703t055051-20210703t055116-038609-048e45-001.xml'  'annotation/calibration/noise-s1a-iw-grd-vh-20210703t055051-20210703t055116-038609-048e45-002.xml'}

Read the image data from the SAR data using the readS1Image helper function. The helper function is attached to this example as a supporting file. Read the VV image from the image data.

imgDataArchive = readS1Image(metaDataArchive);
imgArchive = imgDataArchive.iw_vv.data;

Get the meta information for the annotation by parsing the annotation file using the parseAnn helper function. The helper function is attached to this example as a supporting file.

annFilesArchive = fullfile(metaDataArchive.rootDir, metaDataArchive.prodAnnFile);
annMetaInfoArchive = parseAnn(annFilesArchive(1,1))

annMetaInfoArchive = struct with fields:
                adsHeader: [1×1 struct]
              qualityData: [1×1 struct]
        generalAnnotation: [1×1 struct]
          imageAnnotation: [1×1 struct]
    geolocationGridPoints: [1×210 struct]

The meta information from the annotation file shows the preprocessing performed on the image.

processingInfoArchive = annMetaInfoArchive.qualityData.processingInformation

processingInfoArchive = struct with fields:
                      rawDataAnalysisUsed: 'true'
                        orbitDataFileUsed: 'false'
                     attitudeDataFileUsed: 'false'
             rxVariationCorrectionApplied: 'true'
           antennaElevationPatternApplied: 'true'
             antennaAzimuthPatternApplied: 'true'
      antennaAzimuthElementPatternApplied: 'true'
                                 dcMethod: 'Data Analysis'
                              dcInputData: 'Range Compressed'
    rangeSpreadingLossCompensationApplied: 'true'
                    srgrConversionApplied: 'true'
                       detectionPerformed: 'true'
          thermalNoiseCorrectionPerformed: 'false'

Select ROI of Image

Select a relevant ROI showing the river region. The rows represent the range direction and the columns represent the azimuth direction.

rowIdx = 1801:4801;
colIdx = 1:1501;
roiImgArchive = imgArchive(rowIdx,colIdx);
figure
image(roiImgArchive)
colormap(gray)
xticks([])
yticks([])
axis equal tight
title("ROI of Archive Image")

Preprocess Archive Image

Perform these preprocessing steps on the ROI from the archive image.

Thermal noise removal
Radiometric calibration
Despeckling
Conversion to decibel (dB) Scale
Thresholding in decibel (dB) scale.

Remove thermal noise present in the archive image using the removeThermalNoise helper function. The helper function is attached to this example as a supporting file.

noiseFile = fullfile(metaDataArchive.rootDir,metaDataArchive.noiseAnnFile{1});
imgTnr = removeThermalNoise(imgArchive,noiseFile);

Apply radiometric calibration to the denoised image to transform the pixel values from integer digital numbers (DN) to backscatter coefficients (sigma nought) using the calibrateImage helper function. The helper function is attached to this example as a supporting file.

calFile = fullfile(metaDataArchive.rootDir,metaDataArchive.calAnnFile{1});
imgTnrCal = calibrateImage(imgTnr,calFile);

Remove the granular noise (speckles) from the calibrated image by using a low-pass adaptive Wiener filter.

imgTnrCalSpk = wiener2(imgTnrCal);

Convert the despeckled image to the decibel (dB) scale. Plot the histogram of the ROI of the logarithmic image.

imgTnrCalSpkdB = 10*log10(imgTnrCalSpk);
roiImgTnrCalSpkdB = imgTnrCalSpkdB(rowIdx,colIdx);
figure
histogram(roiImgTnrCalSpkdB)

Using the histogram as reference, remove everything except the water from the image by thresholding the image with a threshold value of -16 obtained from the histogram to have only the water bodies in the processed image.

imgProc1 = roiImgTnrCalSpkdB;
imgProc1(imgProc1 >= -16) = 1;

Display the ROI of the archive image after each stage of preprocessing.

figure
tiledlayout(1,4)
nexttile
roiImgTnr = imgTnr(rowIdx,colIdx);
roiImgTnr = rescale(roiImgTnr);
imshow(roiImgTnr,[0 0.1])
title({"Thermal Noise", ...
    "Removal"})
nexttile
roiImgTnrCal = imgTnrCal(rowIdx,colIdx);
imshow(roiImgTnrCal)
title({"Radiometric", ...
    "Calibration"})
nexttile
roiImgTnrCalSpk = imgTnrCalSpk(rowIdx,colIdx);
imshow(roiImgTnrCalSpk)
title("Despeckling")
nexttile
imshow(imgProc1)
title("Thresholding")

Read Sentinel-1 SAR Data from Crisis Image

Read the meta information for the SAR data of the crisis image, acquired during the flood, by parsing the manifest.safe file in the data folder using the parseSafeFile helper function.

crisisDataDir = fullfile(dataDir,"S1A_IW_GRDH_1SDV_20210715T055052_20210715T055117_038784_049389_2074.SAFE");
metaDataCrisis = parseSafeFile(crisisDataDir)

metaDataCrisis = struct with fields:
            rootDir: "data\S1A_IW_GRDH_1SDV_20210715T055052_20210715T055117_038784_049389_2074.SAFE"
           platform: 'Sentinel-1A'
          orbitPass: 'descending'
        productType: 'GRD'
            acqMode: 'IW'
           numSwath: 1
            swathID: {'IW'}
    numPolarization: 2
       polarization: {'VV'  'VH'}
       acqStartTime: '2021-07-15T05:50:52.398820'
        acqStopTime: '2021-07-15T05:51:17.396178'
        swathCoords: [4×2 double]
           dataFile: {'measurement/s1a-iw-grd-vv-20210715t055052-20210715t055117-038784-049389-001.tiff'  'measurement/s1a-iw-grd-vh-20210715t055052-20210715t055117-038784-049389-002.tiff'}
        prodAnnFile: {'annotation/s1a-iw-grd-vv-20210715t055052-20210715t055117-038784-049389-001.xml'  'annotation/s1a-iw-grd-vh-20210715t055052-20210715t055117-038784-049389-002.xml'}
         calAnnFile: {'annotation/calibration/calibration-s1a-iw-grd-vv-20210715t055052-20210715t055117-038784-049389-001.xml'  'annotation/calibration/calibration-s1a-iw-grd-vh-20210715t055052-20210715t055117-038784-049389-002.xml'}
       noiseAnnFile: {'annotation/calibration/noise-s1a-iw-grd-vv-20210715t055052-20210715t055117-038784-049389-001.xml'  'annotation/calibration/noise-s1a-iw-grd-vh-20210715t055052-20210715t055117-038784-049389-002.xml'}

Read the image data from the SAR data using the readS1Image helper function. Read the VV image from the image data.

imageDataCrisis= readS1Image(metaDataCrisis);
imgCrisis = imageDataCrisis.iw_vv.data;

Get the meta information for the annotation by parsing the annotation file using the parseAnn helper function.

annFilesCrisis = fullfile(metaDataCrisis.rootDir,metaDataCrisis.prodAnnFile);
annMetaInfoCrisis = parseAnn(annFilesCrisis(1,1))

annMetaInfoCrisis = struct with fields:
                adsHeader: [1×1 struct]
              qualityData: [1×1 struct]
        generalAnnotation: [1×1 struct]
          imageAnnotation: [1×1 struct]
    geolocationGridPoints: [1×210 struct]

The meta information from the annotation file shows the preprocessing performed on the image.

processingInfoCrisis = annMetaInfoCrisis.qualityData.processingInformation

processingInfoCrisis = struct with fields:
                      rawDataAnalysisUsed: 'true'
                        orbitDataFileUsed: 'false'
                     attitudeDataFileUsed: 'false'
             rxVariationCorrectionApplied: 'true'
           antennaElevationPatternApplied: 'true'
             antennaAzimuthPatternApplied: 'true'
      antennaAzimuthElementPatternApplied: 'true'
                                 dcMethod: 'Data Analysis'
                              dcInputData: 'Range Compressed'
    rangeSpreadingLossCompensationApplied: 'true'
                    srgrConversionApplied: 'true'
                       detectionPerformed: 'true'
          thermalNoiseCorrectionPerformed: 'false'

Select ROI of Image

Select an ROI in the crisis image that corresponds to the ROI of the archive image.

roiImgCrisis = imgCrisis(rowIdx,colIdx);
figure
image(roiImgCrisis)
colormap(gray)
xticks([])
yticks([])
axis equal tight
title("ROI of Crisis Image")

Preprocess Crisis Image

Preprocess the crisis image using the same steps used in the preprocessing of the archive image.

noiseFile = fullfile(metaDataCrisis.rootDir,metaDataCrisis.noiseAnnFile{1});
imgTnr = removeThermalNoise(imgCrisis,noiseFile);
calFile = fullfile(metaDataCrisis.rootDir,metaDataCrisis.calAnnFile{1});
imgTnrCal = calibrateImage(imgTnr,calFile);
imgTnrCalSpk = wiener2(imgTnrCal);
imgTnrCalSpkdB = 10*log10(imgTnrCalSpk);
roiImgTnrCalSpkdB = imgTnrCalSpkdB(rowIdx,colIdx);
imgProc2 = roiImgTnrCalSpkdB;
imgProc2(imgProc2 >= -16) = 1;

Display the ROI of the crisis image after each stage of preprocessing.

figure
tiledlayout(1,4)
nexttile
roiImgTnr = imgTnr(rowIdx,colIdx);
roiImgTnr = rescale(roiImgTnr);
imshow(roiImgTnr,[0 0.1]) 
title({"Thermal Noise", ...
    "Removal"})
nexttile
roiImgTnrCal = imgTnrCal(rowIdx,colIdx);
imshow(roiImgTnrCal)
title({"Radiometric", ...
    "Calibration"})
nexttile
roiImgTnrCalSpk = imgTnrCalSpk(rowIdx,colIdx);
imshow(roiImgTnrCalSpk)
title("Despeckling")
nexttile
imshow(imgProc2)
title("Thresholding")

Identify Flooded Areas

Identify the flooded areas by subtracting the thresholded archive image from the thresholded crisis image.

floodedAreas = imgProc2 - imgProc1;

Threshold the difference to get a flood area mask.

floodAreaMask = floodedAreas;
floodAreaMask(floodAreaMask==0) = 1;

Overlay the flood area mask on the ROI of the archive image to show the flooded areas.

roiImgArchive = rescale(roiImgArchive);
floodAreaMaskInv = imcomplement(floodAreaMask);
overlayImg = imoverlay(roiImgArchive,floodAreaMaskInv,"red");

Display the flood area mask and the flooded areas.

figure
tiledlayout(1,2)
nexttile
imshow(floodAreaMask)
title("Flood Area Mask")
nexttile
imshow(overlayImg)
title("Flooded Areas")

Project Flood Area Mask on Map

Display a scatter plot using the latitude and longitude coordinates taken from the geolocationGridPoints field of the meta information to visualize the location that the images represent.

gdpoints = annMetaInfoArchive.geolocationGridPoints;
gdpoints = struct2table(gdpoints);
figure
geoscatter(gdpoints,"latitude","longitude")

Fit a similarity geometric transformation to the control point pairs using the fitgeotform2d function. A similarity geometric transformation estimates a transformation that can correct isotropic scaling, rotation, and translation distortions between the images. Specify the projected CRS as WGS 84 / UTM Zone 31N, which has the EPSG authority code 32631.

p = projcrs(32631);
[x,y] = projfwd(p,gdpoints.latitude,gdpoints.longitude);
movingPoints = [x y];
line = gdpoints.line;
pixel = gdpoints.pixel;
fixedPoints = [line pixel];
tform = fitgeotform2d(movingPoints,fixedPoints,"similarity");

Apply the transformation to the image.

floodAreaMaskTransformed = imwarp(floodAreaMaskInv,tform);

Create a spatial reference object using the latitude and longitude information of the ROI obtained from the European Space Agency (ESA).

latROI = [50.7391 50.7906 50.6564 50.6048 50.7391];
lonROI = [7.1766 6.7590 6.7155 7.1323 7.1766];
[xb,yb] = projfwd(p,latROI,lonROI);
xlimits = [min(xb) max(xb)];
ylimits = [min(yb) max(yb)];
R = maprefcells(xlimits,ylimits,size(floodAreaMaskTransformed),RowsStartFrom="east");
R.ProjectedCRS = p;

Represent the geographic extent of the ROI using a mappolyshape (Mapping Toolbox) object.

dataRegion = mappolyshape(xb,yb);
dataRegion.ProjectedCRS = p;

Obtain the x and y world coordinates of the pixels of the flood area in the polyX and polyY cell arrays using the intrinsicToWorld (Mapping Toolbox) function. Apply a threshold of 16 on the transformed mask to consider only the flooded area pixels.

[rowsVec,colsVec] = find(floodAreaMaskTransformed > 16);
[polyX,polyY] = intrinsicToWorld(R,colsVec,rowsVec);

Create a mappointshape (Mapping Toolbox) object using the polyX and polyY cell arrays containing the world coordinates for the flood area pixels.

shape = mappointshape(polyX,polyY);
shape.ProjectedCRS = p;

Visualize the flood area in a geographic axes with the data region.

figure
geoplot(dataRegion,LineWidth=2)
hold on
geoplot(shape,"r.")
geobasemap satellite

References

[1] Copernicus Open Access Hub https://scihub.copernicus.eu/.

[2] Selmi, Luigi. "Flood Mapping Using the Sentinel-1 Imagery and the ESA SNAP S1 Toolbox." Version 1.0, October 5, 2021. https://github.com/luigiselmi/luigiselmi.github.io/blob/master/assets/sentinel-1/flood_mapping_using_sentinel-1_imagery_v1.pdf.

Map Flood Areas Using Sentinel-1 SAR Imagery

Download the Data Set

Read Sentinel-1 SAR Data from Archive Image

Select ROI of Image

Preprocess Archive Image

Read Sentinel-1 SAR Data from Crisis Image

Select ROI of Image

Preprocess Crisis Image

Identify Flooded Areas

Project Flood Area Mask on Map

References

Related Topics