Landuse classification

Download this project.


Image Classification

Written by Philipp Rudiger
Created: November 30, 2018
Last updated: August 3, 2021

Satellite images often need to be classified (assigned to a fixed set of types) or to be used for detection of various features of interest. Here we will look at the classification case, using labelled satellite images from various categories from the UCMerced LandUse dataset. scikit-learn is useful for general numeric data types, but it doesn't have significant support for working with images. Luckily, there are various deep-learning and convolutional-network libraries that do support images well, including Keras (backed by TensorFlow) as we will use here.

In [1]:
import intake

import numpy as np
import holoviews as hv
import pandas as pd
import random

from holoviews import opts
hv.extension('bokeh')
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
/tmp/ipykernel_6273/3533684752.py in <module>
----> 1 import intake
      2 
      3 import numpy as np
      4 import holoviews as hv
      5 import pandas as pd

ModuleNotFoundError: No module named 'intake'

Get the classes and files

All of the labeled image classification data is in a public bucket on s3 with a corresponding intake catalog. This catalog provides several ways to access the data. You can access one image by landuse and id, access all the images for a given landuse, or access all the images.

In [ ]:
cat = intake.open_catalog('https://s3.amazonaws.com/earth-data/UCMerced_LandUse/catalog.yml')
list(cat)

The first time you run the cell below it will download all the images which takes about 3 minutes on my machine. After that, the images are cached and it'll take about 100 ms.

In [ ]:
%time da = cat.UCMerced_LandUse_all().to_dask()
da

We can see what's going on more easily if we convert the data to a dataset with each data variable representing a different landuse.

In [ ]:
ds = da.to_dataset(dim='landuse')
ds

Split files into train and test sets

In order to accurately test the performance of the classifier we are building we will split the data into training and test sets with an 80/20 split. We randomly sample the images in each category and assign them either to the training or test set:

In [ ]:
train_set = np.random.choice(ds.id, 80, False)
test_set = np.setdiff1d(ds.id, train_set)

Define function to sample from train or test set

In [ ]:
landuses = da.landuse.data
landuse_list = list(landuses)

Next we define a function that randomly samples an image, either from the training or test set and ant

In [ ]:
def get_sample(landuse=None, set='training'):
    landuse = landuse or np.random.choice(landuses)
    i = random.choice(train_set if set == 'training' else train_set)
    return ds[landuse].sel(id=i)

def plot(data):
    options = opts.RGB(xaxis=None, yaxis=None)
    title = '{}: {}'.format(data.name, data.id.item())
    plot = hv.RGB(data.data.compute())
    return plot.options(options).relabel(title)

We can inspect the data on one of these samples to see that the data is loaded as an xarray.DataArray.

In [ ]:
data = get_sample()
data

We can plot this array as a holoviews RGB image so we can visualize it:

In [ ]:
plot(data)
In [ ]:
hv.Layout(list(map(plot, map(get_sample, np.random.choice(landuses, 4))))).cols(2)

Define the model

Now it's time to define a model. The code snippet below defines a convolutional neural network containing a stack of 3 convolution layers with a ReLU activation and followed by max-pooling layers. This is very similar to the architectures that Yann LeCun advocated in the 1990s for image classification (with the exception of ReLU).

In [ ]:
from keras.models import Sequential
from keras.layers import Conv2D, MaxPooling2D
from keras.layers import Activation, Dropout, Flatten, Dense

size = (150, 150)

model = Sequential()
model.add(Conv2D(32, (3, 3), input_shape=(*size, 3)))
model.add(Activation('relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))

model.add(Conv2D(32, (3, 3)))
model.add(Activation('relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))

model.add(Conv2D(64, (3, 3)))
model.add(Activation('relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))

model.add(Flatten())  # this converts our 3D feature maps to 1D feature vectors
model.add(Dense(64))
model.add(Activation('relu'))
model.add(Dropout(0.5))
model.add(Dense(21))
model.add(Activation('softmax'))

model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])

Declare the data

The dataset of images we have is relatively small so to avoid overfitting we will want to apply some augmentation to it. The code below defines a generator that randomly samples from either the training or test set and generates a randomly rotated and cropped 150x150 window onto the 256x256 images along with a label array:

In [ ]:
def get_array(array, size=(150, 150)):
    # Randomly flip
    if random.getrandbits(1):               # equivalent to random.choice([True, False]) but much faster
        array = array.transpose(1, 0, 2)
    if random.getrandbits(1):
        array = array[::-1]
    if random.getrandbits(1):
        array = array[:, ::-1]
    # Randomly crop
    sh, sw = size
    h, w = array.shape[:2]
    b = np.random.randint(h-sh)
    l = np.random.randint(w-sw)
    array = array[b:b+sh, l:l+sw]
    return array/255.

# set up an mapping to an identity matrix to use for one-hot encoding
one_hot_mapping = dict(zip(landuses, np.eye(21)))

def gen_samples(set='training', labels=None):
    "Generates random arrays along with landuse labels"
    while True:
        choice = get_sample(set=set)
        if labels is not None:
            labels.append(choice.name)
        one_hot = one_hot_mapping[choice.name]
        data = choice.data.compute()
        yield get_array(data)[np.newaxis, :], one_hot[np.newaxis, :]        

Run the model

Before we start running the model let's set up a keras Callback to build a dashboard to monitor the accuracy and loss during training:

In [ ]:
import time
from keras.callbacks import Callback

class MonitorCallback(Callback):
    """
    Builds a streaming dashboard to monitor the accuracy
    and loss during training using HoloViews streams.
    """
    
    _format = '%s - Epoch: %d - Elapsed time: %.2fs'
    
    def __init__(self, metrics=['acc', 'loss']):
        super().__init__()
        sample = {'Epoch': np.array([])}
        for metric in metrics:
            sample[metric] = np.array([])
        self.buffer = hv.streams.Buffer(sample)
        dmaps = []
        for metric in metrics:
            def cb(data, metric=metric):
                return hv.Curve(
                    data, 'Epoch', metric, label=self._format
                    % (metric, self.epoch, self.elapsed_time))
            dmap = hv.DynamicMap(cb, streams=[self.buffer])
            dmaps.append(dmap)
        self.layout = hv.Layout(dmaps)
        self.metrics = metrics
        self.start_time = None
        self.epoch = 0

    def on_train_begin(self, logs={}):
        self.start_time = time.time()
        
    @property
    def elapsed_time(self):
        if self.start_time is None:
            return 0
        else:
            return time.time() - self.start_time

    def on_epoch_end(self, epoch, logs=None):
        self.epoch += 1
        data = {'Epoch': [self.epoch]}
        for metric in self.metrics:
            data[metric] = [logs.get(metric)]
        self.buffer.send(data)

Now we can create an instance of the callback and display the dashboard, which will at first appear blank.

In [ ]:
monitor = MonitorCallback()
monitor.layout.options('Curve', width=450)

Next we will fit the model with our training data generator, as the model is running it will update the dashboard above:

In [ ]:
%%time
history = model.fit_generator(gen_samples('training'), steps_per_epoch=50, epochs=500, verbose=False, callbacks=[monitor])

Evaluate the model

First we will have a look at the monitoring output but smooth it slightly so we can make out the overall trend:

In [ ]:
from holoviews.operation.timeseries import rolling
monitor.layout.options('Curve', width=400).map(rolling, hv.Curve)

Now let us test the predictions on the test set, first visually:

In [ ]:
def set_title_color(color, *args):
    """Helper function to set title color"""
    args[0].handles['plot'].title.text_color = color
In [ ]:
def get_prediction(cls):
    sample = get_sample(cls, 'test')
    array = get_array(sample.data.compute())[np.newaxis, ...]
    p = model.predict(array).argmax()
    p = landuses[p]
    return (plot(sample)
            .relabel('Predicted: %s - Actual: %s' % (p, cls))
            .options(hooks=[lambda *args: set_title_color('red' if p!=cls else 'blue', *args)]))

options = dict(fontsize={'title': '8pt'}, width=250, height=250)
hv.Layout([get_prediction(landuse).options(**options) for landuse in landuses]).cols(3)

And now numerically by running 1,000 predictions on the test set:

In [ ]:
ntesting = 1000
labels = []
test_gen = gen_samples('test', labels)
prediction = model.predict_generator(test_gen, steps=ntesting)
y_pred = landuses[prediction.argmax(axis=1)]
y_true = np.array(labels[:ntesting])

accuracy = (y_pred==y_true).sum()/ntesting

print(f'Accuracy on test set {accuracy}')

Next we can see how well the classifier performs on the different categories. To see how it performs we will make 50 predictions on each category and record both the accuracy and the predictions:

In [ ]:
def predict(cls, iterations=50):
    accurate, predictions = [], []
    for i in range(iterations):
        sample = get_sample(cls, 'test')
        array = get_array(sample.data.compute())[np.newaxis, ...]
        p = model.predict(array).argmax()
        p = landuses[p]
        predictions.append(p)
        accurate.append(p == cls)
    return np.sum(accurate)/float(iterations), predictions

accuracies = [(c, *predict(c)) for c in landuses]

We can now break down the accuracy by landuse category:

In [ ]:
df = pd.DataFrame(accuracies, columns=['landuse', 'accuracy', 'predictions'])

hv.Bars(df, 'landuse', 'accuracy', label='Accuracy by Landuse Category').options(
    width=700, xrotation=45, color='landuse', 
    cmap='Category20', show_legend=False)

Another interesting way of viewing this data is to look at which categories the classifier got confused on. We will count how many times the classifier classified one category as another category and visualize the result as a Chord graph where each edge is colored by the actual category. By clicking on a node we can reveal which other categories incorrectly identified an image as being of that category:

In [ ]:
pdf = pd.DataFrame([(p, l) for (_, l, _, ps) in df.itertuples() for p in ps], columns=['Prediction', 'Actual'])
graph = pdf.groupby(['Prediction', 'Actual']).size().to_frame().reset_index()
confusion = graph.rename(columns={0: 'Count'})

hv.Chord(confusion).relabel('Confusion Graph').options(
    node_color='index', cmap='Category20', edge_color='Actual', labels='index',
    width=600, height=600)

Clicking on buildings, for instance, reveals a lot of confusion about overpasses, mediumresidential, and intersections, all of which do share visual features in common. Conversely, number of buildings were misidentified as parklots, which is also reasonable. As we saw in the bar chart above, forests on the other hand, have lots of edges leading back to itself, demonstrating the high accuracy observed for that category of images.

In [ ]:
confusion.Count /= 50
hv.HeatMap(confusion, label='Confusion Matrix').sort().options(
    xrotation=45, width=500, height=500, cmap='blues', tools=['hover'], invert_yaxis=True, zlim=(0,1))
This web page was generated from a Jupyter notebook and not all interactivity will work on this website. Right click to download and run locally for full Python-backed interactivity.

Download this project.