Landsat 8, MATLAB, and the California Drought (Part II)

In the first part of this two part series we looked at how to measure the water level of Cachuma Lake using the Image Processing Toolbox in MATLAB and a Landsat 8 image. We were able to produce a graph showing the decrease in surface area during the California Drought of 2012–2017. In this part I want to produce a time-lapse movie directly visualizing this decrease to really bring home just how bad the drought was, and to show that we are once again experiencing below average rainfall.

We start with our cropped sequence of images; here are two samples:

Two sample images from 2014 (left) and 2016 (right) from our sequence.

There are a couple of issues we need to address with these two images before we would be able to combine them into a movie. First, the brightness and contrast do not match, and secondly the lakes are not aligned in the same position.

I first selected a single image to be my “fixed” reference image. There is nothing really special about this image other than that the lake is quite centered on it. I then go through the process of adjusting and aligning all other images to this reference image.

We can use imhistmatch to autocorrect the brightness and contrast issues:

With histograms matched to provide consistent brightness and contrast.

For alignment I first tried imregister, but had only very limited success. Often the images would not be anywhere close to aligned, despite trying a variety of options. Fortunately, the overall scale of the images is always the same, and there never seemed to be rotation, so we really only need to correct for x-y translation. I decided to just compute the 2D cross-correlation between the corresponding binarized images (as created in Part I), find the point at maximum intensity, then translate the original image by this amount.

% First create binarized images
BWFixed = im2single(cachumabinarize(fixedImage));
BWMoving = im2single(cachumabinarize(movingImage));
% Find 2d cross correlation
XC = xcorr2(BWMoving, BWFixed);
% Find the linear index for the max val
[C n] = max(XC(:));
% Find row, column indices
[y x] = ind2sub(size(XC), n);
% Translate the color image by the computed offset
I = imtranslate(I, -([x y] - size(BWFixed)));

The cross-correlation slides the moving image across the fixed image and measures the amount of overlapping white pixels. It is most intense at the position of maximum overlap.

Finally we crop and resize each image:

I = imresize(imcrop(I, [144 256 511 287]), 2);

Our final images:

We can use MATLAB’s VideoWriter to assemble the images into a movie

(Pause the movie and use Shift-Arrow Key to move frame by frame)