Working with Sentinel 1 - RADAR : Case Study in Jajarkot EarthQuake

The goal of this tutorial is to share my experience and steps I followed when I tried to work with Sentinel-1 Radar Data to findout the displacement caused by the earthquake in Nepal.

Table of Contents

1. Study Area
2. Preparation
3. Download Image
4. Data Preparation
5. Image Preprocessing
6. Generation of Topographic Interferogram
7. Generation of Differential Interferogram
8. Creation of Displacement Map
9. Terrain Correction
10. Export Displacement Map
11. Sources and Credit

Study Area

JajarKot Nepal : https://www.openstreetmap.org/relation/4586357#map=10/28.9246/82.3151

Source : https://edition.cnn.com/2023/11/03/asia/nepal-earthquake-northwest-hnk-intl/index.html

There was a deadly earthquake in Early November 2023 which killed more than 157 people as per CNN . I wanted to figure out how much displacement it has cased in the area , Did topography went up ? or went down. As we know due to tectonic plates movement it will always result in displacement !

Preparation

Lets first verify that event exists , get the event details from USGS . USGS mainatains the eartquake catalog which will give the event date, magnitude and approx coordinate for the event .

https://earthquake.usgs.gov/earthquakes/search/

Prepare your filter accordingly to your event , For us as CNN reported the incident in november 3 , I am keeping filter from Novemeber 2 to 4 in Jajarkot Area . This is how my search filter result looks like : https://earthquake.usgs.gov/earthquakes/map/?currentFeatureId=us7000l8p5&extent=22.2586%2C75.81665&extent=31.72817%2C91.70288&range=search&sort=largest&search=%7B%22name%22%3A%22Search+Results%22%2C%22params%22%3A%7B%22starttime%22%3A%222022-01-01+00%3A00%3A00%22%2C%22endtime%22%3A%222024-12-08+23%3A59%3A59%22%2C%22maxlatitude%22%3A32.319%2C%22minlatitude%22%3A24.121%2C%22maxlongitude%22%3A89.868%2C%22minlongitude%22%3A78.354%2C%22minmagnitude%22%3A5%2C%22orderby%22%3A%22magnitude%22%7D%7D

We found event details , it's coordinate and confirmation from USGS , Now we can move forward to downloading the Sentinel Image

Download Image

Filters

I am using the ASF Data Search by NASA : https://search.asf.alaska.edu/#/

Click on filters Nearby search

Here are the parameters for my Search filters :

• Area of Interest : POINT(82.153 28.858) (You can convert your coords to WKT with https://wktmap.com/ )
• Start Date : 2023/10/01
• End Date : 2023/11/30 (To be on safe side I am trying to find image within two months frame to filter out the best image before and after the event )
• File Type : L1 Single Look Complex (SLC) ( SLC has both amplitude and phase , in GRD we have only amplitude . Here in displacement we need phase hence SLC )
• Beam Mode : IW ( Because IW- Interforemetric covers Orbit and EW- extra wide swath covers Polar region )

Source : My professor Karima Hadj-Rabah & Zahra Dabiri

• Polarizations : VV ( Polarization doesn't really matter in phase differences , It matters in classification , Here we are selecting VV just to reduce the image size )
• Direction : Ascending ( It doesn't matter as well because at the end we will do geocoding and correct the image so image going to be reotated anyway )

Search Results :

I am selecting two image based on my region of interest :

• After the event :: Novemeber 10, 2023
• Before the event :: October 17,2023

if you look closely to the footprint image top left has my region of interest specifically

You need to create account in order to download the image.

Note :
As some of the required steps are computationally intensive, it is good to store the data at a location which offers good reading and writing speed. If your computer has an internal SSD, processing should be done there to ensure best performance. (Credit : ESA Tutorial)

Data preparation

For data preparation and image processing and analysis I am going to use a tool called SNAP . SNAP (Sentinel Application Platform) is a earth observation processing open-source software developed by ESA

Download SNAP :

https://step.esa.int/main/download/snap-download/

SNAP is available for windows , mac and linux
I will be using linux , I downloaded all tollboxes

Setting up SNAP

chmod +x esa-snap_all_linux-11.0.0.sh

./esa-snap_all_linux-11.0.0.sh

SNAP Operations

Open Snap it , Interface looks like this

Load image

Our downloaded image would be in .zip folder , Don't unzip them simply open them in SNAP

File > Open Product > Select your both before after images > Open

Split data to region of interest

Goto

Radar > Sentinel-1 TOPS > S-1 Tops Split

Navigate to Processing Parameters

• Subswath : IW1
• Polarisation : VV
• Bursts : 6 to 9 (For my area of interest bursts 6 to 9 covered the area , You should select the subswath and bursts accordingly . If you remember my are of intrest was in top left of the image , hence I selected the subswath that covers that area and bursts to minimize it )

Lets rename the output file for readability

Hit Run

You should have the splitted file and now

Visualize the intensity

Expand the image > Bands > subswath name > intensity_* Double click

Now repeat the process to split another image as well , Remember to click on the image that you are splitting or change the path

Visualize the image side by side from

Window > Tile Vertically or horizontally

Image preprocessing

Orbit Correction

Technically we suppose that satellite is following perfect line , But in reality there is vibration in satellite and is not really stable . Hence we need to apply correction on the orbit path

Radar > Apply Orbit File

Go to Processing Parameters and Check Do not fail if new orbit file is not found to make sure our process doesn't fail even though orbit file not found ( Because orbit file is only available after two weeks of the acqusition at the moment )

Repeat the process for another image as well . You should have those two new files in your window

Generation of Topographic Interferogram

A topographic interferogram highlights differences in phase caused by terrain elevation. It is a key step in extracting the displacement component. You must align the pixels between two different image ( Before and after one ) to make the comparison hence coregistration

Source ( Prof : Karima )

Cleanup

Lets Clean up the file we don't need

Control select the split file and original image and Right click > Close Products

No need to save Hit NO

Back Geocoding and Enhanced Spectral Diversity

The S-1 Back Geocoding operator coregisters the two split products based on the orbit information added in the previous step and information from a digital elevation model (DEM) which is downloaded by SNAP.

Go to

Radar > Coregistration > S-1 Tops Coregistration > S1- Back Geocoding

Click on Add Opened

Make sure before image is on the top and after image is at the bottom

You should see your image there

Go to Back-Geocoding

Select the DEM source : For me I am selecting SRTM 1Sec

Press Run

You shuold have new stack image with both of your image stacked together , You shouldbe able to see your before and after image intensity

Enhanced Spectral diversity

In order to improve the quality of the after image as related to before , Inorder to remove inospheric error this step is used .

Go to :

Radar > Coregistration > S1 TOPS Coregistration > S-1 Enhanced Spectral Diversity

Hit Run

It should add a file with filename_*stack_esd

You can change the visualization parameter from Color Palette

Interferogram Formation

This is where we actually calculate the phase difference between two images along with coherence

Go To

Radar > Interferometric > Products > Interferogram Formation

if you check on Processing Parameters :

You should see subtract flat-earth phase ticked

Reminder : Make sure input image for the new step is always from the previous step

You should have new filenamewith_ifg extension

TOPSAR Deburst

To remove the seamlines between the single bursts we use TOPSAR Deburst

Go TO :

RADAR > Sentinel 1 - TOPS > S-1 TOPS-Deburst

Click on Run

if you visualize the result , it will remove the small black lines between bursts in the image

Generation of Differential Interferogram

By generating a differential interferogram, we remove the topographic effects, extracting the phase
changes caused by ground displacement.

Source : Prof Karima

Interferogram formation

Topographic Phase Removal

GOTO :

Radar > Interferometric > Products > Topographic Phase Removal

Navigate to Processing Parameters

Choose the DEM source : For me I chose SRTM 1sec

Check Output Topographic Phase band and the elevation band ( Lets visualize the elevation as well )

Hit Run

I get the output like this along with elevation band

Multilooking

We need to create square pixels to make the resolutions same in range(approx 3m) and azimuth (approx 7m). we are used to see pixels in square

GOTO

RADAR > SAR utilities > Multilooking

Navigate to processing parameters and change number of looks for range and azimuth.

For me I am using

• Range Looks : 8
• Azimuth Looks : 2

Because it results approx 30m resolution which is enough for me because usualy displacement will happen in larger level, You can see the resolution after the value changes

This is how output looks like

Goldstein Phase Filtering

Lets filter the noise so that we can improve visual interpretability of differential interferogram.

GOTO :

Radar > Interferometric > Filtering > Goldstein Phase Filtering.

Navigate to processing parameters and play with the filter value , For me I am using 0.9 as There is quite some noise in my image

Hit Run

This is what I got as output , I have some noise in the right bottom part of image : for now it is okay but you can also play with it and increase or decrease the filter

Creation of Displacement Map

Dowload snaphu plugin

The wrapped phase values need to be unwrapped to obtain continuous displacement
measurements.

For this output : We need this plugin to be installed : https://step.esa.int/main/snap-supported-plugins/snaphu/

For Debian based distributions, You can install it directly from the apt

sudo apt get update
sudo apt install snaphu

Export env for snaphu

Lets create env needed for snaphu

GO TO
Radar > Interferometric > Unwrapping > SNAPHU Export

Navigate to SnaphuExport

Select the folder where output would be stored ( Recommended to create new folder ) , For me I created snaphu-export folder with snaphu_output filename

Change No of tile rows and tile columns and number of processer based on your laptop spec , For me I used 8 processers and 10 rows/column . Bearing in mind if you set 1 tile per time it will disable the multiprocesser option ! It is advised to process tiles in parallel

Hit Run

Now navigate to the folder you specified you should see following structure

Now open snaphu.conf in one of your editor ( Notepad ++ , VSCode , Kate any you prefer )

Now Here :

• Comment out logfile
• Find the name of coherence file from the folder and Specify it in CORRFILE

Your final conf should look like this

Now open terminal in the folder dir , Remember in the folder dir

Now you will find the command on top of the file, Copy & run the command in terminal

The time differs based on the performance of laptop & the number of tiles you are processing based on your parameter in snaphu export , For me it took 1 min 43 sec

Import from Snaphu

Now lets import the processing output from snaphu
Goto :

Radar > Interferometric > Unwrapping > Snaphu Import

Select the output file , Remember OUTFILE information is in snaphu.conf

Navigate to : 2 Read-Unwrapped-Phase & Open the result from snaphu

Select .hdr file format

Navigate to Snaphu Import and Tick on Don't save wrapped inferogram in the target Product

Navigate to Write tab

And edit the name add _unwrapped in the name for disntinction

Hit Run

This is how my output looks like

Generate of Displacement Map

This step converts the unwrapped phase into displacement values, enabling the analysis of
ground movement.

Go to

Radar > Interferometric > Products > Phase to Displacement

Hit Run

This is how my displacement map looks like

Terrain Correction

Terrain correction georeferences the displacement map to a coordinate system and removes
distortions caused by topography, ensuring accurate spatial analysis.

Goto :

Radar > Geometric > Terrain Correction > Range-Doppler Terrain Correction

Navigate to Processing Parameters and Select DEM , I selected SRTM 1Sec

Hit Run

This is what I get as output

Lets understand this displacement Map

I changed the visualization from left bottom window and changed editor to slider

If you look into the value , Blue area was the area where terrain went down by approx 16 cm and nearby area terrain raised by 23.7 cm

Export displacement Map

Right click on the image

And Click on Export View as Google Earth KMZ

Lets visualize this KMZ in QGIS

Set layer transparency to 50% so that we can visualize area where terrain went down and terrain went up after earthquake

here is the output and as you can see area is nearly earthquake area and hence we can verify this output makes sense !

You can also export image as geotiff and other formats as well

Thank you for reading !

Sources and Credit :

• https://www.youtube.com/watch?v=ZPMRaztNbVU&t=393s
• Prof. Karima Hadj-Rabah
• My classmate Sahar for her help during blog creation
• Paris - Lodron University , Z-GIS Department
• ESA Snaphu Tutorial