What happened in Rishi Ganga Valley?

On the morning of February 7th 2021, in high regions of mount Nanda Devi in the Himalayan range, Uttarakhand, India, a huge chunk of rock and ice fell, triggering a landslide and wreaking havoc in the downstream valley. This mass movement of rocks, ice, mud, and debris travelled 1,500 metres and eventually cascaded into the rivers down below in the Rishi Ganga valley. The sudden increase in volume in the rivers produced deadly flash floods that damaged and destroyed everything in their path, including critical infrastructure like powerplants.

As shown in the figure below, the deluge came from south of village called Raini and surged westward towards Tapovan dam, washing away a bridge and the 13.2 megawatt Rishiganga project on the banks a tributary of the Dhauliganga, as well as causing major damage to the Tapovan Vishnugad project

Satellite images: PLEIADES © CNES 2021, Distribution Airbus DS / produced by Earthrise and sourced from Reuters.

Analysis of seismological data from nearby observatories ruled out earthquakes as the trigger for the landslide (Pandey, et. al 2021). Remote sensing data and other geospatial tech have contributed significantly to explaining the likely cause of the landslide, long before actually the area could be physically visited. Using scalable infrastructure, on-demand imagery, and off-the-shelf analytics empowers decision-makers to act fast and could even enable risk-managers to plan ahead.

So you may be wondering, what geospatial data and analytics can be used to assess landslides like the one in Rishi Ganga valley:

• Radar pairs or stacks (SLC, pre-event)
  • Differential Interferometrics (DinSAR) techniques
• Additional interpretation with optical data or DTM
  • Visual analysis of before and after (optical)
  • Slope analysis for further risks (DTM)

What is DinSAR?

Advanced Differential SAR Interferometry or DinSAR is a satellite-based remote sensing technique used to monitor ground deformation with high resolution over short periods of time. The basis of SAR interferometry is the comparison of the phase of two or more radar images of a given scene acquired from slightly different positions or at different times. One key result of this comparison is a flattened interferogram, which is essentially an image with coloured gradients that depicts the relative displacement between the first and second acquisition, where the different colours represent the different phase shifts.

. Displacement is derived from the phase shift of two SAR acquisitions. Image from Schindler et al., 2016.

In the cases where there’s observable surface movement (land deformation) between the two radar acquisitions, the flattened interferogram contains information about both the topography and the deformation. Differential InSAR is the process of removing the topography from the flattened interferogram in order to get quantitative results about past ground displacements or movement. The benefits of this method compared to field methods is unmatched scalability and high levels of accuracy, both in detecting the precise location where movement has occurred and the degree of movement.

An interferogram created by using two SAR images that are acquired in two different times and maps the phase shift caused by a movement of the surface between the two acquisitions. Image: Geoscience Australia.

Identifying Precursors for Landslides with SAR and DinSAR

Since SAR data can cover a very large spatial area with centimetre-scale accuracy, DinSAR can be useful in monitoring large areas of landslide risk especially for locations difficult to access. One main limitation to consider can be the difficulty of seeing through vegetation to the ground beneath.

DinSAR has greatly increased the range of landslide research and monitoring applications since their development in the early 2000s. By monitoring the surface deformation of mountainous regions, the movement of rock and sediment can be tracked and that can lead to the uncovering of risks, which can help prevent disaster. DinSAR-detected deformation for landslides can be both vertical and horizontal, as it involves rock and debris moving along a slope.

DinSAR steps: how does UP42 help the landslides use case?

Landslides occur in every corner of the world. The following example shows how the UP42 platform can be used to identify the precursors for any landslide event. This example leverages both the UI for a simple 1-off approach and the UP42 SDK, with which enables quick, scalable access and analysis of SAR and other geospatial data.

NASA Landslide Viewer, global landslide catalog. Each red dot represents a reported landslide over the past 10 years.

Using SAR data within the UI

For this example, we will look at the time period just leading up to the landslide in Rishi Ganga Valley. In order to do this we need to:

1. Easily click together blocks in UI, saving the workflow
2. Find two(+) SLC datasets before the event, taken under the same conditions**
3. Specify AOI and data IDs in the job parameters
4. Run the workflow, download results, and visualize

In the UI, we create a project>workflow and add the S1 SLC SAR data. Catalyst, a PCI Geomatics brand that provides cutting-edge algorithms, has a DinSAR algorithm on the UP42 marketplace called Ground Displacement that calculates cumulative ground displacement as well as deformation velocity and is perfect for applying to this land use case. This is the second block to add and you can click save to save the workflow.

Now you can navigate to the catalog search. Note that when using this algorithm, the number of pixels in either direction must be at least 2000. To be on the safe side, we select an AOI of ~350 km2.

This would look like: Left- Selecting SLC pairs and and AOI with the UP42 Catalog Search. Right- Clicking together the Sentinel-1 SLC block and Catalyst's Ground Displacement block for the DinSAR calculation using the UP42 UI.

You might be asking, how do we ensure the SLC datasets were taken under the same conditions? There’s a couple ways this can be done. First, you can ensure all key search criteria are the same (see the SDK example in the next section) or second, you can check that the footprints of the SLC data match exactly via the catalog search as shown in the video below:

Understanding Rishi Ganga Valley’s landslides with SAR data via the API

Alternatively, we could leverage the API to scale and leverage the UP42 Python SDK to gain a deeper understanding and directly execute some visualizations. For a quick and easy way to replicate this use case, please see the prepared jupyter notebook by Catalyst or for further information using the API, consult the documentation.

We start by installing the UP42 Python SDK & importing it: Create a virtual environment:

mkvirtualenv --python=$(which python3) up42-py

Activate the environment:

workon up42-py

Install the module:

pip install up42-py

Install Jupyter Lab:

pip install jupyterlab

Import the UP42 Python library:

import up42

Feel free to import any other libraries you’re interested in using, especially for your preferred visualization. To get started with UP42 SDK, we will need to create a project via the browser console and get the project authentication keys. These keys need to be stored as config.json in the current working directory. A config file should be structured as follows:

{
    "project_id": "abc",
    "project_api_key": "xyz"
}

To know how to obtain the project API key and project ID, please refer to this page in our documentation.

Now that we have all the prerequisites covered, we can initialize the project, run our first job and assess the Sentinel-1 AOI clipped images over the area of interest.

from pathlib import Path

# Reading the credentials from proj_default.conf, authenticating, and initializing the project
project_config_file = Path.home().joinpath(".up42/proj_default.conf") 
up42.authenticate(cfg_file = project_config_file)
aoi =  up42.read_vector_file("data/rishi_ganga_valley.geojson")
my_catalog = up42.initialize_catalog()

Alternatively, you might initialize your project directly with the project ID and API key, like below.

up42.authenticate(project_id="1234", project_api_key="ABCD")
aoi =  up42.read_vector_file("data/rishi_ganga_valley.geojson")
catalog = up42.initialize_catalog()

In order to run the job, we need to provide a few search parameters. Note the catalog follows STAC specs. For this example we will use Sentinel-1 SLC AOI clipped data block. A full list of data blocks available through UP42 can be found in the UP42 Marketplace.

# Construct search or input parameters
search_params = my_catalog.construct_parameters(geometry = aoi, 
                                   start_date = "2021-01-15",
                                   end_date = "2021-02-05",
                                   sortby = "acquisitionDate",
                                   ascending = False,
                                   sensors = [ "sentinel1" ],a
                                   limit = 100)
#Check data blocks to be considered (all S1 data as specified by sensors = [“sentinel1”]
data_block_names = search_params["query"]["dataBlock"]["in"]

#Select only the SLC S1 data type, excluding the GRD types
from re import compile, match
re_slc = compile(".*-slc-.*")
search_params["query"]["dataBlock"]["in"] = [b for b in search_params["query"]["dataBlock"]["in"] if re_slc.match(b)]

#Run the search
search_results = my_catalog.search(search_parameters = search_params);

When selecting the images for DinSAR, it is important to ensure that all images, within the images series being assessed, are taken under the same conditions. For interferometric algorithms, there are three components of the search criteria that we can easily delve further into with the SDK, in order to guarantee the algorithm works properly. Namely, we need to look at the sceneIds, the relative orbit number, and the look direction.

# Compose the result list from the sceneIds, relative orbit number and direction.
rl = [(sid, acq_date, rp["orbit"]["relativeNumber"], rp["orbit"]["direction"]) for sid, acq_date, rp in zip(search_results.sceneId, search_results.acquisitionDate, search_results.providerProperties)]
# Import Pandas to convert above into a data frame
from pandas import DataFrame as df
# Compose scene list
scene_list = df(rl, columns = ["sceneId", "acquisitionDate", "orbitRelativeNumber", "direction"])
# Assess scenes to find compatible pairs or series
scene_list[scene_list.direction == "ASCENDING"]
scene_list[scene_list.direction == "DESCENDING"]

The output should look like the following: Table of Sentinel-1 SLC Scene IDs, acquisition dates, orbit number, and direction for DinSAR analysis.

From the tables generated above, we can see that it is the images from February 3rd and January 22nd fit our criteria, recalling:

• Same relative orbit number.
• Same direction.
• One of the acquisition dates is just before the landslide (Feb 7th 2021).

SAR and DinSAR workflow with the API

From back on the console, we can easily access the workflow ID of the DinSAR workflow we've already created. With this ID, we'll initialize the workflow, set the parameters, and ultimately run the job.

insar_workflow_id = "1234"
insar_workflow = up42.initialize_workflow(workflow_id = insar_workflow_id)

Once the DinSAR workflow is initialized, we need to pass through the AOI again into the workflow and specify the job parameters, such as the Scene IDs we chose from before.

insar_job_params = insar_workflow.construct_parameters(geometry = aoi["features"][0]["geometry"], 
                                                       geometry_operation = "intersects", 
                                                       scene_ids = final_scene_id_list)
insar_job_params["catalystpro-insstack:1"].update(aoi_geojson = aoi["features"][0]["geometry"])

Now we're ready to run the job. But before we run the job, we can get a cost estimate of the job run.

#Estimate the cost
insar_workflow.estimate_job(input_parameters = insar_job_params)
#Run the job
insar_job = insar_workflow.test_job(input_parameters = insar_job_params, 
                                name="Rishi Ganga Valley 2021 Landslide (SDK)", 
                                track_status = True)
#Download the results
my_path = Path.home().joinpath("up42/catalyst-demo/results-rishi-ganga")
insar_job_results  = insar_job.download_results(output_directory = my_path)

Finally, we can plot some of the DinSAR results, specifically a corrected cumulative displacement geotiff, directly with the matplotlib and rasterio libraries.

from matplotlib import pyplot
import rasterio as rio
displacement_stack_raster = insar_job_results[-3]
with rio.open(displacement_stack_raster, 'r+') as ds:
    ds.nodata = -32768
displacement_stack_raster = rio.open(displacement_stack_raster)
# Set axes and plot size
fig, ax = mpl.pyplot.subplots(1, figsize = (18,18))
# Visualize with a color map
color_map = mpl.cm.get_cmap('RdYlGn')
color_map = color_map.reversed()
pyplot.colorbar(pyplot.imshow(displacement_stack_raster.read(1), 
                              cmap = color_map, vmin = -20, vmax =100),
                              shrink = 0.85,
                              label = "Cumulative Displacement (mm)")
pyplot.title("Cumulative Displacement (mm) EPSG:32644")
rio.plot.show(displacement_stack_raster, cmap = color_map)

If we zoom into these results, we can identify that there was severe movement where the landslide occurred the week prior to when the landslide occurred.

Mapping landslide susceptibility can be critical for risk mitigation and regional planning. DinSAR is a great way to get situational awareness of land movement from a distance. There are even more components that can be derived from geospatial data that can help indicate landslide risk like slope and aspect from DEMs, NDVI from optical imagery, and even weather data. All of these contribute to a comprehensive picture of landslide risk.

For more information, check out CATALYST.Earth's capabilities on the UP42 marketplace, as well as the many previously mentioned data sources like Meteomatics weather or SPOT.