Using Google Earth Engine to examine the Haines landslides of 2020
Erin Dawn Trochim, Edjoyce, Iraluq Bond, Noelle Helder
Part 1
On December 2, 2020, the residents of Haines, Alaska were struck by a devastating landslide that tore through Beach Road, destroying homes and claiming the lives of two individuals. The disaster left a huge scar on the side of Mt. Riley, and it’s hard to imagine the force that could have caused such destruction.
The Beach Road landslide was just one of many that occurred in the traditional territory of the Tlingit, where Haines is known as Deishú in the Lingít language. These landslides were not isolated incidents, but rather a widespread phenomenon in the region.
Examining the Southeast Alaska Geology below is a good starting place. There is also existing landslide inventory, where it is important to note the area to the east of Haines was not surveyed. Over the main townsite of Haines, landslides occurred in modern surficial deposits, composed of sediment recently transported and deposited in channels and washes, on alluvial fans and plains, and on hillslopes. The deposits (Qs) are undifferentiated alluvial deposits of the late Holocene, made up of sand and pebbles to small cobble gravel that have not been assigned to any specific surficial unit. The Beach Road Landslide occurred over ultramafic rocks (Kum), and smaller landslides were also found in the gabbro (Kdb) unit.
Beyond the geology, there is a substantial amount of earth observation data over the region. In this series of blogs, we will use the power of Google Earth Engine to examine the extent of the landslides in the region and understand the spatial impact of the landslides on the community. We will see how the landslides, driven by unusually heavy rain, affected the area and the residents.
A significant increase in Arctic precipitation since the 1950s is now detectable across all seasons due to climate change. Trend analysis indicates increased precipitation in the Haines region. Landslides are likely to become more frequent as a result of this increased precipitation. The trend towards higher levels of precipitation in the Haines region is a clear indication of the impact of climate change on the region’s geohazards, making it imperative for local communities to adapt and prepare for the increased risk of landslides.
While this analysis will insight into the earth observation data available to examine this event, it is important to recognize the impact this event had on residents of the area. Over 50 households had to evacuate, and relief efforts were complicated by the ongoing COVID-19 pandemic. Recovery is expected to take up to a decade. Keeping this in mind, addressing individual and social ethics are a critical component of modern remote sensing. The goal of this effort is to increase data accessibility and better understand data now available to study these types of events, while remembering that recovery is still ongoing.
A team of experts analyzed the Beach Road Landslide and published their findings in the journal Landslides (Darrow et al. 2022). The study revealed that the landslide was approximately 690 meters long, 170 meters wide, with a volume of around 187,000 cubic meters, and an average slope angle of 20 degrees. The cause of the landslide was determined to be an “atmospheric river” that generated rain-on-snow, which is a dangerous combination for landslide initiation. The failure began as a debris avalanche and then transitioned into a debris flow. One of the most intriguing aspects of the study is the history of the ridge from which the failure originated indicating recent paleo-landslides. This points to a persistent failure risk.
Join us as we explore the findings of this important study, gain a deeper understanding of the Haines landslide and use Earth Engine to see the extent of landslides in the region.
First, we’ll start by examining the magnitude and distribution of precipitation during the period around the event using the ERA5-Land reanalysis dataset.
The ERA5-Land dataset has a spatial resolution of about 11 km with hourly data going back to 1981 of 50 different climate variables.
If you need some instructions getting started in Earth Engine, check out the first few chapters in the new textbook Cloud-Based Remote Sensing with Google Earth Engine.
Start by copying in the study area polygon over Haines and loading the ERA5-Land dataset into GEE. Then we will filter the ERA5-Land image collection spatially to the study area and temporally to the week before and after the Haines landslides on December 2, 2020.
// Load the study area over Haines
var studyArea = ee.Geometry.Polygon(
[[[-135.78061942511923, 59.404719431545125],
[-135.78061942511923, 58.951587213737405],
[-134.93467215949423, 58.951587213737405],
[-134.93467215949423, 59.404719431545125]]], null, false);
// Visualize the study area on a map
Map.addLayer(studyArea)
Map.centerObject(studyArea, 9);
// Load the ERA5-Land dataset
var era5Land = ee.ImageCollection("ECMWF/ERA5_LAND/HOURLY")
.filterBounds(studyArea)
// December 2, 2020 is DOY 336
// Filter dates for one week before and after event
.filter(ee.Filter.dayOfYear(336-7, 336+7));Check out the band descriptions on the ERA5-Land reanalysis dataset. There are several options that describe precipitation. All of them are in meters.
In this exercise, we’ll use hourly snowfall and total precipitation, including both rain and snow, produced by the Earth Engine team from accumulated daily versions resetting at midnight. We’ll write a function to convert the precipitation data from meters to centimeters. After filtering for 2020 data, we’ll apply the conversion function to by mapping it to the image collection.
// Convert from meters to centimeters
var convertCM = function(image){
// Multiply the image by 100 to convert from m to cm
var imageCM = image.multiply(100);
// Return the image with the correct date information
return imageCM.copyProperties(image, ['system:index', 'system:time_start']);
}
// Select the hourly snowfall and total precipitation data
var precip2020 = era5Land.select(['snowfall_hourly', 'total_precipitation_hourly'])
// Filter to the year 2020
.filter(ee.Filter.calendarRange(2020, 2020, 'year'))
// Convert to cm
.map(convertCM)Now let’s look at the results. We’ll start by creating an interactive chart comparing snowfall and total precipitation. Using the code below, select the precip2020 image collection, limit the analysis to the study area, apply a mean reducer (ee.Reducer.mean()) using the scale set to the spatial resolution of the ERA-5 Land dataset for each image. The xProperty is the date property of the hourly images.
// Define the chart and print it to the console
var chart =
ui.Chart.image
.series({
// Select the correct image collection
imageCollection: precip2020,
// Limit the analysis to the study area
region: studyArea,
// Reduce the images to the mean of the study area
reducer: ee.Reducer.mean(),
// Set the scale as the spatial resolution of the ERA5-Land dataset
scale: 11132,
// Display the data by hour using the date property
// of the dataset (system:time_start)
xProperty: 'system:time_start'
})
// Set the names of the two different bands
.setSeriesNames(['snowfall', 'total_precipitation'])
.setOptions({
title: 'Hourly Preciptation by Date',
hAxis: {title: 'Date', titleTextStyle: {italic: false, bold: true}},
vAxis: {
title: 'Precipitation type (cm)',
titleTextStyle: {italic: false, bold: true}
},
// Alter these settings to control the line properties
// 0 is snowfall and 1 is total preciptiation
series: {
0: {lineWidth: 3, color: '57e1de'},
1: {lineWidth: 3, color: '1d6b99'}
},
curveType: 'function'
});
print(chart);You will see a chart like the one below in your console. If you select the pop out icon on the upper right of the chart, then you have the additional options to download the data as a csv, or the chart as svg or png. Check additional chart styling options to further customize the look.
Using your cursor to check the dates, notice that late in the day Dec. 1, 2020 the snowfall values diverge from total precipitation. This indicates the ERA-5 Land model producing a mix of rain and snow. This strong trend continued until the early morning of Dec. 3.
We need to remember these are the mean values over our study area for each hour. A reducer aggregates over properties, temporal and/or spatial dimensions. In our case, we are computing the average total snowfall and total precipitation for the study area using a version of the reduceRegion() function embedded in the series function.
What happens if you change the reducer in the chart function to ee.Reducer.max()? Do you notice the values slightly increasing? What does this indicate?
Let’s check by creating a map of the precipitation data. We’ll transform the data from hourly to daily in order to simplify the interpretation.
We will do this by creating a new date variable called day of year (DOY).
The code below creates a function called addDOY, which takes an image as an input and calculates the day-of-year (DOY) property using the image’s date. The precip2020 image collection is then passed through the function to create the new precip2020DOY image collection with the DOY property added. Next, a list of distinct days is created from precip2020DOY using the aggregate_array and distinct methods. Then we define a new function called map_to_sequence_days which filters the precip2020DOY collection to each day, sums all values and sets the DOY. Lastly, a new image collection precip2020Daily is created which contains the daily precipitation sum for each day.
// Function to add day-of-year property to each image
var addDOY = function(image){
var DOY = ee.Date(image.date()).getRelative('day', 'year')
return image.set('doy', DOY);
}
// Apply DOY function to precip2020 image collection
var precip2020DOY = precip2020.map(addDOY)
// Create list of days
var days = precip2020DOY.aggregate_array('doy').distinct();
// Print list to check days
print(days)
// Function to filter to each DOY, sum all values and set DOY
var map_to_sequence_days = function(day){
return precip2020DOY.filter(ee.Filter.eq('doy', day)).reduce(ee.Reducer.sum()).set('doy', day);
}
// Sum daily precipitation
var precip2020Daily = ee.ImageCollection.fromImages(
days.map(map_to_sequence_days)
);Now we are ready to create visualizations of snowfall and total precipitation data. The code below sets a color palette and parameters for the visualizations. Next it adds two layers to the map, one for snowfall and the other for total precipitation. It also creates a color bar and legend to show the scale of the data in the visualizations.
// Set palette for visualizations
var palette = ['#ffffcc','#a1dab4','#41b6c4','#2c7fb8','#253494']
// Set visualization parameters
var vis = {min: 0, max: 10, palette: palette};
// Set number of steps between max and min for legend
var nSteps = 10
// Add maps to console to interactively explore timeseries
Map.addLayer(precip2020Daily.select('snowfall_hourly_sum'),
vis, 'Snowfall')
Map.addLayer(precip2020Daily.select('total_precipitation_hourly_sum'),
vis, 'Total precipitation')
// Creates a color bar thumbnail image for use in legend from the given color palette
function makeColorBarParams(palette) {
return {
bbox: [0, 0, nSteps, 0.1],
dimensions: '100x10',
format: 'png',
min: 0,
max: nSteps,
palette: palette,
};
}
// Create the colour bar for the legend
var colorBar = ui.Thumbnail({
image: ee.Image.pixelLonLat().select(0).int(),
params: makeColorBarParams(vis.palette),
style: {stretch: 'horizontal', margin: '0px 8px', maxHeight: '24px'},
});
// Create a panel with three numbers for the legend
var legendLabels = ui.Panel({
widgets: [
ui.Label(vis.min, {margin: '4px 8px'}),
ui.Label(
((vis.max-vis.min) / 2+vis.min),
{margin: '4px 8px', textAlign: 'center', stretch: 'horizontal'}),
ui.Label(vis.max, {margin: '4px 8px'})
],
layout: ui.Panel.Layout.flow('horizontal')
});
// Legend title
var legendTitle = ui.Label({
value: 'Snowfall / Total Precipitation (cm)',
style: {fontWeight: 'bold'}
});
// Add the legendPanel to the map
var legendPanel = ui.Panel([legendTitle, colorBar, legendLabels]);
Map.add(legendPanel);You can adjust the transparency of the different image collections using the sliders beside their names in Layers tab to make your visualization look like the one below (left). Explore it further by selecting the Inspector tab and drop down the Series menus under the different image collections. These show the individual pixel values for each day at the selected point.
Earth Engine provides better ways to visualize the data. One of these produces a filmstrip of each individual day.
Use the following code to define the parameters for the getFilmstripThumbURL function, which is used to generate a filmstrip of the data. The arguments include the dimensions of the images in the filmstrip, the region of the data being studied, the coordinate reference system, the minimum and maximum values, and the color palette. The code then prints a URL for each of the two data layers, snowfall_hourly_sum and total_precipitation_hourly_sum. Accessing the URL will generate a filmstrip of images that can be downloaded.
// Define arguments for the getFilmstripThumbURL function parameters
var filmArgs = {
dimensions: 128,
region: studyArea,
crs: 'EPSG:3338',
min: 0,
max: 10,
palette: palette
};
// Print a URL that will produce the filmstrip when accessed
print(precip2020Daily.select('snowfall_hourly_sum').getFilmstripThumbURL(filmArgs));
// Print a URL that will produce the filmstrip when accessed
print(precip2020Daily.select('total_precipitation_hourly_sum').getFilmstripThumbURL(filmArgs));If you download the filmstrip and arrange it manually, the daily snowfall values look like the figure below. I created this version in Google Drawing. You can copy the legend from the console. If you want a more automated process that automatically annotates each image, you can replicate the code in the python API using a similar approach to this example which chops up the filmstrip into an animated gif.
Let’s examine the cutout section in the data over the ocean. The name of the dataset gives a hint to what is happening — ERA5-Land — it indicates the original ERA5 dataset which is global has been resampled for land areas. This improved the spatial resolution from ~22 km to ~11 km in the land model. However, the data doesn’t go over ocean areas.
We will create a point feature the specific longitude and latitude of the Beach Road Landslide. You can do this using the ee.Geometry.Point() function, which takes a list of coordinates as an argument. Then we add the point feature to the map as a red dot on the map.
// Create a point at the Beach Road landslide
var beachRoad = ee.Geometry.Point([-135.412, 59.219])
// Add it to the map
Map.addLayer(beachRoad, {color: 'FF0000'}, 'Beach Road Landslide')
Unfortunately the Beach Road Landslide is within the ocean mask in the ERA5-Land dataset. Does this mean this data is not useful — no. You could interpolate the raster data, or use the nearest valid pixel for features like the Beach Road Landslide.
The other option is to look at other datasets. Earth Engine makes is easy to trade, examine and compare multiple datasets. A new option for Alaska is the Daymet V4: Daily Surface Weather and Climatological Summaries dataset. We need to switch our code in a few places: A) the input dataset; B) the band selection and C) the spatial resolution of the reducer. This dataset has two bands describing precipitation. We are only going to look at prcp, which is daily total precipitation in millimeters which sums all forms of water to an equivalent value.
The Daymet data has a spatial resolution of 1 km and as a daily product, it has a temporal resolution of 1 day.
The following code demonstrates how to adapt the previous code to the Daymet data and looks at only the Beach Road Landslide site.
// Load the Daymet dataset
var daymet = ee.ImageCollection("NASA/ORNL/DAYMET_V4")
.filterBounds(studyArea)
// December 2, 2020 is DOY 336
// Filter dates for one week before and after event
.filter(ee.Filter.dayOfYear(336-7, 336+7))
// Select only prcp and swe
.select(['prcp']);
// Convert from millimeters to centimeters
var convertMM = function(image){
// Divide the image by 10 to convert from mm to cm
// Return the image with the correct image properties
return image.divide(10).copyProperties(image, ['system:index', 'system:time_start']);
}
// Select the hourly snowfall and total precipitation data
var precipDay2020 = daymet.filter(ee.Filter.calendarRange(2020, 2020, 'year'))
// Convert to cm
.map(convertMM)
// Define the chart and print it to the console
var chart =
ui.Chart.image
.series({
// Select the correct image collection
imageCollection: precipDay2020,
// Limit the analysis to the Beach Road Landslide
region: beachRoad,
// Reduce the images to the max of the Beach Road Landslide
reducer: ee.Reducer.max(),
// Set the scale as the spatial resolution of the Daymet dataset
scale: 1000,
// Display the data by hour using the date property
// of the dataset (system:time_start)
xProperty: 'system:time_start'
})
// Set the names of the band
.setSeriesNames(['prcp'])
.setOptions({
title: 'Daily Daymet Preciptation at Beach Road Landslide by Date',
hAxis: {title: 'Date', titleTextStyle: {italic: false, bold: true}},
vAxis: {
title: 'Precipitation (cm)',
titleTextStyle: {italic: false, bold: true}
},
// Alter these settings to control the line properties
// 0 is total preciptiation
series: {
0: {lineWidth: 3, color: '1d6b99'}
},
curveType: 'function'
});
print(chart);Comparing the overall results ERA5-Land results to the Daymet data below demonstrates differences. At the Beach Road Landslide, the Daymet data indicates maximum precipitation on Dec. 1, 2020. The ERA5-Land data showed high amounts of precipitation on both Dec. 1 and Dec. 2, 2020. The ERA5-Land data is a reanalysis dataset, while the Daymet product is a gridded weather dataset.
The easiest way to directly compare the datasets is to filter each dataset for an individual day and then subtract the Daymet dataset from the ERA5-Land. Find the code below to complete this task and compare the output images.
// Create Daymet Dec 2 data
var daymetDec2 = precipDay2020.filter(ee.Filter.dayOfYear(336, 336))
.mean()
// Create ERA5-land Dec 2 data
var eraDec2 = precip2020Daily.filter(ee.Filter.equals('doy', 336))
.mean()
.select('total_precipitation_hourly_sum')
// Difference ERA5-Land from Daymet data for Dec 2
var diffDec2 = eraDec2.subtract(daymetDec2)
// Create Daymet Dec 1 data
var daymetDec1 = precipDay2020.filter(ee.Filter.dayOfYear(335, 335))
.mean()
// Create ERA5-land Dec 1 data
var eraDec1 = precip2020Daily.filter(ee.Filter.equals('doy', 335))
.mean()
.select('total_precipitation_hourly_sum')
// Difference ERA5-Land from Daymet data for Dec 1
var diffDec1 = eraDec1.subtract(daymetDec1)
// Set palette for visualizations
var palette = ['blue','grey','red']
// Set visualization parameters
var vis = {min: -5, max: 5, palette: palette};
// Set number of steps between max and min for legend
var nSteps = 10
// Add maps to console to interactively explore timeseries
Map.addLayer(diffDec1,
vis, 'Difference Dec. 1, 2020')
Map.addLayer(diffDec2,
vis, 'Difference Dec. 2, 2020')
// Creates a color bar thumbnail image for use in legend from the given color palette
function makeColorBarParams(palette) {
return {
bbox: [0, 0, nSteps, 0.1],
dimensions: '100x10',
format: 'png',
min: 0,
max: nSteps,
palette: palette,
};
}
// Create the colour bar for the legend
var colorBar = ui.Thumbnail({
image: ee.Image.pixelLonLat().select(0).int(),
params: makeColorBarParams(vis.palette),
style: {stretch: 'horizontal', margin: '0px 8px', maxHeight: '24px'},
});
// Create a panel with three numbers for the legend
var legendLabels = ui.Panel({
widgets: [
ui.Label(vis.min, {margin: '4px 8px'}),
ui.Label(
((vis.max-vis.min) / 2+vis.min),
{margin: '4px 8px', textAlign: 'center', stretch: 'horizontal'}),
ui.Label(vis.max, {margin: '4px 8px'})
],
layout: ui.Panel.Layout.flow('horizontal')
});
// Legend title
var legendTitle = ui.Label({
value: 'Difference (cm)',
style: {fontWeight: 'bold'}
});
// Add the legendPanel to the map
var legendPanel = ui.Panel([legendTitle, colorBar, legendLabels]);
Map.add(legendPanel);There is an obvious difference between the ERA5-Land and Daymet data between the timing of the maximum precipitation. Daymet shows a clear trend towards maximum precipitation on Dec. 1, 2020 while the ERA-5 Land model indicates it occurs the next day. Looking at a map of this analysis indicates that the pattern is common to both Haines and the area surrounding the Beach Road Landslide.
The choice of dataset to use in further analysis depends on a variety of factors. The Daymet data directly covers the Beach Road Landslide. An alternative ERA5-Land data is the ERA5 Daily Aggregates dataset, which is available globally over both land and ocean at approximately 22 km spatial resolution. We can validate the timing by comparing our results against the NOAA weather station at the Haines Airport as shown below. The double peak of precipitation on both Dec. 2 and 3 can be seen in the measured data, so the ERA5-Land data is more accurately capturing the daily variability.
In terms of comparing absolute magnitude, you could rerun the ERA-5 Land code to look at only the location of the weather station. Its location is Lat/Lon: 59.24290/-135.51140. Right now it is important to remember we were looking at a slightly larger region summing multiple pixels. Overall though, the magnitude looks similar.
Examining the climate data this way is a key part of the research process, it is exploratory data analysis. In our modern, data-rich era there are often multiple datasets available. Taking the time to understand the strengths and variability of each dataset is important for informing further analysis. We can also use the spatial characteristics of the climate data to examine landscape dynamics further within packages like Landlab in python.
Earth Engine makes it easy to examine a much larger area in SE Alaska. All this requires is a larger study area. The other type of analysis that is easily enabled is comparing daily precipitation values to historical norms. This simply involves creating a mean reducer where the years are set appropriate.
Follow along the next steps in Part 2 in this series where we explore the distribution of landslides further in the Haines area using Dynamic World data.
