Machine Learning for Change Detection: A Technical Deep Dive

This blog post is part four of a four-part series, “Machine Learning for Change Detection”. 

Our series, “Machine Learning for Change Detection” explores changes to the Greater Austin, Texas area with a combination of:

• Planet’s satellite data 
• Regrid’s parcel data 
• Our change detection maps

Part four is the final part in this series. 

In this technical deep dive, we’ll walk through our workflow and give some insight into how these Open Source tools can be used within the context of Esri software. 

To further explore this demo, you can:

• Follow the YouTube video below which begins at Step 3,
• Follow our condensed instructions in this blog post,
• Access the entire change detection demo on GitHub. 

Let’s dive in.

Data Preparation

1. Downloading the Planet Data

To start this project, we access Planet’s developer resources and API to find images of the Greater Austin area. 

If you want to run the demo with the same imagery we used, you need a Planet API key with permissions to the study area.

Run download_data.py to download the imagery: 

This is broken down into two calls because we don’t want to exceed the rate limit.

2. Creating Mosaics

Next, using GDAL on each directory of rasters, we can produce a VRT and subsequently a mosaic for each year. 

3. Loading in Data (YouTube Video Starts Here)

To load the rasters, we use rasterio and then chunk the data into 8,192 x 8,192 chunks using rioxarray.

We drop the infrared band and preserve the RGB channels.

We then scale the data between 0 and 1 by dividing each channel by 255 (since the data type is INT8).

Next, we align and stack the two rasters together to create a new multi-dimensional xarray DataArray that holds all the data in the following shape:

(2 (time) x 3 (RGB) x columns (x) x rows (y))

4. Normalizing the Data

To normalize the data, we subtract the mean for each band and divide by the standard deviation. 

This is a standard Machine Learning technique that scales data to have zero mean and zero unit variance. The model can then use this new range to predict meaningful color differences across time.

Next, we check if the chosen spatial size will fit properly onto our GPU (we used a NVIDIA A100 GPU) and rechunk the data accordingly. 

This way each block of the data consists of data from both rasters and can be forward propagated through the pre-trained PyTorch change detection model.

Greater Austin area, 2017 and 2022 rasters. Satellite imagery: Planet Labs PBC.

Predicting Geographic Changes for the Entire Scene

Now that our data is normalized, we can pass the data into our pre-trained model. 

For each Dask chunk, we split it into small chips and pass each chip to the GPU for prediction. 

Our output is a binary mask for each chip where values close to 255 indicate areas of change and values close to 0 indicate stationary geographic regions. 

We stack these chips together to make a single mask for the entire scene— this mask is our predictions.

5. Predicting Chips

We load one chip at a time using Dask’s map_blocks (more details on this function in Step 6, Dask map_blocks Function), 

→ Forward propagate that sample through the model, 

→ Offload the data back to the CPU, 

→ Convert it from a Torch Tensor to a NumPy array representation to be stored within an xarray DataArray. 

How copy_and_predict_chunked Works

This function handles the passing of the pre-trained change detection model as well as a one chunk of our data at a time to a Dask-CUDA worker for prediction. 

We split up each chunk into chips that are the input size required by the pre-trained change detection model. 

Each NumPy array is converted to a Torch Tensor. 

A dictionary, gpu_chip, holds the required input data for the change detection model and is instantiated where: 

• “A” is the 2017 sample, 
• “B” is the 2022 sample, 
• “L” is an empty Torch Tensor (since we don’t have any labels for changes that have occurred). 

This chip is placed on the GPU, input to the predict_chips function, and forward propagated through the change detection model. 

The returned NumPy array is the predicted change mask at the same spatial window as the input.

6. Parallelizing Using Dask 

To parallelize our analysis using Dask, we initialize an array of predictions (predictions_array) which calls Dask’s map_blocks function on each chunk of the data within our combined dataset (ds_comb).

Dask map_blocks Function

Runs our copy_and_predict_chunked function on each chunk of our dataset in parallel. 

The output is initialized as a NumPy array of data type UINT8. The size of the output will be the same spatial size as that of the input, drop_axis=[0,1]. 

We input our loaded model that has been waiting on the client, and then wrap this NumPy array as an xarray DataArray.

Visualizing Results

The entire Planet dataset contains approximately 140 km x 120 km of data. 

This takes several minutes to process. 

For our demo, we use a small 4,000 px x 4,000 px snippet of the data (12 km x 12 km) and visualize change detection results. 

To run the full dataset, simply run:

7. Predicting a subset of the data

Next, let’s set the sample size to 4,000 pixels and grab a patch around the center of the dataset.

We will run dask.compute on this subset.

Dask Dashboard Running

We make a figure 1 row x 3 columns that shows our 2017 sample data, predicted change mask, and 2022 sample data.

8. Visual Spot Checking

Let’s visually check that the change mask is corresponding to visual changes in the data.

We can zoom in on areas with the largest changes by leveraging connected component analysis. 

First, let’s find continuous regions within the change mask that are areas of predicted change, predictions_to_regions.

Next, we’ll use skimage functions to extract the regions of interest and perform bookkeeping with Pandas to store the centroids of large areas of change. 

We can then use the coordinates for each region to return shapely Points in case we want to plot where the changes occurred using localize_regions.

We grab patches of our raster data based on the centroids of the regions in grab_patches.

Next, we’ll plot the 2017 data, change maps, and 2022 data from our subset using visual_spot_check.

9. Classifying Changes

Now that we verified our pre-trained change detection model is accurately detecting changes that occurred, we’ll use a classifier to predict a land cover type for each patch of interest. 

We trained a ResNet50 model via fast.ai on the EuroSAT dataset and we’ll use this model to describe the changes.

Let’s load our trained ResNet50 model using load_learner.

Next, we will predict a label for each patch of area where changes were detected for 2017 and 2022 using classify_patches.

• Each patch is transformed to a PIL image of a Torch Tensor to be input to the model.
• The most confident prediction is retained for each sample.
• The label of the class is returned.
• The results are sorted by largest area. This way we can visualize large changes rather than any small changes, which may stem from slight color differences rather than significant geographic changes.
• Predictions are visualized by plotting the areas where changes were detected in 2017, the change mask, and the area from 2022.

This summarizes in broad strokes how we approached this analysis at makepath.

If you’re interested in getting into deeper technical issues or have an idea of another use case for this workflow, please contact us. 

How to Leverage Machine Learning for Change Detection using Esri Software

We completed this project using Open Source Machine Learning and Open Source GIS tools.

But what if you wanted to repeat this workflow using Esri products? 

You can find Esri tools similar to what we used in a few different places:

• ArcGIS API for Python: Through this API you can install deep learning frameworks for ArcGIS by using the arcgis.learn module. This API can also help you classify land cover using deep learning. 
• ArcGIS API for Javascript: This API can help you create interactive web apps like our time sliders using the Swipe widget.

If you’re interested in learning about more ways to use Python within ArcGIS, check out our blog post on Superpowered GIS: Esri’s ArcGIS + Open Source Spatial Analysis Tools.

Series Summary

Part One: Change detection map for an area of 6,750 sq. mi. using Planet’s satellite imagery.

Part Two: Zoomed in time sliders to visualize changes to specific areas.

Part Three: Predictions of future development in Austin by comparing our change maps with Regrid’s parcel data.

Part Four: Technical deep dive on how to use Open Source Machine Learning tools for land cover change detection.

Share Your Projects with Us

If you replicate our project successfully, or have questions about our workflow:

• Share your results with us on Twitter @makepathGIS
• Email us at contact@makepath.com

References

1. https://planetarycomputer.microsoft.com/docs/tutorials/landcover/ – A lot of our inspiration came from this tutorial located on the Planetary Computer.
2. Satellite imagery: Planet Labs PBC.
3. Hao Chen, Z., & Zhenwei Shi (2021). Remote Sensing Image Change Detection with Transformers. IEEE Transactions on Geoscience and Remote Sensing, 1-14.
4. Helber, P., Bischke, B., Dengel, A., & Borth, D. (2018). Introducing EuroSAT: A Novel Dataset and Deep Learning Benchmark for Land Use and Land Cover Classification. In IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium (pp. 204–207).
5. Helber, P., Bischke, B., Dengel, A., & Borth, D. (2019). Eurosat: A novel dataset and deep learning benchmark for land use and land cover classification. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing.