Generally, multispectral imagery is preferred for Landuse Landcover (LULC) classification, due to its high temporal resolution and high spatial coverage. With the advances in remote sensing technologies, hyperspectral images are now also a good option for LULC classification, due to their high spectral resolution. The main difference between multispectral and hyperspectral imagery, is the number of bands and how narrow those bands are. One of the advantages hyperspectral sensors have over multispectral sensors is the ability to differentiate within classes. For instance, due to the high spectral information content that creates a unique spectral curve for each class, hyperspectral sensors can distinguish by tree or crop species.
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.
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.
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.
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:
False raster = 0
Cellsize Type = Max Of
Extent Type = Intersection Of
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:
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.
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('1bb227a836094127b530cf8db3a7fa4d') pca_raster
The following feature layer will be used as label for training the model.
lulc_polygons = ent_gis.content.get('f69011a01c39494bade9b30708bce064') lulc_polygons
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.
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
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
This step will utilize Jupyter Notebooks, and documentation on how to install and setup the environment is available here.
import os from pathlib import Path import arcgis from arcgis.learn import prepare_data, UnetClassifier
We have already exported the data, and it can be directly downloaded using the following steps:
training_data = gis.content.get('0a8ff4ce9f734bb5bb72104aea8526af') training_data
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)))
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.
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)
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.
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.
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.
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.
|NoData||Deciduous Forest||Developed, High Intensity||Developed, Low Intensity||Developed, Medium Intensity||Evergreen Forest||Mixed Forest|
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...
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('df7781fe3d904afa914ef8804b08c802') inference_raster
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")
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('29a95588665245359d3787e48fc501cc') inferenced_results
rgb_imagery = ent_gis.content.get('94d505a8fd084e7f896ba46247b12739') rgb_imagery
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.
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
The resulting predictions are provided as a map for better visualization.
Inferenced results legend
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.
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.
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.