Landsat 8 Color Image Processing

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.



In [1]:

    
# To install watermark extension execute:
# %install_ext https://raw.githubusercontent.com/HyperionAnalytics/watermark/master/watermark.py
%load_ext watermark



In [2]:

    
%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



In [3]:

    
import numpy as np
from skimage import io, exposure
from matplotlib import pyplot as plt
from IPython.display import Image



In [4]:

    
%matplotlib inline



In [5]:

    
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!')



In [6]:

    
def color_image_show(img, title):
    """
    Show image
    Input:
    img - 3D array of uint16 type
    title - string
    """
    fig = plt.figure(figsize=(10, 10))
    fig.set_facecolor('white')
    plt.imshow(img/65535)
    plt.title(title)
    plt.show()



In [7]:

    
cd /Users/kronos/gis/l8/guinea_bissau/









    



/Users/kronos/gis/l8/guinea_bissau

Load RGB Bands

Landsat 8 bands corresponding to RGB are 4-3-2. Data is loaded as 2D uint16 arrays, then stacked into NumPy 3D array of the same type, and imaged.
Please note that matplotlib cannot display RGB uint16 arrays. Data is transformed to float64 in range 0-1.



In [8]:

    
b2 = read_band(2)
b3 = read_band(3)
b4 = read_band(4)

img432 = np.dstack((b4, b3, b2))
%reset_selective -f b



In [9]:

    
color_image_show(img432, '4-3-2 image, data set LC82040522013123LGN01')

Histogram Equalization

Histograms of RGB colors corresponding to raw data show that the data is not utilizing full 16-bit range (0-65535) afforded by the detector. Limits for all three colors are picked and data is rescaled. This results in apparent brightening of the image.



In [10]:

    
fig = plt.figure(figsize=(10, 7))
fig.set_facecolor('white')

for color, channel in zip('rgb', np.rollaxis(img432, axis=-1)):
    counts, centers = exposure.histogram(channel)
    plt.plot(centers[1::], counts[1::], color=color)

plt.show()



In [11]:

    
img432_ha = np.empty(img432.shape, dtype='uint16')
lims = [(7100,14500), (8200, 14000), (9200,13500)]
for lim, channel in zip(lims, range(3)):
    img432_ha[:, :, channel] = exposure.rescale_intensity(img432[:, :, channel], lim)

color_image_show(img432_ha, '4-3-2 image, histogram equilized')

Contrast and Color Balance Adjustment

Adjusting contrast and color balance of an image to their "right" levels is an iterative process. It involves the knowledge of how images at the same latitude and longitude look at a particular time of year. Users can draw on previous experience or search for images of the similar scenes on the Internet to get appropriate contrast and color balance. In this case, we will adjust balance of green and blue colors by shifting them to partially "brighter" values using gamma adjustment.



In [12]:

    
img432_ha[:, :, 1] = exposure.adjust_gamma(img432_ha[:, :, 1], 0.65)
img432_ha[:, :, 2] = exposure.adjust_gamma(img432_ha[:, :, 2], 0.75)
color_image_show(img432_ha, '4-3-2 image, histogram equilized, color gamma adjusted')

Comparison to European Space Agency Image

For comparison, here is image processed by the professionals at ESA, showing the same satellite scene.



In [13]:

    
Image('http://www.esa.int/var/esa/storage/images/esa_multimedia/images/2014/01/'
      'guinea-bissau_and_the_bissagos_islands/13474677-1-eng-GB/'
      'Guinea-Bissau_and_the_Bissagos_islands_node_full_image_2.jpg')









    Out[13]: