Using Google Earth Engine to examine the Haines landslides of 2020
Erin Dawn Trochim, Edjoyce, Iraluq Bond, Noelle Helder
Part 2
In this exercise, we will continue to examine the Beach Road Landslide and surrounding landslides in Haines, Alaska. In Part 1 of this blog series, we used Google Earth Engine to examine precipitation patterns surrounding this event using ERA5-Land and Daymet climate data available in the Earth Engine Data Catalog. We conducted exploratory data analysis to compare spatial and temporal variability in the datasets.
Here, we will use a powerful dataset known as Dynamic World. It is a 10m near-real-time (NRT) Land Use/Land Cover (LULC) dataset that includes class probabilities and label information for nine different classes. It offers a wealth of information on land use and cover, which can be incredibly useful for understanding the impact of landslides on an area.
Dynamic World predictions are available for the Sentinel-2 L1C collection from 2015–06–27 to present. The revisit frequency of Sentinel-2 is between 2–5 days depending on the latitude. Dynamic World predictions are generated for Sentinel-2 L1C images with a cloudy pixel percentage of less than or equal to 35%. The predictions are also masked to remove clouds and cloud shadows using a combination of S2 Cloud Probability, Cloud Displacement Index, and Directional Distance Transform. All probability bands except the “label” band collectively sum to 1.
If you’re new to the Dynamic World dataset, an introduction can be found here. Parts of this exercise are adapted from the part 3 of the series.
The nine different Dynamic World classes are:
- Water: Permanent and seasonal water bodies
- Trees: Includes primary and secondary forests, as well as large-scale plantations
- Grass: Natural grasslands, livestock pastures, and parks
- Flooded vegetation: Mangroves and other inundated ecosystems
- Crops: Include row crops and paddy crops
- Shrub & scrub: Sparse to dense open vegetation consisting of shrubs
- Built area: Low- and high-density buildings, roads, and urban open space
- Bare ground: Deserts and exposed rock
- Snow & ice: Permanent and seasonal snow cover
For looking at landslides, an obvious place to start is with the Bare ground class as it includes deserts and exposed rock which is often associated with landslide-prone areas. Also, the Built area class could be useful as it includes low- and high-density buildings, roads and urban open space which can provide information on the impact of landslides on the infrastructure.
We will start by creating a simple map and chart to visualize data from the Dynamic World dataset. The code below creates a point at a specific location, the Beach Road landslide, and adds it to the map. The data is then filtered by a specific date range ~ a month before and after the Dec. 2, 2020 landslides. It is also filtered by the location of the Beach Road Landslide. Then we create a chart to display the probabilities of different classes (e.g. water, trees, grass) at the landslide over time. The chart is visualized by printing it to the console.
// 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')
Map.centerObject(beachRoad, 16);
// Set start and end dates
var startDate = '2020-11-01';
var endDate = '2021-02-01';
// Load the Dynamic World dataset
// Filter for start and end dates and the Beach Road Landslide
var dw = ee.ImageCollection('GOOGLE/DYNAMICWORLD/V1')
.filterDate(startDate, endDate)
.filterBounds(beachRoad);
// Create a list to select all probability bands
var probabilityBands = [
'water', 'trees', 'grass', 'flooded_vegetation', 'crops', 'shrub_and_scrub',
'built', 'bare', 'snow_and_ice'
];
// Select these bands
var dwTimeSeries = dw.select(probabilityBands);
// Create a function to make a style dictionary for each probability band
var lineStyle = function(label, color) {
var styleDict =
{labelInLegend: label, color: color, lineWidth: 2, pointSize: 3};
return styleDict;
};
// Define the chart and print it to the console.
var chart =
ui.Chart.image
.series({
// Select the correct image collection
imageCollection: dwTimeSeries,
// Limit the analysis to the Beach Road Landslide
region: beachRoad,
// Set the scale as the spatial resolution of the Dynamic World dataset
scale: 10,
// 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: 'Dynamic World data at Beach Road Landslide by Date',
hAxis: {title: 'Date', titleTextStyle: {italic: false, bold: true}},
vAxis: {
title: 'Class probabilities',
titleTextStyle: {italic: false, bold: true},
viewWindow: {min: 0, max: 1}
},
interpolateNulls: true,
series: {
0: lineStyle('Bare', '#A59B8F'),
1: lineStyle('Built', '#C4281B'),
2: lineStyle('Crops', '#E49635'),
3: lineStyle('Flooded vegetation', '#7A87C6'),
4: lineStyle('Grass', '#88B053'),
5: lineStyle('Shrub and scrub', '#DFC35A'),
6: lineStyle('Snow and ice', '#B39FE1'),
7: lineStyle('Trees', '#397D49'),
8: lineStyle('Water', '#419BDF')}
});
print(chart);The Beach Road Landslide is not obviously detectable in this chart. One challenge is the lack of data. The last image before the landslide was collected Nov. 17, 2020 while the next successful collection was over two months later on January 21, 2021.
The issue with satellites like Sentinel-2 is they are unable to penetrate clouds. It appears the period before and after the landslide was cloudy enough to prevent data collection.
We can use the pop out icon on the upper right of the chart to interactively examine the properties on the chart. Note the class probability for trees is 0.092 before the landslide and 0.046 after. In comparison, bare ground went from 0.057 to 0.059. Snow and ice increased while water decreased.
We will try a different approach. In Earth Engine, we’ll use our larger study area around Haines and filter two collections of data, one from Dynamic World and another from Sentinel-2. The code creates two filters, one called colBEFORE and another called colAFTER. The colBEFORE filter is defined by the study area and a date range of June 1st, 2020 to September 30th, 2020. The colAFTER filter is defined by the same study area and a date range of June 1st, 2021 to September 30th, 2021. These filters are used to filter the Dynamic World and Sentinel-2 data, so that only data from the specified date range and within the study area is used in further analysis.
// Examine a larger study area around Haines
var studyArea = ee.Geometry.Polygon(
[[[-135.78061942511923, 59.404719431545125],
[-135.78061942511923, 58.951587213737405],
[-134.93467215949423, 58.951587213737405],
[-134.93467215949423, 59.404719431545125]]], null, false);
// Construct a collection of corresponding Dynamic World and Sentinel-2 for
// inspection. Filter the DW by region and date
var colBEFORE = ee.Filter.and(
ee.Filter.bounds(studyArea),
ee.Filter.date('2020-06-01', '2020-09-30'));
var colAFTER = ee.Filter.and(
ee.Filter.bounds(studyArea),
ee.Filter.date('2021-06-01', '2021-09-30'));Next we will create two median composites of the bare class from the Dynamic World dataset, one for before and one for after the landslide. We use the colBEFORE and colAFTER filters defined in the previous code. Since we will want to use this process several time, this code constructs a function to select the desired bare class from the dataset and calculates the mosiac value of the data using the specified ee.Reducer. We create a binary image using the greater than operator where pixels that have a value greater than the given value are set to 1 and the rest to 0. This binary image is then added to the map twice, once for the before landslide data and once for the after landslide data, each with different color schemes to distinguish between the two. Then we calculate a difference image between the before and after images.
This code applied the function to the bare class using a median reducer and a threshold of 0.3
// Function to create composite of class before and after the landslides
var compositeBeforeAfter = function (pClass, mosiac, threshold){
// Before landslide
var dwBefore = ee.ImageCollection('GOOGLE/DYNAMICWORLD/V1')
.filter(colBEFORE)
// Only look at tree class
.select(pClass)
.reduce(mosiac)
.gte(threshold);
// After landslide
var dwAfter = ee.ImageCollection('GOOGLE/DYNAMICWORLD/V1')
.filter(colAFTER)
.select(pClass)
.reduce(mosiac)
.gte(threshold);
return ee.Image.cat([dwBefore, dwAfter]).rename(['before', 'after'])
}
// Apply the function to the bare class probabilities
var dwBare = compositeBeforeAfter('bare', ee.Reducer.median(), 0.3)
Map.addLayer(dwBare.select('before').updateMask(dwBare.select('before').gt(0)), {min: 0, max: 1, palette:['black', '#d2cdc7']}, "Before landslides bare probability")
Map.addLayer(dwBare.select('after').updateMask(dwBare.select('after').gt(0)), {min: 0, max: 1, palette:['black', '#A59B8e']}, "After landslides bare probability")
// Create a difference image of before and after bare class probabilties
var dwDiffBare = dwBare.select('before').subtract(dwBare.select('after'))
Map.addLayer(dwDiffBare.updateMask(dwDiffBare.gt(0)), {min: 0, max: 1, palette:['black', 'yellow']}, "Difference bare probability")
// Make the background map dark
Map.setOptions('Dark',{
'dark': [
{ 'elementType': 'labels', 'stylers': [ { 'visibility': 'off' } ] },
{ 'elementType': 'geometry', 'stylers': [ { 'color': '#808080' } ] },
{ 'featureType': 'water', 'stylers': [ { 'color': '#404040' } ] }
]
})This analysis produces the following figure where the medium grey outline north of the red dot and Beach Road Landslide shows bare ground that is present in the summer 2021, but not during the same period in 2020. The difference image is empty because the bare ground is not in the before image.
This approach only shows the largest landslide, one on Beach Road. The other smaller landslides found around Haines do not appear using this approach. The main other features in the bare class are found by the airport on the west side of the town along the Haines Highway. In this area, there is change in bare ground class probabilities in the gravel pit and out in the river channel.
Let’s try this approach again, but instead we will examine the tree class probabilities. Adapt your code as shown below to only select the trees class. This time we will create a mosaic using maximum values thresholded so they are greater or equal to 0.5. We will add these layers to the map before calculating the difference between the before and after images to better highlight the changes.
// Apply the function to the tree class probabilities
var dwTree = compositeBeforeAfter('trees', ee.Reducer.max(), 0.5)
Map.addLayer(dwTree.select('before').updateMask(dwTree.select('before').gt(0)), {min: 0, max: 1, palette:['black', '#9cbea4']}, "Before landslides tree probability")
Map.addLayer(dwTree.select('after').updateMask(dwTree.select('after').gt(0)), {min: 0, max: 1, palette:['black', '#397D49']}, "After landslides tree probability")
// Create a difference image of before and after tree class probabilties
var dwDiffTree = dwTree.select('before').subtract(dwTree.select('after'))
Map.addLayer(dwDiffTree.updateMask(dwDiffTree.gt(0)), {min: 0, max: 1, palette:['black', 'red']}, "Difference tree probability")As shown below, the difference between the maximum summer/fall tree probabilities in 2020 and 2021 which were greater than 0.5 (as shown in red below) identify a number of probable landslides in the area. This includes a larger outline for the Beach Road Landslide.
We will check the change built class before moving on. Use the code as adapted below to explore.
// Apply the function to the built class probabilities
var dwBuilt = compositeBeforeAfter('built', ee.Reducer.median(), 0.3)
Map.addLayer(dwBuilt.select('before').updateMask(dwBuilt.select('before').gt(0)), {min: 0, max: 1, palette:['black', '#9cbea4']}, "Before landslides built probability")
Map.addLayer(dwBuilt.select('after').updateMask(dwBuilt.select('after').gt(0)), {min: 0, max: 1, palette:['black', '#397D49']}, "After landslides built probability")
// Create a difference image of before and after built class probabilties
var dwDiffBuilt = dwBuilt.select('before').subtract(dwBuilt.select('after'))
Map.addLayer(dwDiffBuilt.updateMask(dwDiffTree.gt(0)), {min: 0, max: 1, palette:['black', 'orange']}, "Difference built probability")Comparing the built class before and after the landslide gives some information, but it is more challenging to consistently interpret. For example, if you look at the landslides along the highway heading north out of Haines, they appear in the tree class difference but not in the built one. This is because the road was quickly rebuilt after it was damaged by the adjacent landslides, so the class probability did not significantly vary between the summer of 2020 and 2021.
Similar to our first exercise examining precipitation, carefully examining and matching the class probabilities to the actual phenomena is critical. We are adding more exploratory data analysis to our knowledge base.
Another way to validate what we are seeing is to examine the visible bands of the Sentinel-2 dataset. We could add these images to our current code editor, this is covered in here in the first Dynamic World tutorial. An alternative is to quickly check using another prebuilt Earth Engine app. Many of them exist, check out the current list compiled by Awesome Earth Engine Apps. We will use the S2-Cloudless app to compare visible images from before and after. Start by putting “Haines, Alaska, USA” in the search bar of the app and then place a point over Haines. Use the target data selector to look at images before and after the landslides. Due to the snow cover and clouds, it is difficult to detect the landslides until there is a clear image January 21, 2021. Can you see them visible as snow covered areas where there are trees missing?
As our final task for this exercise, we will do a comparison of landslide distribution to topography. Based on our initial work, it appears the before/after tree class probability gives the most consistent results. We will use another resource, the Awesome Earth Engine Community Catalog to find an appropriate topography product for this analysis. The Community Catalog contains peer-reviewed research datasets ready to use just like those in the main Earth Engine Data Catalog. If you look under the Elevation and Bathymetry category on the left-hand side of the website.
We will use the FABDEM (Forest And Buildings removed Copernicus 30m DEM) produced using imagery between 2014 and 2018. All details about the dataset can be found in the description in the Community Catalog including an examples of how to use the data and the correct peer-reviewed citations.
Going from the initial FABDEM to the zonal analysis is a three step process.
Start by reviewing the zonal analysis chapter in the Cloud-Based Remote Sensing with Google Earth Engine textbook. Notice in Section 3. Neighborhood Statistic Examples that this approach is very similar to our desired output. They also show how calculate slope from the DEM. This example is extracting slope from different points. Currently our tree class difference representing landslides is a binary image. For this example we are going to take a more of a grouped reduction approach to zonal analysis through image to image comparisons. Since we only have one class, we will make the analysis simpler by masking out areas that do not have probably landslides and analyzing the results with a histogram. However, you could do a raster/vector conversion, make our landslide image into a polygon feature collection representing individual landslides, and follow along with the approach in the textbook.
For our three step grouped reduction approach, start by adding the FABDEM to the existing code we’ve been using, so you can reuse the tree class difference image. Then calculate slope. Use the tree class difference image to updateMask (or mask out) areas of slope that do not have landslides. Then create the histogram chart.
// Import the FABDEM dataset
var fabdem = ee.ImageCollection("projects/sat-io/open-datasets/FABDEM");
// Explanation on setting default Projection here https://twitter.com/jstnbraaten/status/1494038930643042309
var elev = fabdem.mosaic().setDefaultProjection('EPSG:3857',null,30)
// Calculate slope from the DEM
var slope = ee.Terrain.slope(elev);
// Add the slope to the difference image of before and after tree class probabilties
var slopeMaskDiffTree = slope.updateMask(dwDiffTree.updateMask(dwDiffTree.gt(0)));
// Define the chart and print it to the console
var chart =
ui.Chart.image.histogram({image: slopeMaskDiffTree, region: studyArea, scale: 10, maxPixels: 1e13})
.setSeriesNames(['Slope'])
.setOptions({
title: 'FABDEM Slope Histogram over Probable Landslides',
hAxis: {
title: 'Slope',
titleTextStyle: {italic: false, bold: true},
},
vAxis:
{title: 'Count', titleTextStyle: {italic: false, bold: true}},
colors: ['cf513e']
});
print(chart);Your results will look like the chart below. It shows the count of 10 m2 elevation pixels over probable landslides for our study area around Haines, Alaska. Notice how the mean value is just over 30 degrees. There are also some lower angle slopes. If you add the slopeMaskDiffTree image to the display, you will notice there are some errors in the data where the classification is incorrectly classifying ocean pixels. We could remove these by using the same masking technique to remove probable landslides in the tree class difference approach by FABDEM elevation less than or equal to 0.
These results conclude this blog. We’ve explored Dynamic World data, learned that examining probabilities over time during winter in Alaska to look at landslides is complex, that a simplified approach looking at tree class probabilities the summer before and after the Dec. 2, 2020 event was the most successful approach, and then examined the relationship of the probable landslides in Haines, Alaska to slope. We could extend our analysis by using the approach in the textbook and the USA Structures dataset in the Community Catalog to look at buildings and their topography impacted by the probable landslides.
The final blog in this series will look at these types of relationships using machine-learning approaches in Earth Engine. We will also add in the precipitation and geology data to better understand the connections between datasets and impacts of this event.
