EO Training

DAMAGE ASSESSMENT WITH SENTINEL-1 & SENTINEL-2

Categories: Fire, Security,

DAMAGE ASSESSMENT WITH SENTINEL-1 & SENTINEL-2

Introduction Print

On 4 August 2020 in Lebanon, at around 17:55 local time (14:55 UTC), a fire broke out in Warehouse 12 at the Port of Beirut, where a cargo of 2,750 tonnes of ammonium nitrate was stored, followed by 2 explosions a few minutes later. The first explosion produced a large cloud of smoke, damaging heavily the warehouse, while the second that was caused by nitrogen dioxide, a byproduct of ammonium nitrate decomposition, produced a red-orange cloud surrounded by a white condensation cloud. The consequences were severe in many fields and levels. It is estimated that more than 200 people lost their lives, 6,500 were injured and 300,000 people went homeless since homes as far as 10 km away were damaged. The blast affected over half of Beirut, with the cost to be above $15 billion in property damages.

Within the port area, the explosion destroyed a section of shoreline leaving a crater around 124m in diameter and 43m in depth. Experts estimated that the explosion was one of the biggest non-nuclear explosions in history, it was detected as a seismic event of magnitude 3.3. It was felt in many countries such as Turkey, Syria, Israel, Palestine and parts of Europe, and was heard in Cyprus, more than 240 km away. The Lebanese government declared a two-week state of emergency in response to the disaster.

Note: This exercise was developed by the Serco EO Training Team within the fraamework of the RUS contract. The final outcome of this tutorial has not been validated and the information contained in this tutorial does not purport to constitute professional advice.

Data used
  • Three Sentinel-1 IW SLC images acquired on 25 July 2020, 31 July 2020 and 6 August 2020 [downloadable at https://scihub.copernicus.eu/]
    S1B_IW_SLC__1SDV_20200725T033445_20200725T033512_022622_02AEEE_A0BA
    S1A_IW_SLC__1SDV_20200731T033531_20200731T033558_033693_03E7B2_D1A9
    S1B_IW_SLC__1SDV_20200806T033445_20200806T033512_022797_02B43D_3A0
    B
  • Two Sentinel-2 Level 2A images (Tile ID: T36SYC) acquired on 24 July 2020 and 18 August 2020 [downloadable at https://scihub.copernicus.eu/]
    S2A_MSIL2A_20200724T081611_N0214_R121_T36SYC_20200724T110803
    S2B_MSIL2A_20200818T081609_N0214_R121_T36SYC_20200818T110855
Software used

Internet browser, SNAP + Sentinel-1 and Sentinel-2 Toolboxes, QGIS.

NOTE: This tutorial was prepared on a Linux Ubuntu 16.04 OS, steps may differ for other OS (Windows version under development)

Reference folder structure

Folder structure referenced in the exercise is as follows:

  • /HAZA08_LebanonDamageAssessment
    • /AuxData – includes auxiliary data needed for the processing – provided below or you can re-create them yourself
    • /Original – should contain your downloaded input products
    • /Processing – should contain all intermediate and final results of the processing organized in subfolders

Data download –ESA SciHUB

In this step, we will download the Sentinel-1 and Sentinel-2 scenes we will use for the exercise, from the Copernicus Open Access Hub using the online interface. Launch your Internet Browser and go to https://scihub.copernicus.eu

Go to “Open HUB”, if you do not have an account please register by going to “Sign-up” in the LOGIN menu in the upper right corner.

If you had to fill-in the registration form, you will receive an activation link by e-mail. Once your account is activated or if you already have an account, “LOGIN”.
Navigate over Lebanon, eastern Mediterranean Sea (approximate areagreen rectangle).

We need to download three Sentinel-1 images over the area of interest, two before the event and one after it. Zoom in a bit more, switch to “drawing mode” and draw a search rectangle approximately as indicated below. Open the search menu by clicking to the left part of the search bar. We will specify at once the parameters for both products before the date of the explosion and for the one afterwards as well:

For Sentinel-1:
Sensing period: From 2020/07/25 to 2020/08/06
Select: Mission: Sentinel-1
Product Type: SLC
Sensor Mode: IW
Relative Orbit Number: 21

Then click on the “Search” icon. The search returns 3 results for the time period we set. We need to download them all, so we click on the “Download Product” icon of each one:
S1B_IW_SLC__1SDV_20200725T033445_20200725T033512_022622_02AEEE_A0BA
S1A_IW_SLC__1SDV_20200731T033531_20200731T033558_033693_03E7B2_D1A9
S1B_IW_SLC__1SDV_20200806T033445_20200806T033512_022797_02B43D_3A0
B

NOTE 1: Please keep in mind that you cannot download more than 2 products at the same time, per account from SciHub.

Return to the search menu and set the parameters for the two Sentinel-2 products we will need, one before and one after the explosion. Deselect Sentinel-1 mission and apply the following steps.

For Sentinel-2:
Sensing period: From 2020/07/24 to 2020/08/18
Select: Mission: Sentinel-2
Product Type: S2MSI2A

Then click on the “Search” icon. The search returns 6 results for the time period and the area we set. If you draw a larger rectangle, you will have more results when searching. Download only the two following images by clicking on the “Download Product” icon:
S2A_MSIL2A_20200724T081611_N0214_R121_T36SYC_20200724T110803
S2B_MSIL2A_20200818T081609_N0214_R121_T36SYC_20200818T110855

The products will be downloaded in the Download folder of your browser. Move the downloaded files to the exercise folder …/HAZA08_LebanonDamageAssessment/Original

SNAP –open and explore data

Open SNAP software from the icon located on the desktop or go to Applications →Processing →SNAP Desktop. Click the Open Product icon , navigate to: …/HAZA08_LebanonDamageAssessment/Original folder and open all Sentinel-1 products (from oldest to the most recent). Alternatively, drag the products from the folder one by one and drop them to the Product Explorer Window (first the 20200725, then the 20200731 and last the 20200806).The opened products will appear in Product Explorer window.

Click + to expand the contents of product [1] from 25 July 2020, then expand Bands folder and double click on Intensity_IW1_VV band to visualize it in the View window. You can go to the World View tab and zoom in to see the location of the opened product on the globe (See NOTE 2).

NOTE 2: The Interferometric Wide (IW) swath mode captures three sub-swaths using Terrain Observation with Progressive Scans SAR (TOPSAR). IW SLC products contain one image per sub-swath and one per polarisation channel, for a total of three (single polarisation) or six (dual polarisation) images in an IW product. Each sub-swath image consists of a series of bursts, where each burst has been processed as a separate SLC image. The images for all bursts in all sub-swaths are resampled to a common pixel spacing grid in range and azimuth while preserving the phase information. (Source: https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/acquisition-modes/interferometric-wide-swath) The Beirut port is located on the IW3 sub-swath of the selected images.

Open the Intensity_IW2_VV and Intensity_IW3_VV bands as well, and go to Window →Tile Horizontally. Go to the Navigation tab and click on the two icons shown within the red rectangular below,to synchronize the views and the cursor position between the views. (See NOTE 3).

NOTE 3: The RADAR instrument onboard Sentinel-1 carries an antenna that is looking always to the right during its pass. These two scenes were acquired during descending pass (the satellite was moving in direction from north to south). That is why we see that the view of the image appears as if “mirrored”, because the view shows the pixels in order of the data acquisition.

Sentinel-1 Processing–Part 1

We will use the GraphBuilder tool, to create a chain with the steps of the processes we want to apply and at the end,only the final product will be physically saved (this way we will also save disk space since the products of the intermediate steps will not be stored). In this part, we will reach up to the step where we will use two pairs of Sentinel-1 images, in order to create two coregistered products. The first pair will consist of the two images before the event, and the second pair will consist of the one image right before the event and the other one after it.

Graph Builder

Go to Tools → GraphBuilder to build our graph. We can see that the graph has only two operators: Read (to read the input) and Write (to write the output). Below there also are the corresponding to the operators’ tabs. Right-click on the Write operator and Delete it. The corresponding tab will be removed as well. This is to avoid confusion to the sequence of the graph. The Write operator will be added again at the end.

For now, we will not define any parameters in the tabs, we will first create the graph and then we will go tab by tab to insert the desired parameters. Follow the instructions written in black to build the graph and once it is ready, go back and insert the parameters written in purple. There will be a reminder and further instructions when you reach that step.

Read the input products

There is only one Read operator in the graph (for the first product of the pair) and we will add one more, for the second product of the pair. Right-click at the empty white space below the Read operator and go to Add → Input-Output → Read. The Read (2) operator will be created. Place it below Read.

In the Read tab, as Source Product, select the product with Name: S1B_IW_SLC__1SDV_20200725T033445_20200725T033512_022622_02AEEE_A0BA

In the Read(2) tab, as Source Product, select the product with Name: S1A_IW_SLC__1SDV_20200731T033531_20200731T033558_033693_03E7B2_D1A9

TOPS Split

Every Interferometric Wide swath (IW) consists of 3 sub-swaths and each one of maximum 9 bursts. In SNAP we can process only one swath at a time until the Deburst step. Our area of interest is located in the IW3 swath and is covered sufficiently by processing 2 bursts. We will use the TOPSAR-Split operator; this way we will reduce the total processing time.

To add the TOPSAR-Split operator for the first product, right-click right of the Read operator and go to Add → Radar →Sentinel-1 TOPS →TOPSAR-Split. Connect the Read operator to it by dragging the red arrow from the right side of Read operator towards the TOPSAR-Split operator.

Repeat the same to add the TOPSAR-Split(2) operator and connect it with the Read(2) operator.

In the TOPSAR-Split tab, Zoom in to the product and select: Subswath: IW3, Polarisations: VV, Bursts: 5 and 6 (drag the two sliders accordingly).

In the TOPSAR-Split(2) tab Zoom in to the product and select: Subswath: IW3, Polarisations: VV, Bursts: 3 to 4 (drag the two sliders accordingly).

Apply Orbit File

Next, we will apply the updated orbits to the products (See NOTE 4). Right-click and go to Add → Radar → Apply-Orbit-File. Connect the TOPSAR-Split operator to it.

NOTE 4: The orbit state vectors provided in the metadata of a SAR product are generally not accurate and can be refined with the precise orbit files which are available days-to-weeks after the generation of the product. The orbit file provides accurate satellite position and velocity information. Based on this information, the orbit state vectors in the abstract metadata of the product are updated. (SNAP Help)

Repeat the same to add the Apply-Orbit-File(2) operator and connect it with the TOPSAR-Split(2) operator.

In both the Apply-Orbit-File tab and the Apply-Orbit-File(2) tab, we will keep the default settings, and make sure that you will select the “Do not fail if new orbit file is not found” option.

Back Geocoding

Now we will coregister the two products. Image coregistration is the alignment of master and slave images, the pixels of the slave images correspond to those of the master and represent an identical area. To add the Back-Geocoding operator right-click and go to Radar → Coregistration → S-1 TOPS Coregistration → Back-Geocoding. Connect both Apply-Orbit-File and Apply-Orbit-File(2) operators to it.

In the Back-Geocoding tab set: Digital Elevation Model: SRTM 1Sec HGT (Auto Download), select the “Output Deramp and Demod Phase” option as well and leave the rest parameters as by default.

Enhanced Spectral Diversity

This operator follows the Back-Geocoding operator, it first estimates a constant range offset for each burst using a small block of data in the center of the burst and then it estimates a constant azimuth offset. Finally, the estimates from all bursts are averaged to get the final constant range and azimuth offset for the whole image.
To add the Enhanced-Spectral-Diversity operator right-click and go to Radar → Coregistration → S-1 TOPS Coregistration → Enhanced-Spectral-Diversity. Connect the Back-Geocoding operator to it.

In the Enhanced-Spectral-Diversity tab keep all the default parameters.

Write – create the output

To add the Write operator right-click and go to Add → Input-Output → Write. Connect the Enhanced-Spectral-Diversity operator to it.

In the Write tab set the following:
Name: S1_IW_SLC_20200725_20200731_Orb_Stack
Directory: …/HAZA08_LebanonDamageAssessment/Processing (Click on the icon to set the appropriate path).

Click on the icon Save to save the graph for future use. Go to the: …/HAZA08_LebanonDamageAssessment/AuxData folder and save it with the name: Graph_S1_part1.

Now go all the way up to the purple parts of the steps “Read the input products” to “Write – create the output”, fill the tabs with the appropriate parameters and finally click Run.
Once the processing is completed, the coregistered product will appear at the Product Explorer Window.

Repeat the steps “Read the inputs products” to “Write – create the ouput”

We will now keep opened the same graph we have built, and we will go again tab by tab to set the new parameters in each one, so that we create the coregistered product of the second pair of images.
In the Read tab, as Source Product, select the product with Name: S1A_IW_SLC__1SDV_20200731T033531_20200731T033558_033693_03E7B2_D1A9

In the Read(2) tab, as Source Product, select the product with Name: S1B_IW_SLC__1SDV_20200806T033445_20200806T033512_022797_02B43D_3A0B

In the TOPSAR-Split tab Zoom in to the product and select: Subswath: IW3, Polarisations: VV, Bursts: 3 to 4 (drag the two sliders accordingly).

In the TOPSAR-Split(2) tab Zoom in to the product and select: Subswath: IW3, Polarisations: VV, Bursts: 4 to 5 (drag the two sliders accordingly).

In both the Apply-Orbit-File tab and the Apply-Orbit-File(2) tab, we will keep the default settings, and make sure that you will select the “Do not fail if new orbit file is not found” option.

In the Back-Geocoding tab set again: Digital Elevation Model: SRTM 1Sec HGT (Auto Download), select the “Output Deramp and Demod Phase” option as well and leave the rest parameters as by default.

In the Enhanced-Spectral-Diversity tab keep all the default parameters.

In the Write tab set the following:
Name: S1_IW_SLC_20200731_20200806_Orb_Stack
Directory: …/HAZA08_LebanonDamageAssessment/Processing

Click Run. Once the processing is completed, the coregistered product will appear at the Product Explorer Window as well.
If we now expand the Bands folder of both the coregistered products, we will see that they contain the information of both the images used.

Sentinel-1 Processing – Part 2

In this part of the Sentinel-1 processing, we will create another graph where we will apply the necessary steps on the two coregistered products, in order to have them ready to extract the information showing the damage caused at the area due to the explosion. First, we will perform it for the coregistered product with the images before the event, and then for the coregistered product with the one image before and the one after the event.

Graph Builder

Let’s go again to Tools → GraphBuilder, right-click on the Write operator and Delete it. This time we can set the parameters in each tab, along with the creation of the corresponding operator.

In the Read tab, as Source Product, select the product with Name: S1_IW_SLC_20200725_20200731_Orb_Stack

Coherence Estimation

Since we are having in the coregistered product the Intensity bands of both products, we will use the Coherence operator, to estimate how well “correlated” the pixels between the master and slave images used, are. Right-click and go to Add → Radar → Interferometric → Products → Coherence. Connect the Read operator to it.

In the Coherence tab, set as Coherence Range Window Size: 20 and the Coherence Azimuth Window Size will automatically turn to 5.

TOPS Deburst

Now we will remove the “black space” between the two bursts (See NOTE 5). To add the TOPSAR-Deburst operator right-click and go to Add → Radar → Sentinel-1 TOPS → TOPSAR-Deburst. Connect the Coherence operator to it.

NOTE 5: There is overlapping information in every burst with its neighbouring ones, both in range and azimuth direction in order to provide contiguous coverage of the ground. Until now each burst has been processed as a separate SLC image We will merge the bursts (in azimuth direction) and preserve the phase information as well. For the overlapping region in range, merging is done between subswaths.

In the TOPSAR-Deburst tab keep the default settings (Polarizations: VV).

Multilooking

By applying this operator, we will reduce the inherent speckle noise that originally appears to the SAR images and we will obtain square pixels. To add the Multilook operator right-click and go to Add → Radar → SAR Utilities → Multilook. Connect the TOPSAR-Deburst operator to it.

In the Multilook tab keep the “GR Square Pixel” option selected and set Number of Range Looks: 8. The Number of Azimuth Looks will change to 2 and the Mean GR Square Pixel to 27.348389.

Geocoding – Terrain Correction

We will apply the Terrain-Correction operator, to convert the RADAR coordinates into geographic. To add the Terrain-Correction operator, right-click and go to Radar → Geometric → Terrain Correction → Terrain-Correction.

In the Processing Parameters tab keep the default parameters.

Subset

We will now create a Subset of the product, since we do not need the whole extent of it. You can set a different area according to your needs, or even skip the Subset step and keep the whole extent. To add the Subset operator, go to Add → Raster → Geometric → Subset.

At the Subset tab, select the Geographic Coordinates option. Copy the WKT (well know text) shown below and paste it to the text window below the map.

POLYGON ((35.46431271053535 33.92680438517669, 35.57418094686884 33.92680438517669, 35.57418094686884 33.85970773162247, 35.46431271053535 33.85970773162247, 35.46431271053535 33.92680438517669))

Click Update and then click the Zoom-in icon to see your subset on the map.

Write – create the output

To create the output of this processing, we will add the Write operator. Right-click and go to Add → Input-Output → Write. Connect the Subset operator to it.

In the Write tab set the following:
Name: S1_IW_SLC_20200725_20200731_Orb_Stack_Coh_Deb_ML_TC_Subset
Directory: …/HAZA08_LebanonDamageAssessment/Processing

Click on the Save icon to save the graph for future use. Go to the: …/HAZA08_LebanonDamageAssessment/AuxData folder and save it with the name: Graph_S1_part2.
Click Run. Once the processing is completed, the product will appear at the Product Explorer Window.

Repeat the steps from “Coherence Estimation” to “Write – create the output”

We will now keep opened the same graph we have just built, and we will go again tab by tab, to set the new parameters in each one, so that we will extract the information showing the damage caused at the area due to the explosion, of the second coregistered product (with the one image before and the one after the event). Let’s go all the way up to the steps “Coherence Estimation” to “Write – create the output”.
In the Read tab, as Source Product, select the product with Name: S1_IW_SLC_20200731_20200806_Orb_Stack

In the Coherence, TOPSAR-Deburst, Multilook and Terrain-Correction tabs, keep the same parameters as before. Check each one to make sure they have not been modified.

In the Subset tab, keep the same parameters as before. Copy again the WKT (well know text) from below and paste it to the text window below the map.

POLYGON ((35.46431271053535 33.92680438517669, 35.57418094686884 33.92680438517669, 35.57418094686884 33.85970773162247, 35.46431271053535 33.85970773162247, 35.46431271053535 33.92680438517669))

Click Update and then click the Zoom-in icon to see your subset on the map.

In the Write tab set the following:
Name: S1_IW_SLC_20200731_20200806_Orb_Stack_Coh_Deb_ML_TC_Subset
Directory: …/HAZA08_LebanonDamageAssessment/Processing

Click Run. Once the processing is completed, the product will appear at the Product Explorer Window.

Sentinel-1 Processing – Part 3

Now that we have created both final products that include the bands with the estimated coherence, we need to put them into one, in order to be allowed to apply the Change Detection operator. To do so, we will use the Create Stack operator.

Create Stack

In the 1-ProductSet-Reader tab, click on Add Opened, and all the products that are opened in the Product Explorer Window will be loaded. We need to keep only the last two. Select all the previous ones, and click on the icon to remove them.
Once there are only the 2 products that we want to stack, click on the Refresh icon. The relative information will appear at the rest fields as shown below.

In the 2-CreateStack tab, set as Initial Offset Method: Product Geolocation, and keep the rest parameters as by default.

In the Write tab set the following:
Name: S1_IW_SLC_20200725_20200731_Orb_Stack_Coh_Deb_ML_TC_Subset_Stack
Directory: …/HAZA08_LebanonDamageAssessment/Processing
Click Run. The Stack product that will contain the coherence bands of both final products, will appear at the Product Explorer Window.

Change Detection

This is the final step, where we will detect the changes between the two coherence images, and the result of it will show us the damage caused due to the explosion.
Go to Radar → SAR Applications → Change Detection.
In the I/O Parameters tab set the following:
Source: S1_IW_SLC_20200725_20200731_Orb_Stack_Coh_Deb_ML_TC_Subset_Stack
Name: S1_IW_SLC_20200725_20200731_Orb_Stack_Coh_Deb_ML_TC_Subset_Stack_change
Directory: …/HAZA08_LebanonDamageAssessment/Processing

In the Processing Parameters tab keep the default parameters and click Run.

The product with the change detection result will appear at the Product Explorer Window. Expand the Bands folder and double click on the ratio band to open it at the View Window. You can also expand the Masks folder and open the ratio_change band.
We will work with the ratio band, but you can apply the following steps to the ratio_change band if you wish as well.

Download AuxData_HAZA08.zip from the link at the start of this tutorial and unzip it to a new folder …/HAZA08_LebanonDamageAssessment/AuxData then go to the Colour Manipulation tab on the left and click on the Import Colour Palette icon. Navigate to the folder, select Damage_S1_change.cpd and click Open. At the window message that appears below, click No.

Repeat the last step. Open the palettes Damage_S1_change_1.cpd and Damage_S1_change_2.cpd.

Export Sentinel-1 product

Select the S1_IW_SLC_20200725_20200731_Orb_Stack_Coh_Deb_ML_TC_Subset_Stack_change product and go to File → Export → GeoTIFF.
The File Name will be set to S1_IW_SLC_20200725_2020_0731_Orb_Stack_Coh_Deb_ML_TC_Subset_Stack_change.tif.

At the Save in field, set the …/HAZA08_LebanonDamageAssessment/Processing path. Click Export Product.

Sentinel-2 Processing – Part 1

Close SNAP and reopen it, to load the Sentinel-2 products this time. First load the product before the event, 2A_MSIL2A_20200724T081611_N0214_R121_T36SYC_20200724T110803 and then the one after the event, 2B_MSIL2A_20200818T081609_N0214_R121_T36SYC_20200818T110855.

Create RGB image

Right click on the first product and select Open RGB Image Window. At the Profile, select: Sentinel 2 MSI Natural Colours. The B4, B3 and B2 bands will be selected. Click OK. Right click at the second product as well and open another RGB image.

Go to Window → Tile Horizontally to see the two RGB images created side by side. If we zoom in over the area of the port, we can clearly see the damages on the right one.

We can close the second image, and right click again on the first one to create another RGB image. This time, at the Profile select: Sentinel-2 MSI False-Colour Infrared. The B8, B4 and B3 bands will be selected. Click OK.
We see the different result B8 (NIR) is giving us and this is the band we will use to detect the changes over the area once we complete all the pre-processing steps.

Graph Builder

We will create again a Graph, where we will insert all the necessary operators, in order to prepare our data for the damage assessment. First, we will apply it for the product before the event and then we will repeat it for the product after the event. Go to Tools → GraphBuilder. Right-click on the Write operator and Delete it.

Read inputs

We will set the parameters in each tab along with the operators.

In the Read tab, as Source Product, select the product with Name: S2A_MSIL2A_20200724T081611_N0214_R121_T36SYC_20200724T110803

Resample

We then need to resample all the Bands so that they have a common resolution, otherwise we cannot proceed with the next steps. To add the Resample operator, right-click and go to Add → Raster → Geometric → Resample. Connect the Read operator to it.

In the Resample tab, select under the “Define size of resampled product” select the option: By reference band from source product and then select B2. This way we will resample all the bands in 20m. If you want to resample them in e.g. 10m, you can select the B1 band.

Subset

We will now create a Subset of the first product, as we did for the Sentinel-1 product before. To add the Subset operator, go to Add → Raster → Geometric → Subset.

At the Subset tab, select the Geographic Coordinates option and copy the WKT (well know text) from below and paste it to the text window below the map.

POLYGON ((35.46431271053535 33.92680438517669, 35.57418094686884 33.92680438517669, 35.57418094686884 33.85970773162247, 35.46431271053535 33.85970773162247, 35.46431271053535 33.92680438517669))

Click Update and then click the Zoom-in icon to see your subset on the map.

Write – create the outputs

To create the output of this processing, we will add the Write operator. Right-click and go to Add → Input-Output → Write. Connect the Subset operator to it.

In the Write tab set the following:
Name: S2_MSIL2A_20200724_resampled_subset
Directory: …/HAZA08_LebanonDamageAssessment/Processing

Click on the icon Save to save the graph for future use. Go to the: …/HAZA08_LebanonDamageAssessment/AuxData folder and save it with the name: Graph_S2.
Click Run. Once the processing is completed, the product will appear at the Product Explorer Window.

Repeat steps from “Read inputs” to “Write – create the outputs”

Now let’s repeat the same steps for the second product. In the Read tab, as Source Product, select the product with Name:
S2B_MSIL2A_20200818T081609_N0214_R121_T36SYC_20200818T110855

In the Resample and Subset tabs, keep the same parameters as before. Check each one to make sure they have not been modified.

In the Write tab set the following:
Name: S2_MSIL2A_20200818_resampled_subset
Directory: …/HAZA08_LebanonDamageAssessment/Processing

Click Run. Once the processing is completed, the product will appear at the Product Explorer Window.
You can open both Subset products at the View Window, to see their extent.

Sentinel-2 Processing – Part 2

Once we have both products subset, we need to put them in one, as we did before with the Sentinel-1 products. For sentinel-2 products, we will use the Collocation operator.

Collocation

Go to Raster → Geometric Operations → Collocation.
As “Master” product, select the S2_MSIL2A_20200724_resampled_subset
As “Slave” Product, click on the Add product(s) and select the S2_MSIL2A_20200724_resampled_subset and click OK.

As “Name”, set S2_20200724_20200818_collocate
Under the “Renaming of Source Product Components”, at “Rename master components” set: ${ORIGINAL_NAME}_20200724 and at “Rename slave components” set: ${ORIGINAL_NAME}_20200818.
Click Run and once the processing is completed, the collocated product will appear in the Product Explorer Window.

Sea Mask

We want to remove the pixels that correspond to the Sea. Go to Raster → Masks → Land/Sea Mask.
In the I/O Parameters tab set as “Source Product”: S2_20200724_20200818_collocate, as “Name”: S2_20200724_20200818_collocate_msk and as directory …/HAZA08_LebanonDamageAssessment/Processing

In the Processing Parameters tab, at the “Source Bands” keep the Ctrl pressed and select only the B8_20200724 and B8_20200818 bands. Check the Mask out the Sea option and click Run.

If we open the B8 bands of the masked, collocated product, we will see that the sea has been removed.

Band Maths

Next, we want to create the final band that will show us the areas that have been damaged. For this, we right-click on the S2_20200724_20200818_collocate_msk product and click on the “Band Maths”. Set as “Name”: S2_change and deselect the “Virtual (save expression only, don’t store data) option”. Then click on ”Edit Expression” and insert the following:

B8_20200724 - B8_20200818

Click OK in both windows. The new band will be created in the product. Double click and open it.

If you want to save this band, right click on the product and click on “Save Product”.

Click on the Import Colour Palette icon. Then, navigate to the /Training/HAZA08_LebanonDamageAssessment/AuxData folder, select the Damage_S2.cpd and click Open.
At the window message that appears below, click No.

In red we can see the areas that have been damaged from the explosion.

Reprojection

The last step is to reproject the product so that we can insert it in QGIS. Go to Raster → Geometric Operations → Reprojection.
In the I/O Parameters tab set the following:
Source: S2_20200724_20200818_collocate_msk
Name: S2_20200724_20200818_collocate_msk_reprojected
Directory: …/HAZA08_LebanonDamageAssessment/Processing

In the Reprojection Parameters tab, keep all the parameters as by default. Click Run.

Export Sentinel-2 product

Select the S2_20200724_20200818_collocate_msk_reprojected product and go to File → Export → GeoTIFF. Click on Subset and select only the S2_change band. Click OK. The File Name will be set to S2_20200724_20200818_collocate_msk_reprojected.tif. At the Save in field, set the …/HAZA08_LebanonDamageAssessment/Processing path. Click Export Product.

QGIS Visualization

Open QGIS and load the S1_IW_SLC_20200725_20200731_Orb_Stack_Coh_Deb_ML_TC_Subset_Stack_change.tif and the S2_20200724_20200818_collocate_msk_reprojected.tif at the Layers panel. Right-click on the S1 product and go to “Properties”.
In the “Symbology” tab, set at “Render Type” Singleband Pseudocolour
On “Style”, click and select “Load Style”. Navigate to the …/HAZA08_LebanonDamageAssessment/AuxData folder, select Damage_S1.qml and click Open to load the palette.
We see that there is a black background in both images, so we will remove it by going in the “Transparency” tab and set as “Additional no data value”: 0

Click OK.

Repeat the same steps for the Sentinel-2 product in the Symbology and the Transparency tabs, and when going at Load Style, select the Damage_S2.qml palette.

Go to Web → OpenLayers plugin → Bing Maps → Bing Aerial with labels (See NOTE 6).

NOTE 6: In case the OpenLayers plugin is not installed, click on Plugins → Manage and Install Plugins. Select the “All” tab on the left side panel and write “OpenLayers plugin” on the search box. If you cannot find it, go to “Settings” and select the “Show also experimental plugins” option. Go back to the “All” tab, select the plugin on the list and click “Install Plugin”. Restart QGIS to finalize the installation.

Now that we have added a basemap at the background, we can select which tif. files we want to visualize on the top of them.

For example, we see the Sentinel-1 product with the Damage_S1_2 palette:

the Sentinel-2 product with the Damage_S2 palette:

and combined damage detection from both S-1 and S-2 products:

Further reading and resources

https://sentinel.esa.int/web/sentinel/missions/sentinel-1 – Sentinel-1 Mission
https://sentinel.esa.int/web/sentinel/missions/sentinel-2 – Sentinel-2 Mission

Search
Browse by topic
Related Tutorials

OIL SPIL MAPPING WITH SENTINEL-1

Reference folder structure Folder structure referenced in the exercise is as follows: /OCEA03_OilSpill_Kuwait /Original – should contain yo…

Doing business with us

EO expertise at your service

We are used to prepare and release training materials, but also to work in partnership with the Earth Observation industry and we may complement your business case with training or consultancy services. Do not hesitate to get in touch with us to get an offer for a customized training event or EO consultancy services…