Introduction
Image segmentation and classification are very important topics in GIS and remote sensing applications. Both are approaches to extracting features from imagery based on objects. Segmentation groups pixels in close proximity and having similar spectral characteristics into a segment, which doesn't need any training data and is considered as unsupervised learning. In contrast, image classification is a type of supervised learning which classifies each pixel to a class in the training data. In this guide, we are going to demonstrate both techniques using ArcGIS API for Python.
import arcgis
from arcgis import GIS
from arcgis.raster.analytics import *
from arcgis.features import FeatureSet, FeatureCollectiongis = GIS(url='https://pythonapi.playground.esri.com/portal', username='arcgis_python', password='amazing_arcgis_123')Classification
In this example, we are going to perform a land cover classification using a Landsat image in Iowa and hand labeled training data. In the training data, there are four classes in total: Developed Area, Forest, Planted/Cultivated, and Water.
First, let's draw the training data on a map and visualize it.
items1 = gis.content.search("training_samples", item_type="Feature Layer")
items1[0]
Feature Layer Collection by api_data_owner
map1 = gis.map('Iowa')
map1
defined_extent = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
'xmin': -10357358.3716515,
'ymin': 5063187.75726388,
'xmax': -10340267.8481094,
'ymax': 5075204.7581259}map1.extent = defined_extent
map1.content.add(items1[0])map1.legend.enabled = TrueWe need to convert the training sample feature class into JSON format, which is the format the image classification requires.
query_result = items1[0].layers[0].query()
training_sample_json = query_result.to_jsonNow let's search for Landsat data that will be used for image classification.
items2 = gis.content.search("title: multispectral landsat", item_type="Imagery Layer", outside_org=True)
items2[0]
Imagery Layer by esri_livingatlas
landsat_layer = items2[0].layers[0]Set the analysis extent to a small region in Iowa so we don't process all landsat images.
arcgis.env.verbose = True
arcgis.env.analysis_extent = defined_extentWe can retrieve a landsat scene by OBJECTID.
filtered_layer = landsat_layer.filter_by('OBJECTID=3134102')With the Landsat layer and the training samples, we are now ready to train a support vector machine (SVM) model.
classifier_definition = train_classifier(input_raster=filtered_layer,
input_training_sample_json=training_sample_json,
classifier_parameters={"method":"svm",
"params":{"maxSampleClass":500}},
gis=gis)Submitted.
Executing...
Start Time: Tuesday, November 19, 2019 10:02:48 AM
Running script TrainClassifier...
Prepare training feature.
DEBUG: {'method': 'svm', 'params': {'maxSampleClass': 500}}
Training...
Succeeded at Tuesday, November 19, 2019 10:03:31 AM (Elapsed Time: 43.48 seconds)
If we print out classifier_definition, we can see that it contains information such as number of classes and coefficients that define the SVM model.
classifier_definition["Definitions"][0]["Classes"][{'ClassValue': 10,
'ClassName': 'Water',
'Red': 84,
'Green': 117,
'Blue': 168},
{'ClassValue': 20,
'ClassName': 'Developed',
'Red': 232,
'Green': 209,
'Blue': 209},
{'ClassValue': 40,
'ClassName': 'Forest',
'Red': 133,
'Green': 199,
'Blue': 126},
{'ClassValue': 80,
'ClassName': 'Planted / Cultivated',
'Red': 251,
'Green': 246,
'Blue': 93}]Now we have the model, we are ready to apply the model to a landsat image to make prediction.
from arcgis.raster.analytics import classify
classified_output = classify(input_raster=filtered_layer,
input_classifier_definition=classifier_definition)
classified_outputSubmitted.
Executing...
Start Time: Tuesday, November 19, 2019 10:04:38 AM
Running script Classify...
Image service {'name': 'Classify_FHAX9B', 'serviceUrl': 'https://datascienceqa.esri.com/server/rest/services/Hosted/Classify_FHAX9B/ImageServer'} already existed.
Classifying...
Finished
Updating service with data store URI.
Getting image service info...
Updating service: https://dev0002085.esri.com:6443/arcgis/admin/services/Hosted/Classify_FHAX9B.ImageServer/edit
Update item service: https://dev0002085.esri.com:6443/arcgis/admin/services/Hosted/Classify_FHAX9B.ImageServer successfully.
Succeeded at Tuesday, November 19, 2019 10:04:49 AM (Elapsed Time: 10.79 seconds)
Imagery Layer by portaladmin
Now let's visualize the classification result and compare it with the original Landsat image.
map2 = gis.map('Iowa') # shows the Landsat layer
map2
map2.extent = defined_extent
map2.content.add(filtered_layer)map2 = gis.map('Iowa') # shows the classification result on the top of Landsat layer
map2
map3.extent = defined_extent
map3.content.add(filtered_layer)
map3.content.add(classified_output.layers[0])We can see that our model is making reasonable prediction results by comparing the prediction with the Landsat image. For example, water has been clearly identified.
Segmentation
The goal of image segmentation is to identify segments in your imagery by grouping adjacent pixels together that have similar spectral and spatial characteristics. arcgis.raster.analytics.segment uses an algorithm called Mean Shift which uses a moving window that calculates an average pixel value to determine which pixels should be included in each segment. As the window moves over the image, it iteratively recomputes the value to make sure that each segment is suitable. The result is a grouping of image pixels into a segment characterized by an average color.
In this example, we are going to use image segmentation to extract center-pivot in a chosen region in Saudi Arabia. First, let's search for a 8-bit landsat layer we've published beforehand. Please note that our image segmentation requires the input raster layer to be 8-bit.
landsat_pivot_items = gis.content.search("title: landsat_pviot_345_8bit", item_type="Imagery Layer")
landsat_pivot_items[0]
Imagery Layer by api_data_owner
landsat_pivot_layer = landsat_pivot_items[0].layers[0]map3 = gis.map('Saudi Arabia')
map3
map2.content.add(landsat_pivot_layer);The characteristics of the image segments depend on three parameters: spectral detail, spatial detail, and minimum segment size. These parameters determine whether a pixel should be part of a segment or not. You can vary the amount of detail that characterizes a feature of interest. For example, if you are more interested in impervious features than in individual buildings, adjust the spatial detail parameter to a small number; a lower number results in more smoothing and less detail.
- The segment function provides multi-band support. You can define which 3 bands are used in segmentation. Default is [0,1,2].
spectral_detailparameter is used to set the level of importance given to the spectral differences of features in your imagery. Default is 15.5.spatial_detail, on the other hand, sets the level of importance given to the proximity between features in your imagery. Default is 15.
from arcgis.raster.analytics import segment
seg_output = segment(input_raster=landsat_pivot_layer) # it uses band [0, 1, 2] by default
seg_outputmap3 = gis.map('Saudi Arabia')
map3
map3.extent = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100},
'type': 'extent',
'xmax': 4735079.27946308,
'xmin': 4720420.87480519,
'ymax': 3226614.75178954,
'ymin': 3214862.51053094}
map3.content.add(seg_output)Now we can see most of the center pivots have been correctly extracted. Notice that by default the function uses band [0, 1, 2] for segmentation and renders the result using the same colormap.