How to detect changes between two GeoTIFF files using NumPy

1   Introduction

Try these exercises first:

• How to read a GeoTIFF file as a NumPy array using GDAL
• Exercise: Plot a GeoTIFF file
• Exercise: Calculate the mean cell value of a GeoTIFF file
• How to extract cell values from a GeoTIFF file using NumPy
• How to classify a GeoTIFF file using NumPy

In this exercise, we’ll try to detect changes between two GeoTIFF files. Use nlcd2001_clipped.tif and nlcd2016_clipped.tif. For np.where(), see this exercise.

2   Python

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

def read_geotiff(filename):
    ds = gdal.Open(filename)
    band = ds.GetRasterBand(1)
    arr = band.ReadAsArray()
    del ds, band
    return arr

nlcd01_arr = read_geotiff("nlcd2001_clipped.tif")
nlcd16_arr = read_geotiff("nlcd2016_clipped.tif")

nlcd_changes = nlcd16_arr - nlcd01_arr

nlcd_changed = np.where(nlcd_changes != 0, 1, 0)

plt.subplot(221)
plt.imshow(nlcd01_arr)

plt.subplot(222)
plt.imshow(nlcd16_arr)

plt.subplot(223)
plt.imshow(nlcd_changes)

plt.subplot(224)
plt.imshow(nlcd_changed)

plt.show()

3   Explanation

Import necessary modules and define the read_geotiff() function.

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt

def read_geotiff(filename):
    ds = gdal.Open(filename)
    band = ds.GetRasterBand(1)
    arr = band.ReadAsArray()
    del ds, band
    return arr

Read data as 2-dimensional NumPy arrays.

nlcd01_arr = read_geotiff("nlcd2001_clipped.tif")
nlcd16_arr = read_geotiff("nlcd2016_clipped.tif")

NumPy supports very simple algebraic operations between arrays. In this case, we want to subtract the old 2001 NLCD class values from the new 2016 ones. If there were no changes between the two years i some locations, the resulting cell values at those no-change locations would be 0 because we subtract class values from class values. If there were changes, those cells would be assigned positive or negative values.

nlcd_changes = nlcd16_arr - nlcd01_arr

Let’s create a binary array that contains 1 for changed cells and 0 for the rest. We use np.where() and the nlcd_changes != 0 condition. What this expression means is that if nlcd_changes is not 0 (NLCD changed), use 1. Otherwise (NLCD didn’t change), use 0.

nlcd_changed = np.where(nlcd_changes != 0, 1, 0)

Since nlcd01_arr, ncd16_arr, nlcd_changes, and nlcd_changed are all 2-dimensional arrays, we can use plt.imshow() to plot them.

plt.subplot(221)
plt.imshow(nlcd01_arr)

plt.subplot(222)
plt.imshow(nlcd16_arr)

plt.subplot(223)
plt.imshow(nlcd_changes)

plt.subplot(224)
plt.imshow(nlcd_changed)

plt.show()