landsat 8 python library

read the true color

def landsat8_rgb_show(file_path):
    """
    show the true color of landsat 8
    file_path: the path of the true color file
    """
    file_name = os.path.join(file_path, os.path.basename(file_path) + ".tif")
    raster_rgb_array = plt.imread(file_name)
    plt.imshow(raster_rgb_array)
    plt.show()

calculate the metadata of Lansat-8:

def landsat8_read_meta(file_path, show_flag=False):
    file_name = os.path.basename(file_path)[:42] + "_MTL" ".txt"
    file_name = os.path.join(file_path, os.path.basename(file_path), file_name)
    i = 0
    metadata = {}
    with open(file_name, "r") as mtl_obj:
        data = mtl_obj.readline()
        cur_first_key = data.split(" = ")[1].split("\n")[0]
        metadata[cur_first_key] = {}
        while data[:9] != "END_GROUP":
            i = i + 1
            data = mtl_obj.readline()
            cur_sec_key = data.split(" = ")[1].split("\n")[0]
            metadata[cur_first_key][cur_sec_key] = {}
            data = mtl_obj.readline()
            if data == "END\n":
                break
            while data[:11] != "  END_GROUP":
                meta_key = data.split("    ")[1].split(" = ")
                if meta_key[1].find("\"") != -1:
                    metadata[cur_first_key][cur_sec_key][meta_key[0]] = meta_key[1].split("\"")[1]
                else:
                    if meta_key[0] == "FILE_DATE":
                        metadata[cur_first_key][cur_sec_key][meta_key[0]] = meta_key[1].split("\n")[0]
                    else:
                        metadata[cur_first_key][cur_sec_key][meta_key[0]] = meta_key[1].split("\n")[0]
                data = mtl_obj.readline()
                i = i + 1
                if i > 3000:
                    return
    if show_flag:
        for key in metadata.keys():
            metadata_sec = metadata[key]
            for key_sec in metadata_sec.keys():
                print("\t" + key_sec)
                metadata_third = metadata_sec[key_sec]
                for key_third in metadata_third.keys():
                    print("\t\t" + key_third + "\t:\t" + metadata_third[key_third])
    return metadata

Download the srtm dem data from the earthdata.nasa website, conduct the topographic correction.
image mosaic with arcgis
Arcmap->Data Management Tools->Raster->Mosaic Dataset->Mosaic To New Raster

Calulate the slope from a dem:
method 1 (using the gdal.DEMProcessing):

from osgeo import gdal
import numpy as np
import rasterio.io
import matplotlib.pyplot as plt
import os

def slope_and_aspect_cal(file_path, show_flag=False):  # uncompleted

    # read the dem
    dem_path = os.path.join(file_path, "dem_RimFire.tif")
    dem_raster = gdal.Open(dem_path)
    dem_array = dem_raster.GetRasterBand(1).ReadAsArray()  # get the dem raster

    def calculate_slope(DEM):
        slope_path = os.path.join(file_path, "slope.tif")
        gdal.DEMProcessing(slope_path, DEM, 'slope')
        with rasterio.open(slope_path) as dataset:
            slope = dataset.read(1)
        return slope

    def calculate_aspect(DEM):
        slope_path = os.path.join(file_path, "slope.tif")
        gdal.DEMProcessing(slope_path, DEM, 'aspect')
        with rasterio.open(slope_path) as dataset:
            aspect = dataset.read(1)
        return aspect

    slope = calculate_slope(dem_path)
    aspect = calculate_aspect(dem_path)

    print(type(slope), slope.dtype, slope.shape)
    print(slope.min(), slope.max())
    print("The number of pixels less than 0: %d" % len(slope < 0))

    plt.hist(slope.flatten(), bins=256)
    plt.title(r'Histogram of slope')
    plt.show()

    slope[slope <0] = 255

    if show_flag is True:
        plt.imshow(slope)
        plt.title("slope")
        plt.show()
    return slope

method 2 (using the richdem.

import os
from osgeo import gdal

def topo_correction_cal(file_path, show_flag=False):  # uncompleted
    # read the dem
    dem_path = os.path.join(file_path, "dem_RimFire.tif")
    dem_raster = gdal.Open(dem_path)
    dem_array = dem_raster.GetRasterBand(1).ReadAsArray()  # get the dem raster
    import richdem as rd
    shasta_dem = rd.LoadGDAL(dem_path, no_data=0)
    slope = rd.TerrainAttribute(shasta_dem, attrib='slope_riserun')
    rd.rdShow(slope, axes=False, cmap='magma', figsize=(8, 5.5))
    plt.show()
    return slope

atmospheric correction with ARCSI

In the view of the atmospheric correction, ARCSI is a interface using the python.