Blogs

DK for Python Supporting NumPy

This example features the Developer Kernel for Python edition with Jupyter Notebook being used with the NumPy (scientific computing) and Matplotlib (data visualization & graphical plotting) packages to perform satellite/aerial imagery analyses.

Import required modules¶

In [1]:

import tatukgis_pdk as pdk
import glob
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from numbers import Number
np.seterr(divide='ignore', invalid='ignore');

In [2]:

band_desc = {
    "B1": "Coastal / Aerosol",
    "B2": "Visible blue",
    "B3": "Visible green",
    "B4": "Visible red",
    "B5": "Near-infrared",
    "B6": "Short wavelength infrared",
    "B7": "Short wavelength infrared",
    "B8": "Panchromatic",
    "B9": "Cirrus",
    "B10": "Long wavelength infrared",
    "B11": "Long wavelength infrared"
}

bT+ibClfzI5TI8jWZQ4OaU9YdiXs+oipKGT/qrwk3/NJ5jUMRvcSk6/mnUk6hZlP4ESfT0z+7l6w=")

Add helper functions for loading raster data to numpy array and normalizing matrices¶

In [4]:

def band_to_array(band_path):
    """Reads a raster from path and returns
    a band name and an array with the band data.

    Args:
        band_path (str): path to a raster image

    Returns:
        str, numpy.ndarray: band name, band data as numpy array
    """
    # get band name from long path
    # "LC08_L1TP_014032_20191024_20191030_01_T1_B1.TIF" 
    band_name = Path(band_path).stem.split("_")[7]
    
    # create a pixel layer
    lp = pdk.TGIS_Utils.GisCreateLayer(band_name, band_path)
    lp.Open()
    
    # interpretate a pixel lauyer as a grid
    lp.Interpretation = pdk.TGIS_LayerPixelInterpretation().Grid
    lp.Params.Pixel.GridBand = 1
    
    # lock a whole grid for reading
    lp_lock = lp.LockPixels(lp.Extent, lp.CS, False)
    try:
        # get grid cells as a numpy array
        band_arr = lp_lock.AsArray()
    finally:
        # remember to unlock!
        lp.UnlockPixels(lp_lock)
    
    return band_name, band_arr


def norm(arr, scale=2**16, factor=1.0):
    """Normalizes arr in range 0..1.

    Args:
        arr (numpy.ndarray): numpy array representation
                             of the grid or the pixel layer
        scale (Number)     : divide an arr by the scale (divisor)
        factor (Number)    : multiply arr by a factor (multiplicand)

    Returns:
        numpy.ndarray: normalized array
    """
    res_arr = (arr / scale * factor)
    res_arr[res_arr>1] = 1
    return res_arr

Load Landsat 8 data¶

In [5]:

landsat_paths = glob.glob("E:\\GisData\\Landsat 8\\Cincinnati\\cut\\*.tif")

# list paths
print(*landsat_paths, sep="\n")

# store all landsat bands in a dictionary with band name as a key
bands = {name: arr for (name, arr) in map(band_to_array, landsat_paths)}

E:\GisData\Landsat 8\Cincinnati\cut\LC08_L1TP_020033_20220604_20220610_02_T1_B1.tif
E:\GisData\Landsat 8\Cincinnati\cut\LC08_L1TP_020033_20220604_20220610_02_T1_B10.tif
E:\GisData\Landsat 8\Cincinnati\cut\LC08_L1TP_020033_20220604_20220610_02_T1_B11.tif
E:\GisData\Landsat 8\Cincinnati\cut\LC08_L1TP_020033_20220604_20220610_02_T1_B2.tif
E:\GisData\Landsat 8\Cincinnati\cut\LC08_L1TP_020033_20220604_20220610_02_T1_B3.tif
E:\GisData\Landsat 8\Cincinnati\cut\LC08_L1TP_020033_20220604_20220610_02_T1_B4.tif
E:\GisData\Landsat 8\Cincinnati\cut\LC08_L1TP_020033_20220604_20220610_02_T1_B5.tif
E:\GisData\Landsat 8\Cincinnati\cut\LC08_L1TP_020033_20220604_20220610_02_T1_B6.tif
E:\GisData\Landsat 8\Cincinnati\cut\LC08_L1TP_020033_20220604_20220610_02_T1_B7.tif
E:\GisData\Landsat 8\Cincinnati\cut\LC08_L1TP_020033_20220604_20220610_02_T1_B8.tif
E:\GisData\Landsat 8\Cincinnati\cut\LC08_L1TP_020033_20220604_20220610_02_T1_B9.tif

Vizualize all bands¶

In [6]:

fig, axes = plt.subplots(4, 3,
                         subplot_kw={'xticks': [], 'yticks': []},
                         figsize=(10,15))
for i in range(len(bands)):
    band_name = f"B{i+1}"
    band_title = f"{band_name}\n{band_desc[band_name]}"
    if not band_name in bands:
        continue
    band_arr = bands[band_name]
    
    axes[i // 3, i % 3].imshow(band_arr, cmap="gist_earth")
    axes[i // 3, i % 3].set_title(band_title)

axes[3, 2].remove()  # don't display empty axis

fig.tight_layout()
plt.show()

Create and display RGB Composite Image¶

In [7]:

plt.rcParams['figure.figsize'] = [10, 10]

rgb = norm(np.stack((bands["B4"], bands["B3"], bands["B2"]), axis=-1),
           factor=2)
plt.imshow(rgb)
plt.axis('off')
plt.title("True Color Composite Image (RGB, Landsat 8 bands: 4, 3, 2")
plt.show()

Create and display False Composite Image¶

In [8]:

false_composite = norm(np.stack((bands["B5"], bands["B4"], bands["B3"]), axis=-1),
                       factor=1.5)
plt.imshow(false_composite)
plt.axis('off')
plt.title("False Color Composite Image (Landsat 8 bands: 5, 4, 3)")
plt.show()

Calculate and display Normalized Difference Vegetation Index (NDVI)¶

In [9]:

ndvi = (bands["B5"] - bands["B4"]) / (bands["B5"] + bands["B4"])
plt.imshow(ndvi, cmap="Greens")
plt.colorbar()
plt.axis('off')
plt.title("Normalized Difference Vegetation Index (NDVI)")
plt.show()

Calculate and display Modified Normalized Difference Water Index (MNDWI)¶

In [10]:

mndwi = (bands["B3"] - bands["B6"]) / (bands["B3"] + bands["B6"])
plt.imshow(mndwi, cmap="Blues")
plt.colorbar()
plt.axis('off')
plt.title("Modified Normalized Difference Water Index (MNDWI)")
plt.show()

Extract water pixels¶

In [11]:

water_mask = mndwi > 0
rgb_water = rgb.copy()
rgb_water[water_mask] = np.array([147, 181, 234]) / 255.0
plt.imshow(rgb_water)
plt.axis('off')
plt.title("Water pixels overlayed on the true color image")
plt.show()

Extract Built-up areas¶

In [12]:

ndbi = (bands["B6"] - bands["B5"]) / (bands["B6"] + bands["B5"])
bu = ndbi - ndvi
builtup_mask = bu > -0.5
bu[builtup_mask] = 1
bu[~builtup_mask] = 0
bu[water_mask] = 0

plt.imshow(bu, cmap="Reds", interpolation="nearest")
plt.colorbar()
plt.axis('off')
plt.title("Built-up areas (BU) extracted from Landsat 8 images")
plt.show()

Posted: September 26, 2022
Filed under: ANALYSES, PYTHON