This notebook is part of the collection accompanying the talk "Analyzing Satellite Images With Python Scientific Stack" by Milos Miljkovic given at PyData NYC 2014. The content is BSD licensed.
Landsat 8 has two onoboard imaging instruments:
OLI collects data from 9 imaging bands. Five bands cover the range from violet to red in the visible spectrum, and 4 cover the range from near infrared to short infrared wavelengths. Detailed information is available at Operational Land Imager
TIRS acquisition is limited to two long wave infrared bands. Detailed information is available at Thermal Infrared Imaging Sensor
All raw images from Landsat 8 satellite have one feature in common - they appear to be dark. The reason for this is that all imaging sensors have a wide dynamic range in order to capture scenes ranging from dark colored oceans to brightly lit snow. Histogram equalization technique is used to process images in order to increase brightness and contrast. Since material presented here is intended for intermediate user, please refer to this article explaining histogram equalizationn.
The following sections will explain name, spectral range, spatial resolution, and use of all 11 bands in a typical Landsat 8 satellite image analysis.
# To install watermark extension execute:
# %install_ext https://raw.githubusercontent.com/HyperionAnalytics/watermark/master/watermark.py
%load_ext watermark
%watermark -v -m -p numpy,scikit-image,matplotlib
CPython 3.4.2 IPython 2.3.1 numpy 1.9.1 scikit-image 0.10.1 matplotlib 1.4.2 compiler : GCC 4.2.1 (Apple Inc. build 5577) system : OS X 10.10.1 release : 14.0.0 machine : x86_64 processor : Intel(R) Core(TM) i7-3740QM CPU @ 2.70GHz CPU cores : 8 interpreter: 64bit
import re
import numpy as np
from skimage import io, exposure
from matplotlib import pyplot as plt
%matplotlib inline
def read_band(n):
"""
Load Landsat 8 band
Input:
n - integer in the range 1-11
Output:
img - 2D array of uint16 type
"""
if n in range(1, 12):
tif_list = !ls *.TIF
band_name = 'B' + str(n) + '.TIF'
img_idx = [idx for idx, band_string in enumerate(tif_list) if band_name in band_string]
img = io.imread(tif_list[img_idx[0]])
return img
else:
print('Band number has to be in the range 1-11!')
def image_show(img, title):
"""
Show image
Input:
img - 2D array of uint16 type
title - string
"""
fig = plt.figure(figsize=(10, 10))
fig.set_facecolor('white')
plt.imshow(img, cmap='gray')
plt.title(title)
plt.show()
def image_histogram(img):
"""
Plot image histogram
Input:
img - 2D array of uint16 type
"""
co, ce = exposure.histogram(img)
fig = plt.figure(figsize=(10, 7))
fig.set_facecolor('white')
plt.plot(ce[1::], co[1::])
plt.show()
def image_adjust_brightness(img, limit_left, limit_right, title):
"""
Adjust image brightness and plot the image
Input:
img - 2D array of uint16 type
limit_left - integer
limit_right - integer
title - string
"""
img_ha = exposure.rescale_intensity(img, (limit_left, limit_right))
fig = plt.figure(figsize=(10, 10))
fig.set_facecolor('white')
plt.imshow(img_ha, cmap='gray')
plt.title(title)
plt.show()
def get_gain_bias(n):
"""
Get band radiance gain and bias
Input:
n - integer in the range 1-11
Output:
gain - float
bias - float
"""
if n in range(1, 12):
n_str = str(n)
s_g = 'RADIANCE_MULT_BAND_' + n_str + ' = '
s_b = 'RADIANCE_ADD_BAND_' + n_str + ' = '
fn = !ls *_MTL.txt
f = open(fn[0], 'r+')
search_str_g = '(?<=' + s_g + ').*'
search_str_b = '(?<=' + s_b + ').*'
for line in f:
s1 = re.search(search_str_g, line)
s2 = re.search(search_str_b, line)
if s1:
gain = float(s1.group(0))
elif s2:
bias = float(s2.group(0))
f.close()
return gain, bias
else:
print('Band number has to be in the range 1-11!')
def get_thermal_const(n):
"""
Get band thermal constants K1 and K2
Input:
n - integer in the range 10-11
Output:
Ks - list of two floats
"""
if n in range(10, 12):
n_str = str(n)
s = 'K._CONSTANT_BAND_' + n_str + ' = '
search_str = '(?<=' + s + ').*'
fn = !ls *_MTL.txt
f = open(fn[0], 'r+')
Ks = []
for line in f:
k = re.search(search_str, line)
if k:
Ks.append(float(k.group(0)))
f.close()
return Ks
else:
print('Band number has to be in the range 10-11!')
This band has a dual use:
atmospheric correction procedures.
Landsat 8 is unique in possessing the capability to acquire data in the violet part of the visible spectrum.
In this example, we'll use data set LC82040522013123LGN01 showing the mouth of the Rio Geba in Guinea Bissau, on the west coast of Africa.
cd /Users/kronos/gis/l8/guinea_bissau/
/Users/kronos/gis/l8/guinea_bissau
b1 = read_band(1)
image_adjust_brightness(b1, 10000, 13500, 'Band 1, data set LC82040522013123LGN01')
These three bands correspond to blue, green, and red colors from the visible spectrum. Processing of natural color images will be discussed in detail in a separate notebook titled "color_image_processing".
This band covers the near infrared (NIR) part of the spectrum. NIR measurements are important for agriculture and ecology because the water in healthy plants efficiently scatters the wavelengths in this range. Used in conjunction with red band, it is possible to estimate the state of green plants using the normalized difference vegetation index (NDVI). NDVI will be discussed in detail in a separate notebook titled "ndvi_calculation".
In this example, we'll use data set LC82040522013123LGN01 showing the mouth of the Rio Geba in Guinea Bissau, on the west coast of Africa. Lush vegetation is visible as bright areas on the river banks and islands.
cd /Users/kronos/gis/l8/guinea_bissau/
/Users/kronos/gis/l8/guinea_bissau
b5 = read_band(5)
image_adjust_brightness(b5, 6000, 23000, 'Band 5, data set LC82040522013123LGN01')
The bands cover two different portions of the shortwave infrared (SWIR) spectrum. Bands 6 and 7 are used to assess soil water content, and also distinguish rocks which look very similar in the images formed using other bands. The two bands can also be used to acquire clear images through smoke.
In this example, we'll use data set LC80430332014262LGN00 showing King Fire, which burned in August 2014 around the Lake Tahoe. Image shows no sign of smoke, which is visible in red color image (Band 4).
cd /Users/kronos/gis/l8/fire/
/Users/kronos/gis/l8/fire
b6 = read_band(6)
b4 = read_band(4)
image_adjust_brightness(b6, 4000, 30000, 'Band 6, data set LC80430332014262LGN00')
image_adjust_brightness(b4, 8000, 22000, 'Band 4, data set LC80430332014262LGN00')
This band is called panchromatic due to the fact it collects all the visible spectrum colors at the same time, combining them into one channel. Panchromatic band has the highest spatial resolution. Use of panchromatic band to produce color pan-sharpened images will be discussed in a separate notebook titled "panchromatic_sharpening".
Unlike all the other bands, this band is using part of the spectrum where the atmospheric absorption is strong due to the presence of water vapor. Cirrus clouds are found at high altitudes, they are made up of ice crystals, and appear as thin, wispy strands. Band 9 enables detection of cirrus clouds which can be hard to spot in satellite images.
In this example, we'll use data set LC81990262013104LGN01 acquired over Paris. The image shows many criss-crossing lines, and all the lines are contrails from heavy air traffic.
cd /Users/kronos/gis/l8/paris/
/Users/kronos/gis/l8/paris
b9 = read_band(9)
image_adjust_brightness(b9, 5000, 7500, 'Band 9, data set LC81990262013104LGN01')
These two bands are used for thermal mapping of the ground and improved estimation of soil moisture. Using formula derived from the Plank's law and two thermal constants provided in dataset metadata file, it is possible to calculate the ground temperature using the formula:
$$T = \frac{K_{2}}{\ln (\frac{K_{1}}{L_{\lambda}} + 1)}$$where:
In this example, we'll use data set LC80080292014065LGN00 acquired over the Bay of Fundu in Nova Scotia. In the bay waters, we can see swirls from warmer bay water mixing with colder water from the ocean.
cd /Users/kronos/gis/l8/bay/
/Users/kronos/gis/l8/bay
b10 = read_band(10)
gain, bias = get_gain_bias(10)
Ks = get_thermal_const(10)
L = gain * b10 + bias # Convert raw pixel readouts to spectral radiance
temp = Ks[1] / np.log(Ks[0] / L + 1)
image_adjust_brightness(temp, 255, 273, 'Band 10, brightness temperature,'
'data set LC80080292014065LGN00')