Classify land cover to measure shrinking lakes

Compare imagery to calculate area change in Lake Aculeo Lagoon, Chile.

Introduction

One of the important freshwater waterbodies of central Chili, Lake Aculeo Lagoon shrank over time due to climatic drought, drying of surface tributaries and streams that provided a continuous flow of water and high existing demand of ground water [1]. The lagoon dried completely in May of 2018. Those whose livelihoods depend on the lake are alarmed, as the shrinking lake changes the land cover of the area and impacts the economy. In order to help them identify and save their waterbodies, we'll generate land cover and compare imagery between 2016 and 2017 to quantify the surface area of the lake and show changes over time.

Manually extracting features from raw data, such as generating land cover maps, is time consuming. Deep learning automates the process and minimizes the manual interaction necessary to complete these tasks. To quantify the change in lake surface area from 2016 to 2017, we'll classify the land cover in both images, identifying the areas covered with water and distinguishing them from other land cover, such as vegetation or urban areas.

This sample aims to demonstrate how ArcGIS pretrained models can be used to generate landcover from imageries of different time periods for further analysis. Once the land covers are created, we will use arcgis.raster submodule of ArcGIS API for Python for calculating surface area reduced from 2016 to 2017.

Note: This sample is supported with ArcGIS Image for ArcGIS Online. For more details read here.

Necessary imports

import arcgis
from arcgis.gis import GIS
from arcgis.map import Map
from arcgis.learn import classify_pixels
from arcgis.raster.functions import equal_to
from arcgis.raster.functions import extract_band


import pandas as pd
from datetime import datetime as dt
from ipywidgets import HBox, VBox, Label, Layout
gis = GIS('home')

Get data and model for classification

Search for Multispectral Sentinel-2 Level 1-C imagery for year 2016 and 2017 on the Copernicus Open Access Hub. Download imagey of the desired time period and publish them in ArcGIS Online. The code below uses the get method to get items.

sentinel2016 = gis.content.get("a264d1eaa8bd4ff2a0a1e27b59daa7d3")
sentinel2017 = gis.content.get("ea192d03d1324c62a2cb91e26e3a0ece")

Search for the Land Cover Classification (Sentinel-2) deep learning package in ArcGIS Online.

model = gis.content.get('afd124844ba84da69c2c533d4af10a58')
model
Land Cover Classification (Sentinel-2)
Deep learning model to perform land cover classification on Sentinel-2 imagery.
Deep Learning Package by esri_analytics
Last Modified: August 01, 2024
40 comments, 59,220 views

This model is trained on Corine Land Cover (CLC) 2018 data that contains 16 classes. The screenshot below shows the color code for each land cover class.

Visually compare Lake Aculeo Lagoon over time

As all the bands cannot be depicted at the same time, we usually pick a combination of three bands that we display through the color channels red, green, and blue, which can be seen by the human eye.

sentinel2016_rgb = extract_band(sentinel2016.layers[0], [4, 3, 2])
sentinel2017_rgb = extract_band(sentinel2017.layers[0], [4, 3, 2])

Let's compare the imagery visually to get a sense of how the shape of the lake has evolved over time.

map1 = Map()
map1.content.add(sentinel2016_rgb)
map2 = Map()
map2.content.add(sentinel2017_rgb)
map1.sync_navigation(map2)
from ipywidgets import *

map1.layout = Layout(flex = '1 1', padding = '10px')
map2.layout = Layout(flex = '1 1', padding = '10px')

# Create VBoxes for each map and label
box1 = VBox([Label("Sentinel-2016"), map1], layout=Layout(width='50%'))
box2 = VBox([Label("Sentinel-2017"), map2], layout=Layout(width='50%'))

# Place the VBoxes side by side using an HBox
hbox = HBox([box1, box2])

# Display the HBox
hbox
<PIL.PngImagePlugin.PngImageFile image mode=RGBA size=1993x642>

The imagery shows a bright and clear distinction between the blue lake and the green vegetation nearby.

map1.zoom_to_layer(sentinel2016_rgb)

Classify land cover in 2016

To quantify the change in lake surface area from 2016 to 2017, we'll classify the land cover in both images, identifying the areas covered with water and distinguishing them from other land cover, such as vegetation or urban areas. In multispectral imagery, such as Sentinel-2, every individual pixel (or cell) in the image has a value for every spectral band. As we can see from the vibrant imagery of Lake Poyang, there are many possible color values for all varieties of shades and hues. However, all the pixels representing the same land cover tend to have somewhat similar spectral values. By classifying the image, we'll identify the pixels that are similar in value and group them together to represent a small number of classes, such as water, vegetation, or urban areas.

ext = map1.extent  # desired extent for generating land cover
classified_out2016 = classify_pixels(
    input_raster=sentinel2016.layers[0],
    model=model,
    model_arguments={
        "padding": 56,
        "batch_size": 64,
        "predict_background": True,
        "tile_size": 224,
        "test_time_augmentation": True,
        "merge_classes": True,
        "sentinel_imagery_level": 1,
    },
    output_name="classified_output_2016" + str(dt.now().microsecond),
    context={"extent": ext, "processorType": "GPU", "cellSize": 10},
    tiles_only=False,
)
classified_out2016
classified_output_2016150954
Analysis Image Service generated from ClassifyPixelsUsingDeepLearning
Imagery Layer by demos_deldev
Last Modified: September 19, 2024
0 comments, 3 views

Let's visualize the 2016 raster and its land cover on the map.

map3 = Map()
map3.content.add(sentinel2016_rgb)
map4 = Map()
map4.content.add(classified_out2016)
map3.sync_navigation(map4)
map3.layout = Layout(flex = '1 1', padding = '10px')
map4.layout = Layout(flex = '1 1', padding = '10px')

# Create VBoxes for each map and label
box1 = VBox([Label("Sentinel-2016"), map3], layout=Layout(width='50%'))
box2 = VBox([Label("Classified Raster 2016"), map4], layout=Layout(width='50%'))

# Place the VBoxes side by side using an HBox
hbox = HBox([box1, box2])

# Display the HBox
hbox
<PIL.PngImagePlugin.PngImageFile image mode=RGBA size=1992x635>
map3.zoom_to_layer(sentinel2016_rgb)

Classify land cover in 2017

Next, We will follow the same process to generate land cover for the 2017 raster.

classified_out2017 = classify_pixels(
    input_raster=sentinel2017.layers[0],
    model=model,
    model_arguments={
        "padding": 56,
        "batch_size": 64,
        "predict_background": True,
        "tile_size": 224,
        "test_time_augmentation": True,
        "merge_classes": True,
        "sentinel_imagery_level": 1,
    },
    output_name="classified_output_2017" + str(dt.now().microsecond),
    context={"extent": ext, "processorType": "GPU", "cellSize": 10},
    tiles_only=False,
)
classified_out2017
classified_output_2017563397
Analysis Image Service generated from ClassifyPixelsUsingDeepLearning
Imagery Layer by demos_deldev
Last Modified: September 19, 2024
0 comments, 4 views
map5 = Map()
map5.content.add(sentinel2017_rgb)
map6 = Map()
map6.content.add(classified_out2017)
map5.sync_navigation(map6)
map5.layout = Layout(flex = '1 1', padding = '10px')
map6.layout = Layout(flex = '1 1', padding = '10px')

# Create VBoxes for each map and label
box1 = VBox([Label("Sentinel-2017"), map5], layout=Layout(width='50%'))
box2 = VBox([Label("Classified Raster 2017"), map6], layout=Layout(width='50%'))

# Place the VBoxes side by side using an HBox
hbox = HBox([box1, box2])

# Display the HBox
hbox
<PIL.PngImagePlugin.PngImageFile image mode=RGBA size=1993x638>
map6.legend.enabled = True
map5.zoom_to_layer(sentinel2017_rgb)

Clean up the classification

We'll now clean up the classified images with generalization analysis tools to remove minor water bodies around the lake. We'll also smooth the lake's boundaries.

We are only interested in pixels that have been clasified as water, so we will run an equal_to tool to produce a raster that contains pixels classified as water.

water_mask2016 = equal_to(
    rasters=[
        classified_out2016.layers[0],
        51,
    ],  # 51 is the code value for Inland Waters class.
    astype="U16",
).save("water_mask2016" + str(dt.now().microsecond))
water_mask2017 = equal_to(
    rasters=[classified_out2017.layers[0], 51], astype="U16"
).save("water_mask2017" + str(dt.now().microsecond))

Next, we will convert the raster into polygon features in order to calculate area of desired lake polygon.

lake_poly_2016 = water_mask2016.layers[0].to_features()
lake_poly_2017 = water_mask2017.layers[0].to_features()

All features with gridcode 0 have nodata values, while gridcode 1 has features with water polygons.

df_2016 = lake_poly_2016.layers[0].query(as_df=True, where="gridcode= 1")
df_2017 = lake_poly_2017.layers[0].query(as_df=True, where="gridcode= 1")

The code below calculates the area of all of the features and then returns the area of largest water polygon, the lake.

Calculate area over time

area_lost = (
    df_2016.SHAPE.geom.area.max() - df_2017.SHAPE.geom.area.max()
)  # area in square meters
area_lost
1300486.944107932
area_lost/1000000 #area in square kilometers
1.300486944107932

Conclusion

In this sample, we compared visually and classified Sentinel-2 imagery of Lake Aculeo Lagoon to understand how much the lake's area has changed over time. These findings indicate a severe problem: the lake has lost around 1.3 square kilometers in only 1 year. These calculations don't reveal the causes of Lake Aculeo Lagoon's reduction, however, they do provide factual evidence of a serious problem and provide a starting point for environmental scientists, and others, to conduct further research.

References

Your browser is no longer supported. Please upgrade your browser for the best experience. See our browser deprecation post for more details.