# 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

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

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

b1 = read_band(1)

image_adjust_brightness(b1, 10000, 13500, 'Band 1, data set LC82040522013123LGN01')

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

b5 = read_band(5)

image_adjust_brightness(b5, 6000, 23000, 'Band 5, data set LC82040522013123LGN01')

cd /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')

cd /Users/kronos/gis/l8/paris/

b9 = read_band(9)

image_adjust_brightness(b9, 5000, 7500, 'Band 9, data set LC81990262013104LGN01')

cd /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')