In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Lecture 25:

  • Learn how to use satellite imagery to understand land usage.
  • Learn how to use patches in matplotlib.
  • Learn about using the scikit-learn python package to classify data.

Using Landsat Satellite Data

The landsat series of satellites provide satellite imagery for the entire Earth's surface. They are frequently used to study land management. Landsat data vary depending on which satellite is being used. Today we will be using data from Landsat Series 7. Landsat satellites record the brightness of reflected light in multiple different wavelength 'bands'. This means that the satellite gets a separate value for different colors of light.

Band Wavelength (micrometers) Color Resolution (m)
1 0.45-0.52 Blue 30
2 0.52-0.60 Green 30
3 0.63-0.69 Red 30
4 0.77-0.90 Near Infrared 30
5 1.55-1.75 Infrared 1 30
6 10.40-12.50 Thermal 60
7 2.09-2.35 Infrared 2 30
8 0.52-0.90 Panchromatic 15

For this lecture, we're going to look at putting labels on the landsat data for the area around the town of Bellingham, Washington. The data come from the USGS Global Land Survey dataset for this region from Sep 25th 2009 https://earthexplorer.usgs.gov/.

In the datasets folder, we have a set of TIF images for each of the bands in the landsat dataset. Each image contains the level of brightness in the given color band for each picture. To get an understanding, let's read in one of our TIF files, and plot the data using plt.imread( ) and plt.imshow( ).

In [2]:
B4=plt.imread('Datasets/Landsat_Data/B4.TIF') #Read in the TIF file #filter the TIF file data for the region of interest
plt.figure(figsize=(16,12)) #make  a big figure object
plt.imshow(B4,cmap='Greys_r') #Plot the data for band 4 as a greyscale image.
plt.axis('Off'); #Turn off the plotting axes with tick marks, etc.

This gives us the brightness seen just in the 'Near Infrared' band. We can definitely see the land here and ocean here, but it's hard to tell what's going on in the landscape and we don't just want our satellite imagery in black and white - we want color!.

For simplicity's sake, let's ignore bands 6 and 8 which have different resolutions and just read in the images for the other 5 color bands.

In [3]:
B1=plt.imread('Datasets/Landsat_Data/B1.TIF')
B2=plt.imread('Datasets/Landsat_Data/B2.TIF')
B3=plt.imread('Datasets/Landsat_Data/B3.TIF')
B5=plt.imread('Datasets/Landsat_Data/B5.TIF')
B7=plt.imread('Datasets/Landsat_Data/B7.TIF')

Let's combine some of these bands into a False color image. We can use different combinations of the bands by making a 1000x1000x3 array and passing it to plt.imshow( ). Each of the 1000x1000 elements in the array represents a pixel, within that we have 3 values for the amount of Red, Green or Blue.

For our False color image, we're going to use band 5 (infrared) for "Red", band 4 (near infrared) for "Green" and band 3 (red) for "Blue". We will use these bands as it looks a bit nicer than a true color image (for which we would use bands 3, 2 and 1 - but you could try it both ways!).

Note that we will at times be using the .T array method which is essentially identical to the .transpose( ) array method in the following.

In [4]:
RGB=np.array([B5.T,B4.T,B3.T]).T #Make an RGB image from bands 5, 4 and 3
plt.figure(figsize=(16,12)) #Plot a big figure
plt.imshow(RGB);#Plot the RGB array
#plt.axis('Off'); # we need to see the coordinates now