Classify land cover to measure shrinking lakes

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

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

Input
import arcgis
from arcgis.gis import GIS
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
Input
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.

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

Search for the Land Cover Classification (Sentinel-2) deep learning package in ArcGIS Online. We can search for content shared by users outside our organization by setting outside_org to True.

Input
model = gis.content.search(
    "Land Cover Classification (Sentinel-2) owner: esri_analytics",
    "Deep Learning Package",
    outside_org=True,
)[0]
model
Output
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: October 26, 2021
13 comments, 5,021 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.

Input
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.

Input
map1 = gis.map()
map1.add_layer(sentinel2016_rgb)
map2 = gis.map()
map2.add_layer(sentinel2017_rgb)
Input
map1.sync_navigation(map2)
Input
hbox_layout = Layout()
hbox_layout.justify_content = "space-around"
hb1 = HBox([Label("Sentinel-2016"), Label("Sentinel-2017")])
hb1.layout = hbox_layout
VBox([hb1, HBox([map1, map2])])
Output

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

Input
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.

Input
ext = map1.extent  # desired extent for generating land cover
Input
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,
    },
    output_name="classified_output_2016" + str(dt.now().microsecond),
    context={"extent": ext, "processorType": "GPU", "cellSize": 10},
    tiles_only=False,
)
classified_out2016
Output
classified_output_2016261561
Analysis Image Service generated from ClassifyPixelsUsingDeepLearningImagery Layer by demos_deldev
Last Modified: November 28, 2021
0 comments, 0 views

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

Input
map3 = gis.map()
map3.add_layer(sentinel2016_rgb)
map4 = gis.map()
map4.add_layer(classified_out2016)
Input
map3.sync_navigation(map4)
Input
hbox_layout = Layout()
hbox_layout.justify_content = "space-around"
hb1 = HBox([Label("Sentinel-2016"), Label("Classified Raster 2016")])
hb1.layout = hbox_layout
VBox([hb1, HBox([map3, map4])])
Output
Input
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.

Input
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,
    },
    output_name="classified_output_2017" + str(dt.now().microsecond),
    context={"extent": ext, "processorType": "GPU", "cellSize": 10},
    tiles_only=False,
)
classified_out2017
Output
classified_output_201752976
Analysis Image Service generated from ClassifyPixelsUsingDeepLearningImagery Layer by demos_deldev
Last Modified: November 28, 2021
0 comments, 0 views
Input
map5 = gis.map()
map5.add_layer(sentinel2017_rgb)
map6 = gis.map()
map6.add_layer(classified_out2017)
Input
map5.sync_navigation(map6)
Input
hbox_layout = Layout()
hbox_layout.justify_content = "space-around"
hb1 = HBox([Label("Sentinel-2017"), Label("Classified Raster 2017")])
hb1.layout = hbox_layout
VBox([hb1, HBox([map5, map6])])
Output
Input
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.

Input
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.

Input
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.

Input
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

Input
area_lost = (
    df_2016.SHAPE.geom.area.max() - df_2017.SHAPE.geom.area.max()
)  # area in square meters
area_lost
Output
1300486.9210555563
Input
area_lost/1000000 #area in square kilometers
Output
1.3004869210555563

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.