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, Layoutgis = 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
Deep Learning Package by esri_analytics
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
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 coverclassified_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
Imagery Layer by demos_deldev
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
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_out2017map5 = 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
map6.legend.enabled = Truemap5.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_lost1300486.944107932
area_lost/1000000 #area in square kilometers1.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.