Landcover mapping using hyperspectral imagery and deep learning

In this notebook, we will use hyperspectral data to train a deep learning model and will see if the model can extract subclasses of two LULC classes: developed areas and forests. Hyperion imagery is used in the current analysis to classify the types of forests and developed areas. The data can be downloaded from USGS earth explorer.

Hyperion data preparation

The Earth Observing-1 (EO-1) satellite was launched November 21, 2000 as a one-year technology demonstration/validation mission. After the initial technology mission was completed, NASA and the USGS agreed to the continuation of the EO-1 program as an Extended Mission. The EO-1 Extended Mission is chartered to collect and distribute Hyperion hyperspectral and Advanced Land Imager (ALI) multispectral products according to customer tasking requests.

Hyperion collects 220 unique spectral channels ranging from 0.357 to 2.576 micrometers with a 10-nm bandwidth. The instrument operates in a pushbroom fashion, with a spatial resolution of 30 meters for all bands. The standard scene width is 7.7 kilometers.

Download hyperion imagery

Login to the earth explorer using the USGS credentials. Then, select the Address/Place option in the Geocoding Method dropdown and write the name or address of the area of interest.

Draw a polygon over the area of interest and change the Cloud Cover Range to 0% - 10%. Click on Results in the bottom left.

Next, expand EO-1, select EO-1 Hyperion, and click on Results.

Download the data from the product list.

Pre processing of hyperion imagery

Remove bad bands

Not all bands are useful for analysis. Bad bands are primarily water vapor bands consists inforation about atmosphere, that cause spikes in the reflectance curve. List of bad bands of the Hyperion sensor, L1R product are: 1 to 7 - Not illuminated

58 to 78 - Overlap region

120 to 132 - Water vapor absorption band

165 to 182 - Water vapor absorption band

185 to 187 - Identified by Hyperion bad band list

221 to 224 - Water vapor absorption band

225 to 242 - Not illuminated

Bands with vertical stripping - (8-9, 56-57, 78-82, 97-99, 133-134, 152-153, 188, 213-216, 219-220)

Create composite rasters

Two sets of composite rasters were created using Composite bands function. The band range of two raster composites are as follows: Two sets of composite rasters were created using the Composite bands function function. The band range of the two raster composites are as follows:

NIR & visible composite - 9 to 55

SWIR composite - 82-96, 100-119, 135-164, 183-184, 189-212, 217-218

Click on Imagery tab -> click on Raster Functions -> Search "Composite Bands" -> Click on browse button and select the bands -> After selecting the bands, click on Create new layer.

Convert DN values to at-sensor radiance

Scale factors will be used to convert the pixel DN value to at-sensor radiance. A scale factor of 40 will be used for the VNIR composite raster and a scale factor of 80 will be used for the SWIR composite. The rasters will be divided by the scale factor using the Raster Calculator tool in ArcGIS Pro.

Create combined composite raster

Combine both the VNIR and SWIR radiance composite rasters using the "Composite Band" function.

Eliminate NoData pixels

Create null raster

The output raster has a 0 pixel value for NoData that is shown as a black pixel. The NoData pixels can be eliminated using the Set Null tool.

Use the following parameters:

Raster: combined_composite

False raster = 0

Cellsize Type = Max Of

Extent Type = Intersection Of

Clip

The Clip tool will be used to eliminate instances of NoData from the combined composite raster using the null raster that was created in the previous step.

The input parameters are as follows:

Raster: combined_composite

Clipping Type: Outside

Clipping Geometry/Raster: null_raster

Principal Component Analysis

Principle component analysis can be used to determine the comparatively higher information bands of Hyperion's 224 total bands. Principal Components performs Principal Component Analysis (PCA) on a set of raster bands and generates a single multiband raster as its output. The value specified for the number of principal components determines the number of principal component bands in the output multiband raster. The number must not be larger than the total number of raster bands in the input.

The input parameters are as follows:

Input raster bands: Clip_combined_composite

Number of Principal components: 12 (12 Principal component were chosen because they cover more than 98% of the total variance)

The output raster will have 12 bands representing 12 Principal components.

Export training data for deep learning model

from arcgis.gis import GIS
gis = GIS('home')
ent_gis = GIS('https://pythonapi.playground.esri.com/portal', 'arcgis_python', 'amazing_arcgis_123')

The raster generated from the Hyperion imagery through the Principal Components tool will be used as the input raster for the training data.

pca_raster = ent_gis.content.get('0d86c3f9dedc4566a7cb3d2c18d1ffef')
pca_raster

input_pca_raster_hyperion
input_pca_raster_hyperionImagery Layer by api_data_owner
Last Modified: March 21, 2021
0 comments, 0 views

The following feature layer will be used as label for training the model.

lulc_polygons = ent_gis.content.get('f69011a01c39494bade9b30708bce064')
lulc_polygons

label_hyperion
lulc_hyperionFeature Layer Collection by api_data_owner
Last Modified: March 20, 2021
0 comments, 6 views

m = gis.map('USA', 4)
m.add_layer(lulc_polygons)
m.legend=True
m

The above layer consists of 7 classes. Developed areas have three sub classes: High, Medium, and Low density areas and Forested areas have 3 types of forests: Evergreen, Mixed, and Deciduous forest. No DATA represents all other LULC classes.

m.zoom_to_layer(lulc_polygons)

A feature layer will be used as a input mask polygon while exporting the training data to define the extent of the study area.

mask_poly = ent_gis.content.get('bfe7bb7ea8454d9384cb8005a5c0f855')
mask_poly

mask_polygon_hyperion
mask_polygon_hyperionFeature Layer Collection by api_data_owner
Last Modified: March 22, 2021
0 comments, 0 views

The Export training data for deep learning tool is used to prepare training data for training a deep learning model. The tool is available in both ArcGIS Pro and ArcGIS Enterprise.

Model training

This step will utilize Jupyter Notebooks, and documentation on how to install and setup the environment is available here.

Necessary Imports

import os
from pathlib import Path

import arcgis
from arcgis.learn import prepare_data, UnetClassifier

Get training data

We have already exported the data, and it can be directly downloaded using the following steps:

training_data = gis.content.get('0a8ff4ce9f734bb5bb72104aea8526af')
training_data

landcover_classification_using_hyperspectral_imagery_and_deep_learning
Image Collection by api_data_owner
Last Modified: March 22, 2021
0 comments, 2 views

filepath = training_data.download(file_name=training_data.name)

import zipfile
with zipfile.ZipFile(filepath, 'r') as zip_ref:
    zip_ref.extractall(Path(filepath).parent)

data_path = Path(os.path.join(os.path.splitext(filepath)[0]))

Visualize training data

The prepare_data function takes a training data path as input and creates a fast.ai databunch with specified transformation, batch size, split percentage, etc.

data = prepare_data(data_path,  
                    batch_size=8)

To get a sense of what the training data looks like, use the show_batch() method to randomly pick a few training chips and visualize them. The chips are overlaid with masks representing the building footprints in each image chip.

data.show_batch(alpha=1)

Load model architecture

arcgis.learn provides the UnetClassifier model for per pixel classification that is based on a pretrained convnet, like ResNet, that acts as the 'backbone'. More details about UnetClassifier can be found here.

model = UnetClassifier(data, pointrend=True)

Train the model

Learning rate is one of the most important hyperparameters in model training. We will use the lr_find() method to find an optimum learning rate at which we can train a robust model.

lr = model.lr_find()

We are using the suggested learning rate above to train the model for 2000 epochs.

model.fit(100, lr=lr)

We have trained the model for a further 1900 epochs to improve model performance. For the sake of time, the cell below is commented out.

# model.fit(1900)

Visualize classification results in validation set

It's a good practice to see results of the model vis-a-vis ground truth. The code below picks random samples and shows us ground truths and model predictions, side by side. This enables us to preview the results of the model within the notebook.

model.show_results()

Evaluate model performance

As we have 7 classes for this segmentation task, we need to perform an accuracy assessment for each class. To achieve this, ArcGIS API for Python provides the per_class_metrics function that calculates a precision, recall, and f1 score for each class.

model.per_class_metrics()

Save model

We will save the model that we trained as a 'Deep Learning Package' ('.dlpk' format). A Deep Learning Package is the standard format used to deploy deep learning models on the ArcGIS platform.

We will use the save() method to save the trained model. By default, it will be saved to the 'models' sub-folder within our training data folder.

model.save('unet_2000e', publish=True, gis=gis)

Computing model metrics...

WindowsPath('C:/data/2021/hyperspectral/hyperion_l8lulc_128px_64strd_30px_mask_frst_bltp_7classes/models/unet_2000e')

Model inference

Using ArcGIS Pro, we can use the trained model on a test image/area to classify different types of built-up areas and forests in the hyperspectral satellite image.

After we have trained the UnetClassifier model and saved the weights for classifying images, we can use the Classify Pixels Using Deep Learning tool, available in both ArcGIS pro and ArcGIS Enterprise, for inferencing at scale.

inference_raster = ent_gis.content.get('a69ec3006c9a45a49b85cc43aa529e44')
inference_raster

pca_raster_for__inferencing
pca_raster_for_inferencingImagery Layer by api_data_owner
Last Modified: March 20, 2021
0 comments, 1 views

with arcpy.EnvManager(processorType="GPU"): out_classified_raster = arcpy.ia.ClassifyPixelsUsingDeepLearning("pca_raster_for_inferencing", r"C:\path\to\model.dlpk", "padding 32;batch_size 2;predict_background True; tile_size 128", "PROCESS_AS_MOSAICKED_IMAGE", None); out_classified_raster.save(r"C:\sample\sample.gdb\inferenced_results")

Results visualization

The classified output raster is generated using ArcGIS Pro. The output raster is published on the portal for visualization.

inferenced_results = ent_gis.content.get('03609691b9e74def954f81b9285dee93')
inferenced_results

inferenced_results_
inferenced_resultsMap Image Layer by api_data_owner
Last Modified: March 20, 2021
0 comments, 0 views

rgb_imagery = ent_gis.content.get('f9a5e23541b845da8846eb54ea561594')
rgb_imagery

fcc_imagery_l8
fcc_imagery_l8Imagery Layer by api_data_owner
Last Modified: March 20, 2021
0 comments, 2 views

Create map widgets

Two map widgets are created showing the RGB imagery and Inferenced results.

m1 = gis.map()
m1.add_layer(rgb_imagery)
#m1.basemap = 'satellite'
m2 = gis.map()
m2.add_layer(inferenced_results)

Synchronize web maps

The maps are synchronized with each other using MapView.sync_navigation functionality. Map view synchronization helps in comparing the inferenced results with the RGB imagery. A detailed description about advanced map widget options can be referenced here.

m1.sync_navigation(m2)

Set the map layout

from ipywidgets import HBox, VBox, Label, Layout

Hbox and Vbox were used to set the layout of map widgets.

hbox_layout = Layout()
hbox_layout.justify_content = 'space-around'

hb1=HBox([Label('FCC imagery'),Label('Results')])
hb1.layout=hbox_layout

Results

The resulting predictions are provided as a map for better visualization.

VBox([hb1,HBox([m1,m2])])

Inferenced results legend

m2.zoom_to_layer(inferenced_results)

In the map widgets, it can be seen that the model is able to differentiate between forests and open spaces (gardens, parks, etc). The different types of developed areas are also accurately classified by the model, with the bright blue-gray patches representing high density developed areas, for instance. The forested areas have also been accurately classified by the model.

Limitation

One of the major limitations in working with hyperspectral data is its lack of availability for many parts of the world. In comparison to multispectral imagery tiles, like those from Landsat-8, one tile of hyperion imagery covers a relatively small area. If a training dataset is prepared with multiple tiles, a model can be trained and more accurate results can be generated.

Conclusion

In this notebook, we have covered a lot of ground. In Part 1, we covered how Hyperion data can be downloaded using USGS Earth Explorer, the pre-processing steps of Hyperion, how to create a Principal Component raster, and how to export training data for deep learning using ArcGIS Pro. In Part 2, we demonstrated the steps to prepare the input data, train a pixel-based classification model, visualize the results, and generate the accuracy metrics. Finally, in Part 3, we demonstrated the process of inferencing results on a test raster/area. The same workflow can be used to identify minerals/rocks, plant species, oil spills, etc. using hyperspectral imageries and arcgis.learn models.

References

• https://developers.arcgis.com/python/guide/how-unet-works/