Learnings from rapid response efforts to remotely detect landslides triggered by the August 2021 Nippes earthquake and Tropical Storm Grace in Haiti

Abstract

On August 14, 2021, a Mw 7.2 earthquake struck the Tiburon Peninsula of western Haiti triggering thousands of landslides. Three days after the earthquake on August 17, 2021, Tropical Storm Grace crossed shallow waters offshore of southern Haiti triggering more landslides worsening the situation. In the aftermath of these events, several organizations with disaster response capabilities or programs activated to provide information on the location of landslides to first responders on the ground. Utilizing remote sensing to support rapid response, one organization manually mapped initiation point of landslides and three automatically detected landslides. The 2021 Haiti event also provided a unique opportunity to test different automated landslide detection methods that utilized both SAR and optical data in a rapid response scenario where rapid situational awareness was critical. As the methods used are highly replicable, the main goal of this study is to summarize the landslide rapid response products released by the organizations, detection methods, quantify accuracy and provide guidelines on how some of the shortcomings encountered in this effort might be addressed in the future. To support this validation, a manually mapped polygon-based landslide inventory covering the entire affected area was created and is also released through this effort.

1 Introduction

Every year, geohazards cause human fatalities and adverse impacts to communities, and are expected to increase in frequency and severity with anticipated climate change scenarios in the twenty-first century (Gariano and Guzzetti 2016; Jakob 2022). Economically developing nations bear the brunt of these impacts, where depressed socioeconomic status is a driving indicator for increased fatalities due to landslides (Kirschbaum et al. 2010, 2015; Dowling and Santi 2014; Froude and Petley 2018; Pollock and Wartman 2020), flooding (Flanagan et al. 2011; De Silva and Kawasaki 2018) and earthquakes (DesRoches et al. 2011).

Local geohazard expertise is often lacking during response operations, especially in developing nations (Collins et al. 2020). The lack of in-country geohazard expertise, or the size of the disaster, may trigger the deployment of international geohazard specialists to the site. Once onsite, the role of the specialists is to provide contextual information and awareness regarding the disaster to first responders performing response operations (e.g., Huang et al. 2015).

In addition to onsite assistance, remote analysis of the geohazard event can also aid in the response effort, especially for geohazards that span large spatial extents. Landslides are a pervasive geohazard causing hundreds to thousands of fatalities and billions of dollars in damage each year (Kirschbaum et al. 2015; Froude and Petley 2018). Recent technological advances in remote sensing, satellite technology, and computational methods have permitted faster and more detailed analysis of landslide events from the global research community. The pace of this type of rapid response is increasing, with many researchers and research organizations providing information in days to weeks following the landslide event (Yun et al. 2015; Kirschbaum and Stanley 2018; Emberson et al. 2021; Stanley et al. 2021; Allstadt et al. 2022).

Many large landslide events are comprised of hundreds to thousands of landslides, which require remote analysis for rapid assessment. This rapid assessment can provide key information regarding the ongoing hazard potential, triage and resource allocation decisions, or estimating potential impact extents. Rapid mapping and monitoring of landslides is commonly conducted to provide situational awareness to first responders (e.g., Manconi et al. 2014; Iverson et al. 2015; Casagli et al. 2016; Gallen et al. 2017; Stock et al. 2018; Williams et al. 2018; Collins et al. 2020; Wooten et al. 2022). Rapid mapping can be done using manual and automated mapping methods (Guzzetti et al. 2012; Scaioni et al. 2014). Manual mapping is a traditional technique where a human digitizes landslides via visual interpretation of aerial or optical satellite imagery. As a human uses their expert judgment for digitizing landslides, this method tends to produce less false positives and can simultaneously identify landslide impacts during mapping. However, this method is labor- and time-intensive and quality is variable and dependent on the experience of the human analyst (Guzzetti et al. 2012; Hölbling et al. 2017). This method becomes cumbersome when mapping is to be done over a large area and timeliness of the results are important (Galli et al. 2008). Automated or semi-automated mapping involves the use of supervised or unsupervised classification techniques to detect landslides from aerial or satellite imagery. Pixel-based (Nichol and Wong 2005; Borghuis et al. 2007; Parker et al. 2011; Burrows et al. 2019, 2020; Jung and Yun 2020), object-based image analyses (OBIA) (Martha et al. 2010, 2016; Lu et al. 2011; Stumpf and Kerle 2011; Hölbling et al. 2012, 2015; Amatya et al. 2019, 2021a, b; Adriano et al. 2020; Esposito et al. 2020; Comert 2021) and deep-learning based (Ghorbanzadeh et al. 2019; Sameen and Pradhan 2019; Meena et al. 2021; Nava et al. 2022; Bhuyan et al. 2023) methods have been extensively used for landslide mapping all around the world. These methods are scalable and can be used to map large areas with less human interaction. However, these methods can result in more false positives compared to a human analyst and require further interpretation (Li et al. 2014).

In this study we revisit multiple remote sensing-based landslide rapid response inventories/maps (Table 1) that were used to identify thousands of landslides triggered by the August 2021 Nippes earthquake and Tropical Storm Grace in Haiti. We aim to quantify the accuracy of automated and semi-automated rapid response landslide inventories/maps produced by various government, academic, and private organizations, and draw useful insights which might be beneficial for rapid response in the future. We also release a comprehensive, manually mapped, landslide inventory in polygon format generated using multiple optical sensors. This inventory enabled validation of rapidly generated landslide inventories/maps and may support future landslide studies in this region.

Table 1 Summary of rapid response landslide products generated by various organizations

1.1 The 2021 Nippes earthquake and Tropical Storm Grace

On August 14, 2021, a Mw 7.2 earthquake struck the Tiburon Peninsula of western Haiti (see U.S. Geological Survey, Earthquake Hazards Program, https://earthquake.usgs.gov/earthquakes/eventpage/us6000f65h/executive). The earthquake occurred along the Enriquillo-Plantain Garden fault (EPGF) zone, west of and on the same fault zone of the 2010 Mw 7.0 earthquake that devastated much of Port-Au-Prince (DesRoches et al. 2011). The rupture started in the center of the Tiburon Peninsula, about 125 km west of the capital, Port-au-Prince, and ruptured primarily westward with reverse and left-lateral strike-slip motion along the EPGF zone. The fault rupture terminated in the vicinity of Pic Macaya National Park, a zone of rugged topography in the interior of the Tiburon Peninsula (Fig. 1). Ground accelerations exceeded 0.5 g in the vicinity of Pic Macaya National Park. The earthquake triggered thousands of landslides in a region stretching approximately 100 km west of the epicenter (Martinez et al. 2021a; Zhao et al. 2022; Havenith et al. 2022). The majority of landslides were concentrated in or near the Pic Macaya National Park.

Fig. 1

Conditions associated with the Nippes earthquake and Tropical Storm Grace sequence, August 2021. Peak Ground Acceleration (PGA, unit of %g) contours illustrate decay of ground shaking away from the epicenter (source: U.S. Geological Survey ShakeMap, Version 15: 2022-08-11). Tropical Storm Grace traveled from east to west in shallow waters to the south of the island (storm path shown in blue). Rainfall estimates are from the NASA Global Precipitation Mission (GPM) and are reported as the cumulative value of rainfall during the storm from August 16–17, 2021

On August 17, 2021, three days after the Nippes earthquake, Tropical Storm Grace crossed shallow waters offshore of southern Haiti with tropical storm intensity winds (Reinhart et al. 2022). Haiti does not have a dense network of in situ rain gauges, however, gauges in the eastern Tiburon Peninsula town of Savane Zombi and in the central peninsula town of Fonds des Negres reported 184 and 329 mm of rainfall during Tropical Storm Grace, respectively (Reinhart et al. 2022). These data are corroborated by satellite-derived rainfall estimates from GPM (Huffman et al. 2019; Fig. 1) which suggests > 250 mm of rain fell across the majority of the Tiburon Peninsula. Rainfall was significant enough to trigger landslides, flash flooding, and river flooding across much of the central and southern Tiburon Peninsula (Reinhart et al. 2022; Whitworth et al. 2022).

Due to the close timing of occurrence of the Nippes earthquake and the passage of Tropical Storm Grace, it remains difficult to ascertain the level of landslide impact from each hazard. The cumulative impacts of the earthquake affected at least 800,000 people and resulted in at least 2,248 fatalities, 12,763 injuries, and 137,500 damaged and destroyed homes (OCHA 2021).

1.2 Disaster response programs and responding organizations

Two government agencies, three universities and one private organization responded on behalf of three separate response programs. The responding entities and response program are introduced in the following sections.

1.2.1 U.S. geological survey (USGS)

The USGS works to monitor, assess, and conduct targeted research on a wide range of natural hazards so that policymakers and the public have the understanding they need to enhance preparedness, response and resilience. For domestic geologic disasters, such as an earthquake, volcanic unrest, or significant landslide event, the USGS typically responds with scientists that provide local officials information for decision making. For international earthquake events, the USGS routinely responds to major earthquakes by providing and updating earthquake information products and models. The response may include remote review of satellite imagery and provide information on what we see when relevant for most major events that rupture the surface, damage buildings, or cause ground failure. However, the USGS does not typically assist foreign countries, such as systematically mapping landslides, unless formal assistance is requested. For this event, USGS Earthquake Hazards Program and Landslide Hazards Program scientists collaborated with the USAID-USGS Landslide Disaster Assistance Team (LDAT) to rapidly map landslides from optical satellite imagery. LDAT typically leads landslide technical capacity building with a country's respective geologic survey or similar government institution; however, LDAT also provides technical assistance for some international landslides and related geologic disasters, providing remote sensing and situational awareness to responding humanitarian organizations. For the landslide event in Haiti, LDAT and USGS manually mapped and compiled initiation point (analyst interpretation of highest point on the landslide scar)-based landslide inventory using satellite optical imagery from Copernicus Sentinel-2 (10 m resolution), Maxar WorldView (0.3–0.5 m resolution), and Planet Labs PlanetScope (3.7–4.1 m resolution) to identify areas of concern, such as landslide dammed streams or impacts to transportation networks. Draft landslide point data and daily summaries were shared with responding humanitarian agencies within days of the earthquake and a final point inventory covering the entire landslide affected region was released online on December 7, 2021 (Martinez et al. 2021b).

1.2.2 NASA disasters program

The NASA Disasters Program supports a range of research activities to further our understanding of hazards and utilizes diverse satellite and ground-based products to analyze the interconnected hazards, exposure and risk for a variety of disaster types. Alongside research activities, the program also supports response activities during disasters, enabling stakeholders to use NASA data and analysis to improve decision-making before, during and after disasters. The program engages with over 300 different stakeholders around the world and provided monitoring and support for over 40 different disasters in 2021. In the same period, over 2000 different data consisting of imagery and derived products were uploaded to the NASA Disasters Portal to support stakeholders in an open-science framework.

During earthquakes and tropical storms, the program utilizes satellites run by government agencies to conduct optical- and Synthetic Aperture Radar (SAR)-based analysis of ground surface deformation, flooding, and landslides, among other hazards. In addition, very-high (< 3 m) and high (3–5 m) resolution commercial imagery is also used as part of the ongoing Commercial Satellite Data Acquisition (CSDA) Program agreement between NASA and several private sector entities. These datasets provide overlapping and complementary hazard analysis.

A combined team consisting of the NASA Goddard Space Flight Center landslides team, NASA Jet Propulsion Laboratory, University of California Los Angeles, and University of Maryland (herein NASA Team) responded on behalf of the NASA Disasters program. The primary landslide product produced and shared by the NASA Disaster program was the Sentinel-1 SAR backscatter change-based landslide heatmap (Handwerger et al. 2022) covering the entire landslide affected area. A heatmap is a data visualization technique that highlights landslide areas by calculating the density of potential landslide pixels in a location over a given radius using a kernel density estimation. Publication occurred two days following the earthquake on August 17, 2021.

1.2.3 International charter space and major disasters

The International Charter Space and Major Disasters (herein the Charter) is a worldwide collaboration of space agencies and disaster monitoring organizations, through which satellite datasets are made available for the benefit of disaster management. By combining Earth observation (EO) assets from 17 space agencies, the Charter allows resources and expertise to be coordinated for rapid response to major disaster situations, thereby helping civil protection authorities and the international humanitarian community. The Charter allows mobilizing agencies around the world and benefits from their know-how and their satellites through a single access point that operates 24 h a day, 7 days a week and at no cost to the user. Since 1999, 770 activations over 131 countries allowed the documentation of many geohazard and hydro-meteorological hazard disasters. Two organizations responded on behalf of the Charter upon request from the Charter member representing their country (French Space Agency and Canadian Space Agency).

1.2.3.1 Centre national de la recherche scientifique (CNRS): ecole et observatoire des sciences de la terre (EOST)

In Haiti, the Charter was activated on August 14, 2021 (ID729, https://disasterscharter.org/en/web/guest/activations/-/article/earthquake-in-haiti-activation-729-) to evaluate impacts from the earthquake. Activation was initiated by three entities (the Direction de la Protection Civile de Haiti, the United Nations Operations and Crisis Center (UNOCC), and the Centro Nacional de Prevención de Desastres (CENAPRED) in Mexico) in order to support international aid and help with disaster recovery through dispatch of supplies and personnel. The French Space Agency (CNES) coordinated the provision of satellite imagery for the activation and requested the support of research centers to create value-added products. In this context, Centre National de la Recherche Scientifique (CNRS) through its research institute École et Observatoire des Sciences de la Terre (EOST) (herein CNRS-EOST) was asked by CNES to generate rapid landslide maps containing information on the spatial location and the size of the landslides for the territory. Specifically tasked, very-high (< 3 m) resolution multispectral imagery from Airbus Pléiades was used for mapping landslides using an automated supervised image classification service called ALADIM ( Automated Landslide Detection and Inventory Mapping; Déprez et al. 2022). The Pléiades imagery was acquired through the CIEST (https://www.poleterresolide.fr/ciest-2-nouvelle-generation-2/) (Cellule d’Intervention d’Expertise Scientifique et Technique) of the French Solid Earth datahub ForM@Ter (https://en.poleterresolide.fr/) with the support of CNES. A draft landslide map in polygon format covering the central zone of the Pic Macaya National Park area was produced and delivered to CNES on September 2, 2021. An updated version of the map was delivered again on October 14, 2021.

1.2.3.2 BGC engineering

BGC Engineering (herein BGC) is a Canadian-based geotechnical engineering firm with a specialization in geohazards and remote sensing. The Canadian Space Agency (CSA), which is a member agency of the Charter, requested assistance from BGC in responding to Charter activation ID730 (https://disasterscharter.org/web/guest/activations/-/article/flood-large-in-haiti-activation-730-) on August 17, 2021. This activation followed the passage of Tropical Storm Grace over western Haiti. The activation intended to provide situational awareness to responders regarding flash flooding, riverine flooding, and landslides as a result of heavy rains. Imagery was provided by the Korean Multi-Purpose Satellite (KOMPSAT) and Canadian RADARSAT constellation mission (RCM). BGC produced flooding inundation maps from backscatter analysis on the RCM products. The 0.5 m resolution KOMPSAT imagery was focused on the area surrounding the Pic Macaya National Park, the zone of most dense landslides. Due to the single time step of the imagery (August 17, 2021) and its limited spatial coverage, BGC also utilized Sentinel-2 imagery for broad-scale landslide scar mapping. The mapping was performed with the Google Earth Engine (GEE) application HazMapper, which is targeted at landslide identification from multispectral satellite imagery (Scheip and Wegmann 2021). Manual interpretation of imagery from KOMPSAT was used as a first order estimate on accuracy of the Sentinel-2 based products. These Sentinel-2 based products were delivered as an estimated landslide scar map covering the entire landslide affected area on September 3, 2021.

2 Data

The rapid response products were derived from both optical and SAR imagery.

2.1 Optical imagery

Optical imagery from different satellites with varying degrees of coverage and cloud cover were used by the different groups. Freely accessible medium (5–10 m) resolution imagery, and commercial very-high resolution (< 3 m) and high-resolution (3–5 m) imagery provided for the Charter was accessed using the governmental data access webtool known as the USGS Hazard Data Distribution System (HDDS):

• Copernicus Sentinel-2 images (10 m) covering the full area of interest.

• Airbus Pléiades imagery (multispectral image at 2 m resolution, and a panchromatic image at 0.5 m resolution) acquired through the Charter. A cloud-free image over the central zone of the Pic Macaya National Park covering a 385 km² area was acquired on August 23, 2021.

• Maxar WorldView images (0.3–0.5 m) covering only the surrounding area of the main affected zone. High cloud cover prevented a full exploitation of the images and regular scattered clouds allowed mappers to use windows in cloud cover to map landslides, so we used any available imagery starting on August 14, 2021.

• Planet Labs PlanetScope images (3.7–4.1 m) covering the main affected zone and surroundings. Mappers used windows in cloud cover to map landslides.

2.2 SAR imagery

SAR data from the C-band (~ 5.6 cm wavelength) Copernicus Sentinel-1 satellite constellation archived on GEE was used to detect landslides. GEE provides Level-1 Sentinel-1 Ground Range Detected (GRD) backscatter intensity coefficient (σ∘) images with a 10 m pixel spacing that have been corrected for thermal noise, calibrated in terms of radiometry but are not terrain flattened. Data collected in vertical transmit/horizontal receive (VH) polarization were used. To reduce noise from poor-quality data, all pixels with values ≤  − 30 dB were removed (based on a recommendation from the GEE S1 Data Catalog).

3 Methods

In the aftermath of the Nippes earthquake and Tropical Storm Grace, four rapid response products were released. These four products were created using several methodologies, and released at different times (Table 1). USGS used manual mapping to generate an initiation point-based landslide inventory covering the entire landslide affected area. NASA, CNRS-EOST and BGC used automated approaches to produce landslide inventories/maps offering varying coverage of the landslide affected area (Fig. 2).

Fig. 2

Area covered by different automated rapid response landslide products

3.1 USGS/LDAT manual mapping

A point-based landslide inventory was initiated for the entirety of the Tiburon Peninsula within hours of the release of the first post-earthquake optical imagery (herein USGS inventory). USGS analysts used imagery from Sentinel-2, WorldView, and PlanetScope as it became available to map landslide points at the interpreted initiation area of each landslide. Landslides were mapped in GIS using two methods. The first involved downloading georeferenced imagery from HDDS and the second used GIS tools by Maxar and Planet that allow limited access to each company’s respective imagery through an online viewing tool. Analysts used a variety of platforms to view imagery and map landslides, including ArcMap, ArcGIS Pro, and QGIS. Analysts did not differentiate landslides triggered by the earthquake from those triggered by rainfall because imagery was not acquired over much of the affected area in the short time period between the two events. Therefore, landslides from the tropical storm were likely captured in the inventory but the group noted that few landslides had characteristics of rainfall-triggered landslides and they thought most were triggered by the earthquake (Martinez et al. 2021a). For each point, analysts assigned a certainty to the presence of the landslide, which was based on analyst observations of imagery and comparisons to lidar digital terrain model (DTM) slope shades to ensure the landslide mobilized downslope. Analysts also identified potential hazards such as a landslide dam in a stream channel, impacted transportation networks, and impacted structures. Because of cloud cover, poor quality imagery, lighting issues, image distortion due to rectification issues, and varying resolution imagery, analysts attempted to evaluate the same area several times to capture the greatest population of landslides. Both the spatial dataset and observations were documented in a daily landslide situation report (sitrep) documenting areas of potential concern, such as landslide dams with observed upstream ponding or landslides impacting transportation networks. The sitrep was shared with the USAID, USGS, other responding United States Government (USG) partners, and non-governmental humanitarian organizations.

3.2 Automated landslide detection and mapping

Three automatically generated landslide rapid response products were created and disseminated:

• NASA released a SAR backscatter change-based landslide heatmap covering the entire landslide affected area one day after the earthquake using Sentinel-1 imagery and a method developed by Handwerger et al. (2022).

• CNRS-EOST released a polygon-based landslide inventory after the earthquake (September 2, 2021) for the central zone of the Pic Macaya National Park area produced using Pléiades imagery and the machine-learning detection method ALADIM (Déprez et al. 2022).

• BGC used HazMapper (Scheip and Wegmann 2021) to generate a relative difference Normalized Difference Vegetation Index (rdNDVI) change-based landslide map covering the entire landslide affected area using Sentinel-2 imagery. BGC released the product to the Charter on September 3, 2021.

The automated products were released without any manual corrections. The three methods used to generate landslide products are described in detail in the following sections.

3.2.1 SAR backscatter change

Satellite SAR data can be used to rapidly detect landslides, often within days of a triggering event, because it penetrates clouds and operates day and night. The ability of SAR to penetrate clouds is particularly important for rapid landslide response as cloud cover can hinder landslide detection from optical imagery for days to weeks, which was the case for this event in Haiti (Martinez et al. 2021a; Whitworth et al. 2022; Zhao et al. 2022; Havenith et al. 2022). Landslides can cause measurable backscatter changes because they alter the ground surface properties. For example, landslides cause changes in hillslope geometry and surface roughness, both of which strongly affect the radar backscatter (Rignot and van Zyl 1993; Adriano et al. 2020; Mondini et al. 2021). Thus, changes in backscatter from before and after a triggering event can be used to detect landslide areas.

For rapid response to the Nippes earthquake-triggered landslides, NASA used a Sentinel-1 SAR backscatter change approach developed in GEE (Handwerger et al. 2022) to detect areas with high landslide density (herein SAR heatmap). This SAR heatmap was also used as a case study in Handwerger et al. (2022). This approach detects potential landslides by calculating the change in the backscatter coefficient before and after the triggering event using the standard log ratio change approach, I_ratio, defined as I_ratio = 10 × log₁₀(I_pre/I_post) (e.g., Mondini et al. 2019; Jung and Yun 2020). Because GEE provides the backscatter coefficient data in dB, the backscatter change was calculated as I_pre − I_post which is equivalent to the I_ratio (logarithm quotient rule).

False positives that may occur in flatter areas such as backscatter change due to flooding or agriculture was removed by using topographic slope from the 1 arcsec (∼30 m) resolution NASADEM (NASA JPL 2020) to mask areas with slopes < 10°. To further reduce noise and improve landslide detection, the team calculated the pre-event backscatter by taking the median of 104 pre-event images acquired between August 1, 2017, and August 13, 2021. This stack was differenced from a single post-event image acquired one day after the earthquake (August 15, 2021) to calculate the backscatter change map. Using stacks of SAR data reduces noise, and furthermore, the pre-event stack provides backscatter data that is more representative of the pre-event ground surface properties. Finally, to detect landslide areas, the backscatter change raster was thresholded using the 99th percentile to highlight strong signals of changes, and the heatmap was calculated. An example web URL for conducting this analysis in GEE is provided in the data availability section.

3.2.2 HazMapper

HazMapper (https://hazmapper.org/) is an open-access application that provides pixel-based multispectral analysis capability in GEE (Scheip and Wegmann 2021). Normalized Difference Vegetation Index (NDVI) is a multispectral index that is well suited to identifying the relative health of green vegetation (Tucker 1979). Landslides result in a decrease in vegetation cover compared to pre-landslide conditions. As such, change in NDVI can be used to detect landslides (Mondini et al. 2011; Lu et al. 2019; Martinez et al. 2021c). HazMapper is designed to compute the relative difference in the NDVI (rdNVDI) across a pre- and post-event time period bracketing a time of interest (e.g., a date of a suspected landslide event):

$$rd{\text{NDVI}} = \left( {{\text{NDVI}}_{{{\text{post}}}} - {\text{NDVI}}_{{{\text{pre}}}} {/}\sqrt {{\text{NDVI}}_{{{\text{pre}}}} + {\text{NDVI}}_{{{\text{post}}}} } } \right) \times 100$$

(1)

The method has shown success at locating landslides, particularly when landslides exceed tens of meters in width or length (Scheip and Wegmann 2021; Abancó et al. 2021; Li et al. 2022; Lindsay et al. 2022; Wooten et al. 2022). The HazMapper source code for this analysis is provided in the data availability section. To fully recreate this landslide inventory, users must download the rdNDVI results and apply a rdNDVI threshold (e.g., − 25%).

HazMapper’s rdNDVI metric was employed for the western Tiburon Peninsula in Haiti. Table 2 lists the input parameters used in analysis. Analysis used 2-weeks of available imagery following the date of the earthquake and provided a rapid assessment of ground conditions. Significant cloud cover remained over the region, unsurprising given the time of year and recent passage of a major storm. The 2-week look window included four images used in the post-event period.

Table 2 HazMapper input parameters

Because the Charter requested map deliveries of likely landslide locations across western Haiti, BGC provided a continuously scaled − 5 to − 100% rdNDVI output. This permitted better visualization of landslide locations even though at the pixel scale, provided significant false positives. The objective of this delivery was not to identify individual landslides, but instead to highlight areas with the most dense landslides as a result of the earthquake and tropical storm. Results were uploaded as a PDF to the Charter’s activation page.

Previous work on landslide mapping in the tropics has suggested a − 30% threshold adequately estimated landslide locations in HazMapper-derived results (e.g., Wegmann et al. 2021). This rdNDVI threshold was used to generate landslide polygons needed for validation effort in this study.

3.2.3 ALADIM

ALADIM (Déprez et al. 2022) combines object-based image analysis (OBIA) and machine learning to detect and map landslides triggered by large forcing events such as earthquakes and heavy rainfall. Instead of working with individual pixels, OBIA groups similar pixels to form homogenous objects using image segmentation (Blaschke 2010). Landslide diagnostic features of these objects (e.g., NDVI, brightness, slope etc.) are used to detect landslides using supervised machine learning. ALADIM can process optical imagery of medium (typically Landsat 8/9 or Sentinel-2), high (PlanetScope) and very-high (typically Pléiades, WorldView) resolution and requires a small dataset of manually mapped landslides to train the machine learning model. ALADIM is an online web service accessible on the Geohazards Exploitation Platform (GEP, https://geohazards-tep.eu/) and is developed, maintained and exploited by CNRS-EOST (Strasbourg, France). The processing chain can be used in two modes, using either pre- and post-event images (change detection mode), or using only a post-event image (monodate mode).

ALADIM builds on the OBIA and machine learning approach described in Stumpf and Kerle (2011) and Stumpf et al. (2013). The processing chain consists of applying Mean-Shift segmentation (Michel et al. 2014) on the post-event imagery to generate objects. Next, high-level image (e.g., spectral, texture, shape) and topographic (using ancillary Digital Elevation model (DEM)) features of these objects are calculated. A training sample is generated by interpreting, digitizing and labeling a few landslides on the post-event imagery. Utilizing this training sample, a Random Forest classifier (Breiman 2001) is trained and objects are classified as belonging to landslide or non-landslide class. The output consists of a polygon shapefile with probability of each object belonging to a landslide class.

In the Haiti rapid response, ALADIM was deployed in monodate mode using a Pléiades satellite image acquired on August 23, 2021. The input image was a pan-sharpened image at 0.5 m resolution. The Copernicus GLO-30 DEM with a resolution of 30 m was used to calculate the topographic features. Mean-Shift segmentation was performed with a spatial radius of 5 pixels, range radius of 15 and minimum segment size of 1000 pixels. The minimum segment size sets the smallest area of landslide to be mapped. Based on a single pixel area of 0.25 m², the minimum segment size of 1000 pixels will identify landslides with areas of 250 m² or greater. Sixty landslide diagnostic features (11 spectral, 8 textural, 9 morphometric, 15 histogram-oriented gradient and 17 topographic) were computed for use in landslide detection. A training sample of 394 landslides was created manually over a very small portion of the affected region (4 km² out of a 385 km² area of the full image). This sample was used to train the Random Forest classifier to predict landslides over the complete coverage offered by the Pléiades imagery (green box in Fig. 2). The ALADIM setup and result can be accessed in GEP using the link provided in the data availability section.

3.3 Validation of rapid response products

The rapid landslide inventories/maps discussed herein were delivered with haste to government and non-governmental humanitarian organizations and/or to the public within days to weeks of the event. It is important for scientists to understand both the accuracy and the limitations of various methods used for providing these rapid inventories. As such, here we validate the methods using a reference manual inventory generated within the months following the event. This collaborative validation exercise amongst different agencies should increase our collective ability to more effectively map landslides following future events. Our validation approaches and the manual reference inventory are described in the sections below.

3.3.1 CNRS-EOST reference manual landslide inventory

A polygon-based manually mapped landslide inventory was created using satellite imagery from Pléiades, WorldView and Sentinel-2 (herein CNRS-EOST reference inventory). The coverage from different satellite imagery is shown in Fig. 3. Priority was given to mapping using imagery from Pléiades and WorldView. The Pléiades imagery was distributed as ortho-rectified products by Airbus, and the WorldView images were ortho-rectified at first order using a lidar DTM (2 m resolution) acquired in 2016 (The World Bank 2021). Sentinel-2 was used in areas that were not covered by the commercial imagery. Relevant indicators, easily distinguishable in the imagery were used to detect the presence of landslides using a pre/post-event approach (e.g., changes in land use marked by changes in band radiometry, variation in vegetation density, presence of bare soils, presence of specific landslide features such as scarps, cracks, and surface disturbances). Interpretation of the imagery was also supported by the analysis of the 2016 lidar DTM to identify if the possible landslides are located on slopes > 10°. The image interpretation on Pléiades and WorldView imagery was carried out at an interpretation scale of 1:1000 allowing detection of landslides of minimum size of 10 m². A total number of 14,482 landslides corresponding to a total area of 45 km² was mapped (Fig. 4). The CNRS-EOST reference inventory can be downloaded from the Cooperative Open Online Landslide Repository (Juang et al. 2019; see data availability section).

Fig. 3

Map illustrating the extents of the satellite imagery used for constructing the CNRS-EOST reference inventory

Fig. 4

a Spatial distribution of landslides present in the CNRS-EOST reference inventory. b Zoomed in view of blue box in a, showing mapped landslides. (Pléiades ^©CNES 2021, Distribution Airbus DS)

During the preparation of this manuscript, two independent studies regarding the landslides triggered by the 2021 Haiti events were published by Havenith et al. (2022) and Zhao et al. (2022). Both of these studies created landslide inventories but only Zhao et al. (2022) shared the inventory publicly. To compare CNRS-EOST reference inventory to Zhao et al. (2022) inventory, we compute the frequency-area distribution (FAD; Fig. 5) for each inventory. The FAD follows characteristics observed by others, namely a logarithmic decay for medium and larger landslides and an apparent lack of smaller landslides (Malamud et al. 2004; Tanyaş et al. 2019). Havenith et al. (2022) also computed a FAD; however, the inventory was recognized to be incomplete for landslides below a size of 3,000 m².

Fig. 5

Frequency-Area distribution (FAD) plot of the CNRS-EOST reference inventory and the Zhao et al. (2022) landslide inventory

Manual mapping of landslides is a subjective process and very dependent on the type of imagery used, its resolution, the time of capture, the procedure used for the mapping process and the experience of the landslide experts (Hölbling et al. 2017; Fan et al. 2019). Based on comparisons of available FADs for this event (Fig. 5), it is apparent that the CNRS-EOST reference inventory captures more smaller landslides than Zhao et al. (2022) inventory. The FAD is based on landslide polygons and so the methods around drawing boundaries between landslides or combining multiple landslide polygons into a single polygon during mapping can increase difficulty in comparing inventories based on FADs alone. Zhao et al. (2022) reports more larger landslides than CNRS-EOST, however, we hypothesize this is due to combining multiple landslides into a single polygon as it is unlikely that CNRS-EOST missed landslides exceeding 100,000 m² in area. The CNRS-EOST inventory also identifies smaller landslides (< 50 m²) and also identifies more landslides than Zhao et al. (2022) in the approximately 50 m² to 2,000 m² range.

3.3.2 Validation strategy

Different remote sensing techniques (e.g., pixel vs. OBIA) and outputs (polygons vs. heatmap) make it both difficult and impractical to focus on a single validation method. The most utilized validation method for polygon-based output is calculation of accuracy metrics based on overlapping areas/number of landslides (Martha et al. 2010; Rau et al. 2014; Hölbling et al. 2015, 2016). However, this overlapping area/number analysis cannot be done for SAR heatmap as the outputs are not available in polygon format. Hence, we utilized the two validation approaches described below.

3.3.2.1 Comparison using accuracy metrics

The CNRS-EOST reference inventory was used to validate the landslides detected by ALADIM and HazMapper. The validation was carried out considering complete coverage offered by the individual products (Fig. 2). Three metrics based on overlapping areas of polygons were calculated: true positive (TP), false negative (FN), and false positive (FP). True positives are areas where landslides were detected by automated mapping methods and were present in the reference inventory. False negatives are areas where the automated mapping methods did not identify a landslide, but a landslide was noted in the reference inventory. False positives are areas where automated mapping methods reported a landslide, but a landslide was not present in the reference inventory. Based on these metrics, the two accuracy indices, Producer Accuracy (PA) and User Accuracy (UA) were calculated. The PA denotes how much of the manual inventory was detected. The UA denotes how much of the detected landslides are actual landslides.

$${\text{PA}} = \frac{{{\text{TP}}}}{{{\text{TP}} + {\text{FN}}}}$$

(2)

$${\text{UA}} = \frac{{{\text{TP}}}}{{{\text{TP}} + {\text{FP}}}}$$

(3)

Validation based on the number of landslides would be fraught with error due to difficulties in the ability of automated methods to determine when to combine or separate amalgamated landslides (Marc and Hovius 2015). Hence, number-based validation was not done.

3.3.2.2 Comparison using density maps

We also used a density map approach to validate the landslide products, especially in the context of a rapid response, where it can be useful to rapidly identify areas significantly impacted by landslides. This density map-based approach is just used for visual/qualitative validation and not shared during rapid response. The objective of this approach is not to compare landslide detection methods at the pixel or landslide scale, but to evaluate how effective the methods were at detecting areas with the highest landslide density, which may identify the most severely landslide-impacted regions following catastrophic events.

For ALADIM and HazMapper we used a density map approach to validate results using a regular grid (500 m pixel spacing) for the complete coverage offered by each method. We compute the percentage coverage of each grid cell by landslides. This illustrates the relative density of landsliding determined by each method and allows comparison with the density map generated using CNRS-EOST reference inventory.

We also compared the SAR heatmap to the CNRS-EOST reference inventory. For SAR data validation, we created a new heatmap using the CNRS-EOST reference inventory to more closely reflect the method used by Handwerger et al. (2022). The CNRS-EOST reference inventory was converted to a raster of 10 m cell size. Each cell in the landslide raster was then converted to points for generation of the heatmap.

4 Results

4.1 USGS inventory

The USGS inventory covered the entire Tiburon Peninsula. A total of 4,893 landslides were mapped during the rapid response effort (Fig. 6). USGS analysts were also able to simultaneously identify hazardous areas primarily landslide dams and damage to roads.

Fig. 6

a Spatial distribution of landslides present in the USGS inventory. b Zoomed in view of blue box in a, showing mapped landslides. (Pléiades ^©CNES 2021, Distribution Airbus DS)

4.2 SAR heatmap

The SAR heatmap covered a 12,126 km² area of the Tiburon Peninsula and was only able to identify areas of dense landslide activity around the Pic Macaya National Park (Fig. 7). It was not able to highlight areas that had less or individual landslides.

Fig. 7

SAR heatmap produced one day after the Nippes earthquake with the CNRS-EOST reference inventory overlaid. Heatmap made from stacks consisting of 104 descending pre-event Sentinel-1 SAR images collected between August 1, 2017 and August 13, 2021 and 1 descending post-event SAR image collected on August 15, 2021

4.3 HazMapper-generated landslide map

The HazMapper analysis covered a 4,940 km² area of the Tiburon Peninsula (Fig. 8) and identified 96,238 landslides with a total area of 68.28 km².

Fig. 8

a Landslides detected by HazMapper. b Zoomed in view of blue box in a, showing mapped landslides overlaid over the Pléiades imagery illustrating the 2-week look during the rapid response effort. Note the large areas of landslides visible in the satellite imagery but absent in the HazMapper inventory. Many of those areas were covered by clouds in the Sentinel-2 imagery. (Pléiades ^©CNES 2021, Distribution Airbus DS)

4.4 ALADIM generated landslide map

ALADIM covered a 385 km² area around the central zone of the Pic Macaya National Park area (Fig. 9) and identified 11,114 landslides with a total area of 43.93 km².

Fig. 9

a Landslides detected by ALADIM. b Zoomed in view of blue box in a, showing mapped landslides overlaid over the Pléiades imagery. ALADIM only covered the central zone of the Pic Macaya National Park area. (Pléiades ^©CNES 2021, Distribution Airbus DS)

4.5 Validation

4.5.1 USGS versus CNRS-EOST reference inventory

The USGS inventory is in point format and the CNRS-EOST reference inventory in polygon format, we used intersection of points with polygons to quantify the difference between these two inventories. The USGS inventory had 1772 landslides that were not present in the CNRS-EOST reference inventory. Moreover, 11,463 landslides present in the CNRS-EOST reference inventory were missing from the USGS inventory. There were 3121 common landslides present in both the inventories. Interpretation of satellite imagery and the difficulty in identifying boundaries between adjacent landslides can make direct numerical comparison of the number of landslides challenging. Therefore, these values should be interpreted cautiously.

4.5.2 HazMapper and ALADIM accuracy metrics

For HazMapper, we calculated accuracy indices twice, with and without cloud cover mask. Cloud cover in Sentinel-2 imagery prevented detection of landslides for at least two months. Hence, we masked CNRS-EOST reference inventory with Sentinel-2 cloud mask and calculated the accuracy indices again (HazMapper with cloud mask).

ALADIM did not map landslides smaller than 250 m² due to the minimum segment size (1000 pixels) used during segmentation. Hence comparison was done only for landslides greater than 250 m². ALADIM was able to map 76.9% of landslides present in CNRS-EOST reference inventory and 52.4% of polygons mapped by ALADIM were actual landslides (Table 3). The reasonable accuracy of ALADIM, in part, stems from the fact that the imagery used to detect landslides by ALADIM is the same as that used for the reference manual mapping. It is, in a sense, a true validation because conditions in the imagery are identical. For HazMapper (which uses Sentinel-2 imagery) we only compared landslides greater than 100 m² to compensate for the 10 m resolution. The Sentinel-2 imagery had very different conditions in terms of area visible due to shadows and cloud cover. As a result, our findings should be considered a relative validation. When validating without accounting for cloud cover present in the Sentinel-2 imagery, the PA and UA of HazMapper was 40.9% and 27.1%, respectively (Table 3, HazMapper). In order to conduct validation taking into consideration actual visible landslides in Sentinel-2 imagery, we masked the CNRS-EOST reference inventory using a cloud mask generated from the Sentinel-2 imagery. The PA improved to 59.8% when we did not consider landslides obstructed by clouds (Table 3, HazMapper with cloud mask).

Table 3 Accuracy indices of ALADIM and HazMapper (polygon-based)

Visual analysis of TP, FP and FN over the central zone of the Pic Macaya National Park (area of dense landslide activity) suggests that ALADIM detects a large number of landslides but also includes large numbers of FP (Fig. 10b). In the case of HazMapper, we observe evidence of large amounts of FN which is largely due to clouds covering large areas where landslides were present in the Sentinel-2 imagery (Fig. 10a). We do observe fewer FP in HazMapper when compared to ALADIM over the Pic Macaya National Park area.

Fig. 10

True positive, false positive and false negative visualized over the central zone of the Pic Macaya National Park area (green box in Fig. 2) for illustrative purposes: a HazMapper and b ALADIM. These metrics were calculated over the complete area covered by each method (Fig. 2). Areas not colored green, yellow, or red are true negatives. (Pléiades ^©CNES 2021, Distribution Airbus DS)

4.5.3 Comparison of density maps

The 500 m regular grid-based density map for HazMapper (Fig. 11) and ALADIM (Fig. 12) both capture the dense landsliding in the center of the peninsula and near Pic Macaya National Park. In comparison to the CNRS-EOST reference inventory density map, the automated density maps show larger areas with presence of landslides. Grids with 1–5% landslide coverage are higher than that of the CNRS-EOST reference inventory density map. The dense landslide areas with landslide coverage of 5% or more are easily identifiable and comparable between the density maps.

Fig. 11

Density map of HazMapper relative to the CNRS-EOST reference inventory in a 4940 km² area (red polygon in Fig. 2)

Fig. 12

Density map of ALADIM relative to the CNRS-EOST reference inventory in a 385 km² area (green polygon in Fig. 2) (Pléiades ^©CNES 2021, Distribution Airbus DS)

Overall, the SAR heatmap worked well for rapid response, capturing the area of highest landslide density within a day of the earthquake, but was not useful in identifying areas with low landslide density. The SAR heatmap does highlight larger areas (red blob; Fig. 13) with high landslide density due to false positives generated during the SAR backscatter change landslide detection process.

Fig. 13

SAR heatmap relative to the CNRS-EOST reference inventory in a 12,126 km² area (yellow polygon in Fig. 2)

5 Discussion

Williams et al. (2018) documented experiences and lessons learned conducting satellite-based emergency landslide detection and mapping using optical imagery for the 2015 Mw 7.8 Gorkha earthquake, Nepal. They found forms of landslide assessment (primarily manual mapping) to be too slow to generate relative to the speed of a humanitarian response and hinted that the use of automated landslide mapping techniques and SAR information might help accelerate this assessment. However, they reported application of semi-automated and automated mapping with variable image characteristics and its quality for humanitarian response was not well documented. Seven years after the Gorkha earthquake, the 2021 Haiti event presents an opportunity to validate different automated and manual products utilizing both optical and radar sensors.

5.1 USGS inventory

The USGS inventory, comprising 4893 initiation points, was shared publicly in preliminary form within a few days of the earthquake and a final version was released formally in December 2021 (Martinez et al. 2021a). Comparison of the USGS inventory with the CNRS-EOST reference inventory revealed both inventories contained different numbers of landslides.

Available WorldView imagery sometimes exhibited poor georeferencing and orthorectification, precluding an accurate spatial location for landslides. These landslide initiation points were generally within tens of meters of the actual landslide location, sometimes placing these points outside of reference inventory polygons and vice versa. For the purpose of the rapid response, these were useful information, however, they make quantitative comparison difficult. In addition, comparison of two inventories produced using different satellite imagery is subjective due to difference in spatial resolution, amalgamation and area visible in satellite imagery due to shadow and off-nadir angle of satellite (Meena and Tavakkoli Piralilou 2019).

Although labor intensive, there were several benefits of the manual inventory for an emergency response. The first was the rapid identification of interpreted landslide-related hazards including landslide dammed streams, blocked transportation networks, and damaged buildings. Analysts could quickly identify potential life/safety hazards and flag those impacts in the daily sitrep for first responders. Additionally, following the earthquake, clouds persisted for weeks over much of the Tiburon Peninsula; however, mappers could identify landslides through windows in the clouds, eventually mapping much of the Peninsula in a little over a week. Finally, analysts mapped landslides by systematically panning optical imagery and identified other geologic hazards such as liquefaction and stream-adjacent lateral spreads that were often identified by a subtle headscarp and a toe in the stream.

5.2 SAR heatmap

Under cloud free scenarios, the utility of semi- or fully-automated mapping from optical imagery is unquestionable. However, cloud cover can prevent rapid landslide mapping from optical imagery for weeks to months. For the 2021 Nippes earthquake- and Tropical Storm Grace-triggered landslide events, cloud cover hindered clear observations of the ground for over two months. This scenario is common when attempting to map rainfall triggered landslides in tropical environments. Landslide detection using SAR helps circumvent issues related to clouds and is well suited for rapid response. For the Nippes earthquake event, a single post-event descending pass image from the Sentinel-1 satellite collected one day after the earthquake enabled generation of the first estimate of landslide affected area. The SAR heatmap highlighted dense landslide activity around the Pic Macaya National Park. However, there are limitations when using SAR data. First, analysis and interpretation of SAR data often requires specialized training and software. GEE has eliminated some of these barriers by making the data easier to access, process, and analyze. SAR data are also limited by geometric shadow and distortion (Adriano et al. 2020; Mondini et al. 2021), particularly in steep mountainous terrain, and by speckle noise. To reduce speckle noise smoothing filters can be applied, but this degrades the spatial resolution of the imagery.

The SAR heatmap method uses the 99th percentile threshold to detect landslide pixels from backscatter change raster. Handwerger et al. (2022) showed that the landslide areas become pronounced as more post-event images are added to the backscatter change analysis. However, in rapid response scenarios the number of post-event images is typically limited as it can take several days for repeat open access satellite acquisitions. The thresholded backscatter change by itself may not be helpful at identifying individual landslides but could provide meaningful information of areas with the most significant change. Hence, a heatmap approach which is essentially a landslide density estimation was used. In Haiti, the Pic Macaya National Park area was clearly highlighted, and had the most landslides, whereas the areas with less landslide activity, and smaller landslides, to the East were not highlighted by the heatmap approach. We note that these smaller landslides were also not clearly observed in the heatmap made using the CNRS-EOST reference inventory (Fig. 13).

Although a topographic mask was applied to remove FPs, there is still the possibility of real surface changes related to natural and human-made causes, including but not limited to mining, construction, deforestation, agriculture, flooding, and changes in reservoir water levels. Hence the SAR heatmap for rapid response scenario should be interpreted with caution. We suggest use of other ancillary data such as landslide susceptibility and recent land cover maps to identify changes are happening in regions where landslides are expected. Gaps in cloud cover in the optical imagery can also be used to interpret if areas highlighted by the SAR heatmap are associated with landslides.

The Sentinel-1 SAR constellation had 2 satellites offering 6 days revisit time combined and 12-day revisit time individually. On December 23, 2021, Sentinel-1B satellite experienced an anomaly and has been de-commissioned by ESA as of March 2022. For this event, there was an overpass right after the earthquake on August 15, 2021, enabling generation of a heatmap. However, for future events the longer revisit time might cause a delay in production of the heatmap until the launch of Sentinel-1C in the near future. Commercial SAR companies (e.g., ICEYE, Capella) can acquire repeat imagery over the same area maintaining precise satellite orbit every 24 h. This extremely rapid revisit time can significantly increase the capability to rapidly detect landslides and provide key information to first responders.

5.3 HazMapper

We used rdNDVI threshold to convert the HazMapper results into polygon-based inventory for validation and further analysis. Selecting a rdNDVI threshold will be event and site specific and based on end user’s preference. Some may want to map all landslides with risk of including high FPs whereas others might want to generate less FPs at the risk of omitting some landslides. Interpretation of the accuracy metrics (Sect. 4.2.1) generated in comparison with the CNRS-EOST reference inventory helps to put the FN (missed detections) and FP detections into context.

False negatives were primarily identified either where cloud cover obscured analysis, or for landslides < 100 m² (the pixel size of HazMapper results). HazMapper returns “no data” pixels where clouds were present during analysis, helping to differentiate between false and true negatives. False positives were typically present in drainage channels where debris flows and debris floods were likely occurring or where unidentified sources of noise or anthropogenic activity (e.g., development or deforestation) caused losses in NDVI. Because the CNRS-EOST reference inventory did not aim to identify mass wasting processes in channels (e.g., debris flows/floods), these apparent FPs in the HazMapper result may be due to a difference in objective (e.g., identifying all slope failure processes types compared to only hillslope landslides).

Following product delivery, we evaluated percentage change in HazMapper rdNDVI against the CNRS-EOST reference inventory. The highest agreement between the initial HazMapper inventory and the reference inventory was noted by using an rdNDVI threshold of − 25%. Using this threshold, landslide pixels were then polygonized to provide landslide polygons for analysis. In addition, we updated the analysis with a 2-month look after the August landslide event. This illustrates how open access imagery-based results can improve through time due to reduced cloud cover. Moving from a 2-week window to a 2-month window increases the agreement in the reference grid coverage between HazMapper and the CNRS-EOST reference inventory (Fig. 14c). A slight decrease in the PA (Fig. 14a) illustrates the omission of some TPs identified in the shorter look window. This occurs likely due to vegetative regrowth in the 2-month timeframe following the landslide event coupled with HazMapper’s comparison method, which compares single images representing the peak phenological cycle for pre- and post-event time windows. As post-event windows increase in duration, any vegetative regrowth will no longer be identified as a landslide, causing a reduction in TP identification through time. An increase in the UA (Fig. 14b) illustrates the reduction in FPs and is of a higher magnitude than changes to the PA. The reduction in the FP rate is likely due to a reduction in low probability cloud cover that was used in the original analysis (pixel-based cloud probability threshold of 30%). An updated best-fit agreement of the grid-based approach (Fig. 14c) illustrates closer overall results between the HazMapper and CNRS-EOST reference inventory when a 2-month window is used, however, the use of 10 m imagery will always preclude the identification of smaller landslides, and so in general, the grid-based curve for HazMapper products tends to underreport landslide coverage areas relative to the CNRS-EOST reference inventory.

Fig. 14

Validation results for varying look-windows of HazMapper. Analyses follow the parameters in Table 2 with the exception of the post-window (2 weeks or 2 months): a Producer accuracy (PA) for both analyses. b User Accuracy (UA) for both analyses and c Agreement between percentage of 500 m grid cells covered by landslides between HazMapper and the CNRS-EOST reference inventory

5.4 ALADIM

Errors observed in the ALADIM product can, as for HazMapper, be categorized as FNs and FPs. False negatives were present in areas affected by clouds and/or by bright radiometry due to near-vertical sun illumination, and for landslides < 250 m² (due to minimum segment size of 1000 pixels). False positives are present in the drainage channels where landslide material and channel sediments may have propagated downstream, and in areas affected by deforestation or specific local agricultural practices.

The initial landslide training sample introduced in the ALADIM used only landslide source areas on the upper part of the hillslopes, in order to generate a shallow landslide inventory. By incorporating the sediment propagation area in the training sample, ALADIM is able to detect both the source and the deposition or runout areas. ALADIM is thus flexible to the purpose of the analysis, as the expert can decide which processes need to be detected, and can create its training sample in agreement with the study’s objectives (Stumpf et al. 2010). From a hazard perspective, detecting sediment runouts in the channels (e.g., processes like debris flows or floods) might be an important component of an emergency landslide response effort. ALADIM offers several options to spatially filter the false positives, through the introduction of landcover/landuse information (e.g., OpenStreetMap classes) or through topographic (slope thresholding) and morphological (no detection in the channels) criteria.

As for all supervised classification approaches, the accuracy and completeness of the landslide inventory is in relation with the size and quality of the input training sample. The more diverse samples are introduced in the training set, the more accurate the product is (Stumpf et al. 2013). For the rapid landslide effort in Haiti, 394 landslide samples were digitized on the Pléiades imagery and introduced in the Random Forest classifier. This digitization labor costs 1 h which we believe is an appropriate time frame taking into account that all the other components of the service are fully automated.

Although covering only a small area of the Tiburon Peninsula, the ALADIM generated landslide map has higher accuracy (compared to HazMapper) with a PA of 76.9% and a UA of 52.5%. This is not surprising as it used very-high resolution Pléiades imagery as opposed to Sentinel-2 imagery used by HazMapper. ALADIM also utilizes OBIA which has been shown to map landslides with higher accuracy than pixel-based methods like HazMapper (Moosavi et al. 2014; Keyport et al. 2018). However, ALADIM requires longer processing times than HazMapper and a more advanced computational environment set up.

After the release of the inventory, further analyses were performed by integrating a pre-event Pléiades image acquired on December 25, 2015, image not made accessible by Airbus during the rapid response period. We used ALADIM in change detection mode, by integrating in the analysis the pre-event (2015) and the post-event (2021) Pléiades images. This allowed for calculating spectral attributes in 2015 and 2021, and differences in the spectral attributes. The ALADIM change-detection mode allowed detection of landslides directly triggered by the Nippes earthquake and by Tropical Storm Grace only with fewer FPs (Fig. 15).

Fig. 15

ALADIM products generated in change detection mode (change detection processing with a pre-event image Pléiades image of December 2015 and a post-event image of August 2021): a Location of the use case in the central part of Pic Macaya National Park (e.g., overlapping area of the 2015 and 2021 Pléiades imagery); b Pre-event (December 21, 2015) Pléiades image; c Post-event (August 23, 2021) Pléiades image; d ALADIM change detection generated landslide product for the whole area; e Zoomed in view of the detected landslides; f Accuracy metrics (Pléiades ^©CNES 2021, Distribution Airbus DS)

5.5 Replicability for disaster response efforts for other landslide events

SAR heatmap, HazMapper, and ALADIM were all designed with an eye toward reproducibility. At their core, each method can be applied at any time to any place on Earth with appropriate satellite imagery. Because SAR heatmap and HazMapper are built within GEE, these algorithms require updating of a few parameters to be applied to different events. ALADIM has a more complex workflow, introducing OBIA and machine learning model, but all components are optimized in the GEP platform, are based on High-Performance Computing (HPC) and are fully automated except for the setup of the training sample.

SAR heatmap and HazMapper methods are reproducible by accessing and executing the scripts located in data availability section. A GEE account is required for access, however, though no special permissions are required. ALADIM is freely accessible on the GEP and requires only the setup of an account; the service will further be available in 2023 on the Charter Mapper application (https://eo4society.esa.int/2022/05/24/advanced-cloud-based-environment-to-support-immediate-disaster-response/) of the International Charter Space and Major Disasters. The Charter Mapper is a processing environment which integrates a suite of online EO-based services. It can access and display any open access or user-uploaded optical and SAR EO data, perform visual analysis, apply GIS functionalities, and perform EO-based processing to extract information from imagery.

We see great utility in the further development of open access methods that rely on operational workflows that do not require significant user input or the use of complicated computing systems or environments. Platforms like GEE, the Microsoft Planetary Computer and the GEP, among others, can provide the basis for these types of systems and should be leveraged whenever possible. These types of workflows also operate on the principle of easy updates when new data become available. For example, HazMapper’s URL-updating function allows analysis parameters to be embedded in a transferable URL. As new data becomes available in GEE, execution of the same original analysis URL can ingest new data and provide updated results. An additional example is the new function within ALADIM to ingest PlanetScope imagery integrated in the GEP catalog using the Spatio-Temporal Asset Catalogs protocol, allowing rapid and easy updates as new imagery comes online. The simplicity and repeatability of this style of cloud- and HPC-based satellite image analysis can greatly increase accessibility to non-specialists.

The USGS response and mapping efforts follow a typical remote sensing mapping technique of mapping landslide landforms identified from optical imagery in a GIS. A benefit of optical mapping includes the ability for an analyst to rely on expert judgement in assessing potential impacts to the built environment and other hazards, such as a landslide dam in a stream forming an upstream pond or a blocked transportation network. Because the USGS was coordinating with humanitarian groups in Haiti, this additional assessment can be beneficial to first responders. For the Haiti response, the USGS had several scientists with experience in earthquake-triggered landslides and remote sensing to contribute to the mapping and interpretation and analysis efforts, which is not always the situation. Optical imagery may continue to be useful for remote mapping in future landslide and ground failure events, especially in those events with a humanitarian response and data is needed for rapid situational awareness. Incorporating automated processes such as those described in this paper into this approach could add additional value.

5.6 Collaborative mapping efforts, scientific coordination, and future outlook

While there have been studies that compare inventories following major events such as the 2015 Gorkha Earthquake (Kargel et al. 2016; Meena and Tavakkoli Piralilou 2019), there is little published on how rapid mapping approaches compare with the goal of providing landslide information on timescales relevant for rapid disaster response (generally days to months). By evaluating these methods in terms of accuracy, completeness, data used, processing speed, and complexity, it presents new guidance to inform how rapid mapping can be undertaken following major geohazard events. The work presented herein is not exhaustive and others have also produced inventories for this event (e.g., Martinez et al. 2021a; Havenith et al. 2022; Zhao et al. 2022).

A critical lesson that has emerged from this response and the collaborative effort in writing this manuscript is the importance of coordination between different agencies before and during the immediate aftermath of the disaster. This is important to provide mutual support (e.g., data sharing, situational awareness learnings) and to ensure that efforts are aligned. For example, in the immediate aftermath of the event, authors of this manuscript were independently producing and submitting or releasing inventory products. Our collective contribution could have been more concise and potentially more accurate had we all been collaborating and sharing what we were learning about the landslide event in the early days after the event. For example, manual inventories generated by the USGS may have been able to assist in training the ALADIM automated mapping methods. Similarly, early estimations of the most concentrated landslide activity identified from SAR heatmap, ALADIM, or HazMapper could have helped USGS analysts pinpoint and conduct detailed manual mapping and confirmation.

To this end of increased coordination, the Committee on Earth Observing Satellites (CEOS) Working Group on Disasters (https://ceos.org/ourwork/workinggroups/disasters/) has a Landslide Demonstrator Group whose mission is to ensure sustained coordination of disaster-related activities undertaken by CEOS agencies. It was this group that brought the authors of this manuscript together. The goal of this working group, and other efforts, such as NASA’s Disasters Program and the implementation of the ESA Charter Mapper application, is to demonstrate the feasibility, utility, and approaches to leveraging remotely sensed data to improve situational awareness of landslide hazards. Through open-source methods to diagnose landslide behavior following major events, it provides new opportunities to empower the scientific community to collaborate and advance our understanding of landslide properties and mechanisms while also providing key data more rapidly to support disaster response.

End users and first responders participating in the relief efforts likely lack the time, focus, and nuanced understanding of each product to effectively gain insights from the data. Submitting the same type of data (e.g., landslide inventories) from multiple methods further confounds the matter. One of the largest challenge associated with automated methods is the presence of FPs. ALADIM and HazMapper did not employ any FP removal step. For optical products, filtering steps like slope threshold and use of land cover and susceptibility information can be used to minimize FPs. However, the landslide's spectral characteristic matches with many other objects in the satellite imagery such as barren areas, agricultural areas, and riverbanks. Hence, a method to remove all FPs without also removing from TPs might truly not exist. As automated methods are able to quickly produce results, it may be worthwhile to also employ manual correction so that the responders have the best data possible to make decisions. This can occur in an iterative loop in the days to months following the disaster. For SAR products this step is very challenging even with human intervention. The SAR heatmap method used here employed a slope threshold to reduce FPs. When available, post-event multitemporal image stacking further reduces FPs by suppressing noise and therefore makes landslide areas more prominent (Handwerger et al. 2022). SAR imagery suffers from geometric distortion (foreshortening, shadow, and layover) and some areas within the imagery does not contain any usable data. Identifying areas which cannot be seen by SAR and masking the final heatmap product will be a very useful step to take to avoid confusion. In the same order, sharing of cloud and shadow mask for optical-based products will be beneficial as well.

Ultimate usability and suitability of the products generated to support rapid response comes from feedback provided by the users. Landslide products were disseminated to different end users by each of the involved groups. Following any disaster response, an after-action review (AAR) should be performed as a process to identify aspects of a response that worked well and what needs improvement. For this event, the outcome of the AAR was limited and challenging to interpret for future planning purposes. This is in part because the bulk of the landslides occurred in Pic Macaya National Park, where there were fewer direct impacts to population and infrastructure. This means the interests of stakeholders was geographically separated from the bulk of the mapping effort. In addition, the novelty of some products means users were unfamiliar with the tools, rendering them unable to gauge relevancy for use in a disaster scenario. Future response efforts should make further efforts to connect data providers with first responders to ensure proper guidance on how best to use the data is communicated from the beginning and an AAR should be conducted after the response has ended with participation of all parties involved.

6 Conclusion

This study was undertaken by seven separate government, academic, and industry partners with the goal of centralizing our knowledge of how different rapid response efforts to remotely detect landslides compare to on-the-ground conditions. We used the August 2021 Nippes earthquake and Tropical Storm Grace as a case study. In the aftermath of these events, both manual and automated methods were used for mapping landslides utilizing SAR and optical imagery. The SAR heatmap can only highlight areas with dense landslide activity and might not be suitable in areas with few or singular landslides. The rdNDVI-based HazMapper is fast but highlights vegetation change as possible landslides, leading to high number of FPs. The GEE-based SAR heatmap and HazMapper methods are easily accessible and produce results relatively quickly depending upon imagery availability. These methods can be very effective at establishing first estimates of landsliding. The OBIA and machine learning-based ALADIM has reasonable accuracy but takes more time and effort as it requires generation of training samples and object-based processing. These open-source algorithms have the potential to be widely utilized by the scientific community for other geohazards to better inform rapid landslide mapping. The main issue with these automated methods is FPs. First responders might not have the time or ability to understand limitations associated with the results. As these methods can generate landslide data quickly, it may be worthwhile for data producers to remove FPs by adopting a post processing filtering step or even manually. Manual mapping resulted in less FPs but was labor- and time-intensive. Future efforts should attempt to balance these advantages and disadvantages of the various methods.

We also acknowledge that there are many other automated methods which can be used other than the ones used in this study. However, testing of these methods is needed in a rapid response scenario where quick delivery is a must as opposed to the development phase where one has the luxury of testing various configurations for optimal mapping. Effective channels of collaboration and communication should be in place such that different groups are collaborating and not replicating efforts already underway.

The real utility of these rapid response landslide products can only be quantified using user feedback. An AAR should be initiated by the response programs for future events connecting data providers to the stakeholders to help improve response, coordination, and communication.