modis_owl.py

import numpy as np
from osgeo import gdal
from PIL import Image
import netCDF4
import json
import os
import sys
import datetime
import argparse
import time

def get_path(col_name, year, julday, h, v):
    d = datetime.datetime(year, 1, 1, 0, 0, 0, 0)
    d = d + datetime.timedelta(days=julday-1)

    data_path = os.path.join(root_path, col_name + ".006", "{}.{:02d}.{:02d}".format(d.year, d.month, d.day))
    tile_path = [os.path.join(data_path, f) for f in os.listdir(data_path) if f.startswith("{}.A{}{:03d}.h{:02d}v{:02d}.006.".format(col_name, d.year, julday, tile_h, tile_v)) and f.endswith(".hdf")]

    if len(tile_path) != 1:
        sys.exit("File not found")

    return tile_path[0]

def compute_owl(tile_path):

    # Get [RED, NIR, SWIR2, SWIR3]
    bands = [gdal.Open('HDF4_EOS:EOS_GRID:"{}":MOD_Grid_BRDF:Nadir_Reflectance_Band{}'.format(tile_path, b)) for b in [1, 2, 6, 7]]

    b_stack = np.empty((2400*2400, 6), dtype=np.float32)

    # Add bands
    for i, b in enumerate(bands):
        b_stack[:, i] = np.nan_to_num(b.GetRasterBand(1).ReadAsArray().reshape((2400*2400, ))*np.float32(.0001))


    # Compute and add NDVI = (NIR-RED)/(NIR+RED)
    b_stack[:, 4] = np.nan_to_num((b_stack[:, 1] - b_stack[:, 0]) / (b_stack[:, 1] + b_stack[:, 0]))

    # Compute and add NDWI = (NIR-SWIR2)/(NIR+SWIR2)
    b_stack[:, 5] = np.nan_to_num((b_stack[:, 1] - b_stack[:, 2]) / (b_stack[:, 1] + b_stack[:, 2]))

    G0 = np.float32(1.00)
    G1 = np.float32(-3.4137561)
    G2 = np.float32(-0.000959735270)
    G3 = np.float32(10.00417955330)
    G4 = np.float32(114.1927990)
    G5 = np.float32(-0.430407140)
    G6 = np.float32(-0.0961932990)
    MrVBF = np.float32(110)

    z = G1 + G2 * b_stack[:, 2] + G3 * b_stack[:, 3] + G4 * b_stack[:, 4] + G5 * b_stack[:, 5] + G6 * MrVBF
    owl = G0 / (1 + np.exp(z))

    return owl.reshape((2400, 2400))


def generate_thumbnail(arr, file_name):
    norm_arr = (arr.clip(0,1.0)*255).astype(np.uint8)
    img = Image.fromarray(norm_arr)
    img.save(file_name + ".png")


if __name__ == "__main__":

    # python owl.py /g/data2/u39/public/data/modis/lpdaac-tiles-c6/ 2011 1 4 31 11 ./
    parser = argparse.ArgumentParser(description="""Modis OWL Analysis argument parser""")
    parser.add_argument(dest="root_path", type=str, help="Root path")
    parser.add_argument(dest="year", type=int, help="Year of a modis tile (HDF-EOS).")
    parser.add_argument(dest="julday", type=int, help="Julian Day of a modis tile (HDF-EOS).")
    parser.add_argument(dest="h", type=int, help="Horizontal Corrrdinate of Modis Tile (HDF-EOS).")
    parser.add_argument(dest="v", type=int, help="Vertical Corrrdinate of Modis Tile (HDF-EOS).")
    parser.add_argument(dest="destination", type=str, help="Full path to destination.")
    args = parser.parse_args()

    root_path = args.root_path
    tile_year = args.year
    tile_julday = args.julday
    tile_h = args.h
    tile_v = args.v
    dest_path = args.destination

    mcd43a4 = get_path("MCD43A4", tile_year, tile_julday, tile_h, tile_v)
    print(mcd43a4)

    stack = compute_owl(mcd43a4)
    print(stack.shape)
    generate_thumbnail(stack, "owl_{}{:03d}_h{}v{}".format(tile_year, tile_julday, tile_h, tile_v))

    #mcd43a2 = get_path("MCD43A2", tile_year, tile_julday, tile_h, tile_v)
    #print(mcd43a2)