Detecting settlements using supervised classification and deep learning

• 🔬 Data Science
• 🥠 Deep Learning and pixel-based classification

Introduction

Deep-Learning methods tend to perform well with high amounts of data as compared to machine learning methods which is one of the drawbacks of these models. These methods have also been used in geospatial domain to detect objects [1,2] and land use classification [3] which showed favourable results, but labelled input satellite data has always been an effortful task. In that regards, in this notebook we have attempted to use the supervised classification approach to generate the required volumes of data which after cleaning was used to come through the requirement of larger training data for Deep Learning model.

For this, we have considered detecting settlements for Saharanpur district in Uttar Pradesh, India. Settlements have their own importance to urban planners and monitoring them temporally can lay the foundation to design urban policies for any government.

Necessary imports

%matplotlib inline

import pandas as pd
from datetime import datetime
import matplotlib.pyplot as plt

import arcgis
from arcgis.gis import GIS
from arcgis.geocoding import geocode
from arcgis.geometry import Geometry
from arcgis.learn import prepare_data, UnetClassifier
from arcgis.raster.functions import clip, apply, extract_band, colormap, mask, stretch
from arcgis.raster.analytics import train_classifier, classify, list_datastore_content

Connect to your GIS

gis = GIS('home')

Get the data for analysis

Search for Multispectral Landsat layer in ArcGIS Online. We can search for content shared by users outside our organization by setting outside_org to True.

landsat_item = gis.content.search('title:Multispectral Landsat', 'Imagery Layer', outside_org=True)[0]
landsat = landsat_item.layers[0]
landsat_item

Multispectral Landsat
Landsat multispectral and multitemporal imagery with on-the-fly renderings and indices for visualization and analysis. The Landsat 8 and 9 imagery in this layer is updated daily and is directly sourced from the USGS Landsat collection on AWS.
Imagery Layer by esri
Last Modified: July 01, 2022
3 comments, 14,15,012 views

Search for India State Boundaries 2018 layer in ArcGIS Online. This layer has all the District boundaries for India at index - 2.

boundaries = gis.content.get('76f5897f7e354d399d261dac03af2ae3')
boundaries

India: District Boundary 2024
This layer shows Indian administrative boundaries up to the District level.
Feature Layer Collection by esri_IN_content0
Last Modified: January 03, 2024
0 comments, 5,463 views

district_boundaries = boundaries.layers[0]
district_boundaries

<MapFeatureLayer url:"https://livingatlas.esri.in/server/rest/services/IAB2024/IAB_District_2024/MapServer/0">

from arcgis.map.symbols import SimpleFillSymbolEsriSFS
from arcgis.map.symbols import SimpleLineSymbolEsriSLS
from arcgis.map.renderers import SimpleRenderer

m = gis.map('Saharanpur, India')

outline = SimpleLineSymbolEsriSLS(type = "esriSLS", style = "esriSLSSolid", color = [255, 0, 0, 255], width = 2)
symbol = SimpleFillSymbolEsriSFS(outline = outline, style = "esriSFSSolid", color = [255, 120, 0, 80])
renderer = SimpleRenderer(symbol = symbol)

m.content.add(district_boundaries, drawing_info = {"renderer":renderer}, options = {"title":"District Boundaries"})

m.legend.enabled = True
m

m.zoom = 7

As this notebook is to detect settlements for Saharanpur district, you can filter the boundary for Saharanpur.

area = geocode("Saharanpur, India", out_sr = landsat.properties.spatialReference)[0]

saharanpur = district_boundaries.query(where = "ID = '9177673'")   # query for Saharanur district boundary
saharanpur_geom = saharanpur.features[0].geometry                  # Extracting geometry of Saharanpur district boundary   

extent = dict(Geometry(saharanpur_geom).envelope)                  # Set the extent
landsat.extent = extent
saharanpur.features[0].extent = extent                     

Geometry(saharanpur_geom)

Get the training points for training the classifier. These are 212 points in total marked against 5 different classes namely Urban, Forest, Agricultural, Water and Wasteland. These points are marked using ArcGIS pro and pulished on the gis server.

data = gis.content.get('0e03348c4d4145fa9bcacc4c800c055e')
data

classificationPointsSaharanpur_Projected

Feature Layer Collection by api_data_owner
Last Modified: October 29, 2024
0 comments, 189 views

Filter satellite Imagery based on cloud cover and time duration

In order to produce good results, it is important to select cloud free imagery from the image collection for a specified time duration. In this example we have selected all the images captured between 1 October 2018 to 31 December 2018 with cloud cover less than or equal to 5% for Saharanpur region.

selected = landsat.filter_by(where="(Category = 1) AND (cloudcover <=0.05)",
                             time=[datetime.datetime(2024, 3, 1), 
                                   datetime.datetime(2024, 7, 31)],
                             geometry=arcgis.geometry.filters.intersects(area['extent']))
df = selected.query(out_fields="AcquisitionDate, GroupName, CloudCover, DayOfYear", 
                    order_by_fields="AcquisitionDate").sdf
df['AcquisitionDate'] = pd.to_datetime(df['AcquisitionDate'], unit='ms')
df

We can now select the image dated 17 April 2024 from the collection using its OBJECTID which is "5066173"

saharanpur_image = landsat.filter_by('OBJECTID=5066173') # 2024-04-17 05:18:42
saharanpur_image

m1 = gis.map('Saharanpur, India')
m1.content.add(apply(saharanpur_image, 'Natural Color with DRA'))
m1.content.add(saharanpur, options = {"title" : 'Saharanpur District'})
m1.legend.enabled = True
m1

m1.zoom = 9

A single tile completely covers our area of interest here - 5066173, if our area of interest is larger and can not be enclosed within the tile we have chosen, then a mosaic of two images can be created to enclose the AOI. Mosaicking two rasters with OBJECTID - 5035681 and 5066173 using predefined mosaic_by function.

landsat.mosaic_by(method = "LockRaster", lock_rasters = [5035681, 5066173])
mosaic = apply(landsat, 'Natural Color with DRA')
mosaic

We will clip the mosaicked raster using the boundary of Saharanpur which is our area of interest using clip function.

saharanpur_raster_clip = arcgis.raster.functions.clip(mosaic, saharanpur_geom)
saharanpur_raster_clip.extent = extent
arcgis.env.analysis_extent = extent

saharanpur_raster_clip

clip_out_name = 'saharanpur_raster_clip' + str(datetime.date.today()).replace("-", "")
saharanpur_raster_clip.save(clip_out_name, gis = gis)

saharanpur_raster_clip20241218
Analysis Image Service generated from GenerateRaster
Tiled Imagery Layer by arcgis_python
Last Modified: December 18, 2024
0 comments, 3 views

Visualize the mosaic and training points on map

m2 = gis.map('Saharanpur, India')
m2.content.add(saharanpur_raster_clip)
m2.content.add(data)
m2.legend.enabled = True
m2

m2.zoom = 9

Classification

With the Landsat layer and the training samples, we are now ready to train a support vector machine (SVM) model using train_classifier function.

saharanpur_raster_clip_search = gis.content.search(clip_out_name)
saharanpur_raster_clip_search = saharanpur_raster_clip_search[0] # Check for the name of the layer if it is the same name with which you have saved it in previous steps or else try by changing indexing
saharanpur_raster_clip_search

saharanpur_raster_clip20241218
Analysis Image Service generated from GenerateRaster
Tiled Imagery Layer by arcgis_python
Last Modified: December 18, 2024
0 comments, 3 views

classifier_definition = train_classifier(input_raster = saharanpur_raster_clip_search, 
                                         input_training_sample_json = data.layers[0].query().to_json, 
                                         classifier_parameters = {"method":"svm",
                                                                "params":{"maxSampleClass":100}},
                                         gis = gis)

Now we have the model, we are ready to apply the model to a landsat image of Saharanpur district to make prediction. This classifier will classify each pixel of the image into 5 different classes namely Urban, Forest, Agricultural, Water and Wasteland. We will use predefined classify function to perform the classification.

classified_output = classify(input_raster = saharanpur_raster_clip, 
                             input_classifier_definition = classifier_definition)
classified_output

Classify_ZVG9NR
Analysis Image Service generated from Classify
Tiled Imagery Layer by arcgis_python
Last Modified: December 18, 2024
0 comments, 4 views

Now let's visualize the classification result based on the colormap defined below.

cmap = colormap(classified_output.layers[0],
                 colormap=[[0, 255, 0, 0],
                          [1, 15, 242, 45],
                          [2, 255, 240, 10],
                          [3, 0, 0, 255],
                          [4, 176, 176, 167]])
cmap

out_name = 'saharanpur_classified_output__' + str(datetime.date.today()).replace("-", "")
cmap.save(out_name, gis = gis)

saharanpur_classified_output__20241218
Analysis Image Service generated from GenerateRaster
Tiled Imagery Layer by arcgis_python
Last Modified: December 18, 2024
0 comments, 3 views

m3 = gis.map('Saharanpur, India')
m3.content.add(cmap)
m3.legend.enabled = True
m3

m3.zoom = 10

Here, '0', '1', '2', '3' and '4' represents pixels from 'Urban', 'Forest', 'Agricultural', 'Water' and 'Wasteland' class respectively.

We can see that our model is performing reasonably well and is able to clearly identify the urban settlements while comparing the prediction with the Landsat image.

Mask settlements out of the classified raster

We have a classified raster which has each of its pixel classified against 5 different classes. But we require the pixels belonging to Urban class only, so we will mask all the pixels belonging to that class.

classified_output = gis.content.search('saharanpur_classified_output', outside_org= True)[0]
classified_output

saharanpur_classified_output__20241120
Analysis Image Service generated from GenerateRaster
Tiled Imagery Layer by arcgis_python
Last Modified: November 20, 2024
0 comments, 7 views

As we are concerned only with detecting settlement pixels, we will mask out others from the classified raster by selecting the pixels with a class value of '0' which represents 'Urban' class.

masked = colormap(mask(classified_output.layers[0],
                       no_data_values =[],
                       included_ranges =[0,0]),
                  colormap=[[0, 255, 0, 0]],
                  astype='u8')
masked.extent = extent
masked

We can save the masked raster containing the settlements to an imagery layer. This uses distributed raster analytics and performs the processing at the source pixel level and creates a new raster product.

masked = masked.save()

We can now Convert the masked pixels to polygons using 'to_features' function.

urban_item = masked.layers[0].to_features(output_name='masked_polygons' + str(datetime.datetime.now().microsecond) ,
                                          output_type='Polygon',
                                          field='Value',
                                          simplify=True,
                                          gis=gis)
urban_layer = urban_item.layers[0]

m4 = gis.map('Saharanpur, India')
m4.legend.enabled = True
m4.content.add(urban_item)
m4

m4.zoom = 11

The settlement polygons have many false positives - mainly there are areas near river bed that are confused with settlements.

These false positives were removed by post processing using ArcGIS Pro. The result obtained after cleaning up the supervised classification are shown below, which are the final detected settlements.

settlement_poly = gis.content.search("'settlement_polygons'", outside_org= True)

Deep Learning

The deep learning approach for classification needs more training data. We will use the results that we obtained after cleaning up false positives from SVM predictions above, and exported training chips against each of those polygons using ArcGIS pro to train a deep learning model.

Export training data

Export training data using 'Export Training data for deep learning' tool, detailed documentation here

• Set Landsat imagery as 'Input Raster'.
• Set a location where you want to export the training data, it can be an existing folder or the tool will create that for you.
• Set the settlement training polygons as input to the 'Input Feature Class Or Classified Raster' parameter.
• In the option 'Input Mask Polygons' we can set a mask layer to limit the tool to export training data for only those areas which have buildings in it, we created one by generating a grid of 200m by 200m on Urban Polygons layer's extent and dissolving only those polygons which contained settlements to a single multipolygon feature.
• 'Tile Size X' & 'Tile Size Y' can be set to 224.
• Select 'Classified Tiles' as the 'Meta Data Format' because we are training an 'Unet Model'.
• In 'Environments' tab set an optimum 'Cell Size'. For this example, as we have performing the analysis on Landsat8 imagery, we used 30 cell size which meant 30m on a project coordinate system.

Prepare data

We would specify the path to our training data and a few hyper parameters.

• path: path of folder containing training data.
• batch_size: No of images your model will train on each step inside an epoch, it directly depends on the memory of your graphic card. 8 worked for us on a 11GB GPU.
• imagery_type: It is a mandatory input to enable model for multispectral data processing. It can be "landsat8", "sentinel2", "naip", "ms" or "multispectral".

import os
from pathlib import Path

gis = GIS('home')

training_data = gis.content.get('b9602a576cd54bc6a8e0279a290537af')
training_data

detecting_settlements_using_supervised_classification_and_deep_learning

Image Collection by api_data_owner
Last Modified: September 03, 2020
0 comments, 41 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(filepath.split('.')[0]))

# Prepare Data
data = prepare_data(data_path,
                    batch_size=8,
                    imagery_type='landsat8')

For our model, we want out of the 11 bands in Landsat 8 image, weights for the RGB bands to be initialized using ImageNet weights whereas the other bands should have their weights initialized using the red band values. We can do this by setting the "type_init" parameter in the environment as "red band".

arcgis.env.type_init = 'red_band'

Visualize training data

The code below shows a few samples of our data with the same symbology as in ArcGIS Pro.

• rows: No. of rows we want to see the results for.
• alpha: controls the opacity of labels (Classified imagery) over the background imagery
• statistics_type: for dynamic range adjustment

data.show_batch(rows=2, alpha=0.7, statistics_type='DRA')

Here, the light blue patches show the labelled settlements.

Load model architecture

We will be using U-net [4], one of the well recognized image segmentation algorithm, for our land cover classification. U-Net is designed like an auto-encoder. It has an encoding path (“contracting”) paired with a decoding path (“expanding”) which gives it the “U” shape. However, in contrast to the autoencoder, U-Net predicts a pixelwise segmentation map of the input image rather than classifying the input image as a whole. For each pixel in the original image, it asks the question: “To which class does this pixel belong?”. U-Net passes the feature maps from each level of the contracting path over to the analogous level in the expanding path. These are similar to residual connections in a ResNet type model, and allow the classifier to consider features at various scales and complexities to make its decision.

# Create Unet Model
model = UnetClassifier(data)

Downloading: "https://download.pytorch.org/models/resnet34-b627a593.pth" to ~/.cache\torch\hub\checkpoints\resnet34-b627a593.pth
100%|█████████████████████████████████████████████████████████████████████| 83.3M/83.3M [00:07<00:00, 11.0MB/s]

Train the model

By default, the U-Net model uses ImageNet weights which are trained on 3 band imagery but we are using the multispectral support of ArcGIS API for Python to train the model for detecting settlements on 11 band Landsat 8 imagery, so unfreezing will make the backbone trainable so that it can learn to extract features from multispectral data.

Before training the model, let's unfreeze it to train it from scratch.

model.unfreeze()

Learning rate is one of the most important hyperparameters in model training. Here, we explore a range of learning rate to guide us to choose the best one.

# Find Learning Rate
model.lr_find()

0.00013182567385564074

Based on the learning rate plot above, we can see that the learning rate suggested by lr_find() for our training data is 0.0002754228703338166. We can use it to train our model.

We would fit the model for 100 epochs. One epoch mean the model will see the complete training set once.

# Training
model.fit(epochs=100, lr=0.00022908676527677726);

model.mIOU()

{'0': 0.9867301422601041, '1': 0.6107132199252355}

Here, we can see loss on both training and validation set decreasing, which shows how well the model generalizes on unseen data thereby preventing overfitting. Also, as we can see our model here is able to classify pixels to settlement or non settlement with an accuracy 0.98. Once we are satisfied with the accuracy, we can save the model.

Save the model

We would now save the model which we just trained as a 'Deep Learning Package' or '.dlpk' format. Deep Learning package is the standard format used to deploy deep learning models on the ArcGIS platform. UnetClassifier models can be deployed by the tool 'Classify Pixels Using Deep Learning' available in ArcGIS Pro as well as ArcGIS Enterprise. For this sample, we will be using this model in ArcGIS Pro to extract settlements.

We will use the save() method to save the model and by default it will be saved to a folder 'models' inside our training data folder itself.

model.save('100e')

Computing model metrics...

WindowsPath('~/AppData/Local/Temp/detecting_settlements_using_supervised_classification_and_deep_learning/models/100e')

Load an intermediate model to train it further

If we need to further train an already saved model we can load it again using the code below and repeat the earlier mentioned steps of training.

# Load Model from previous saved files 
# model.load('100e')

Preview results

The code below will pick a few random samples and show us ground truth and respective model predictions side by side. This allows us to validate the results of your model in the notebook itself. Once satisfied we can save the model and use it further in our workflow.

model.show_results(rows=2)

Deploy model and extract settlements

The model saved in previous step can be used to extract classified raster using 'Classify Pixels Using Deep Learning' tool. Further the classified raster is regularised and finally converted to a vector Polygon layer. The regularisation step uses advanced ArcGIS geoprocessing tools to remove unwanted artifacts in the output.

Generate a classified raster

In this step, we will generate a classified raster using 'Classify Pixels Using Deep Learning' tool available in both ArcGIS Pro and ArcGIS Enterprise.

• Input Raster: The raster layer from which we want to extract settlements from.
• Model Definition: It will be located inside the saved model in 'models' folder in '.emd' format.
• Padding: The 'Input Raster' is tiled and the deep learning model classifies each individual tile separately before producing the final 'Output Classified Raster'. This may lead to unwanted artifacts along the edges of each tile as the model has little context to predict accurately. Padding as the name suggests allows us to supply some extra information along the tile edges, this helps the model to predict better.
• Cell Size: Should be close to at which we trained the model, we specified that at the Export training data step .
• Processor Type: This allows to control whether the system's 'GPU' or 'CPU' would be used in to classify pixels, by 'default GPU' will be used if available.
• Parallel Processing Factor: This allows us to scale this tool, this tool can be scaled on both 'CPU' and 'GPU'. It specifies that the operation would be spread across how many 'cpu cores' in case of cpu based operation or 'no of GPU's' in case of GPU based operation.

Output of this tool will be in form of a 'classified raster' containing both background and Settlements. The pixels having background class are shown in green color and pixels belonging to settlement class are in purple color.

                                        A subset of preditions by our Model

Change the symbology of the raster by setting the second class of the classified pixels to 'No color'.

After this step, we will get a masked raster which will look like the screenshot below.

We can now run 'Raster to Polygon' tool in ArcGIS Pro to convert this masked raster to polygons.

Results obtained from Deep Learning approach are ready to be analyzed.

             Deep Learning Results

Conclusion

This notebook shows how traditional algorithms like Support Vector Machines can be used to generate the necessary quality data for Deep Learning models which after cleaning can serve as the training data required by a DL model for generating accurate results at a much larger extent.

References

[1] https://developers.arcgis.com/python/sample-notebooks/detecting-swimming-pools-using-satellite-image-and-deep-learning/

[2] https://developers.arcgis.com/python/sample-notebooks/extracting-building-footprints-from-drone-data/

[3] https://developers.arcgis.com/python/sample-notebooks/land-cover-classification-using-unet/

[4] Olaf Ronneberger, Philipp Fischer, Thomas Brox: U-Net: Convolutional Networks for Biomedical Image Segmentation, 2015; arXiv:1505.04597.