modis_edison.py

import numpy as np
from osgeo import gdal
import scipy.optimize as opt
import scipy.ndimage
from PIL import Image
import netCDF4
import json
import os
import sys
import datetime
import argparse

def get_nc_name(hdf_tile_path, dest, suffix=""):
    file_name = hdf_tile_path.split("/")[-1][:-4]
    file_name_parts = file_name.split(".")
    folder =  ".".join([file_name_parts[2], file_name_parts[1][1:-3]])
    dest_fold = os.path.join(dest, folder)
    if not os.path.exists(dest_fold):
        os.makedirs(dest_fold)
    short_file_name = file_name.split(".")[:-1]
    return os.path.join(dest_fold, "FC" + suffix + ".v302." + ".".join(short_file_name) + ".nc")


def pack_data(hdf_file, arr, dest):

    file_name = get_nc_name(hdf_file, dest)

    with netCDF4.Dataset(file_name, 'w', format='NETCDF4') as dest:
        with open('nc_metadata.json') as data_file:
            attrs = json.load(data_file)
            for key in attrs:
                setattr(dest, key, attrs[key])
        setattr(dest, "date_created", datetime.datetime.now().strftime("%Y%m%dT%H%M%S"))

        ds = gdal.Open('HDF4_EOS:EOS_GRID:"{}":MOD_Grid_BRDF:Nadir_Reflectance_Band1'.format(hdf_file))
        proj_wkt = ds.GetProjection()
        geot = ds.GetGeoTransform()

        x_dim = dest.createDimension("x", ds.RasterXSize)
        y_dim = dest.createDimension("y", ds.RasterYSize)
        b_dim = dest.createDimension("b", 3)

        var = dest.createVariable("x", "f8", ("x",))
        var.units = "m"
        var.long_name = "x coordinate of projection"
        var.standard_name = "projection_x_coordinate"
        var[:] = np.linspace(geot[0], geot[0]+(geot[1]*ds.RasterXSize), ds.RasterXSize)

        var = dest.createVariable("y", "f8", ("y",))
        var.units = "m"
        var.long_name = "y coordinate of projection"
        var.standard_name = "projection_y_coordinate"
        var[:] = np.linspace(geot[3], geot[3]+(geot[5]*ds.RasterYSize), ds.RasterYSize)

        var = dest.createVariable("phot_veg", "u1", ("y", "x"), fill_value=255)
        var.long_name = "Photosynthetic Vegetation"
        var.units = '%'
        var.grid_mapping = "sinusoidal"
        var[:] = arr[:, :, 0]

        var = dest.createVariable("nphot_veg", "u1", ("y", "x"), fill_value=255)
        var.long_name = "Non Photosynthetic Vegetation"
        var.units = '%'
        var.grid_mapping = "sinusoidal"
        var[:] = arr[:, :, 1]

        var = dest.createVariable("bare_soil", "u1", ("y", "x"), fill_value=255)
        var.long_name = "Bare Soil"
        var.units = '%'
        var.grid_mapping = "sinusoidal"
        var[:] = arr[:, :, 2]

        generate_thumbnail(arr, file_name)

        var = dest.createVariable("sinusoidal", 'S1', ())

        var.grid_mapping_name = "sinusoidal"
        var.false_easting = 0.0
        var.false_northing = 0.0
        var.longitude_of_central_meridian = 0.0
        var.longitude_of_prime_meridian = 0.0
        var.semi_major_axis = 6371007.181
        var.inverse_flattening = 0.0
        var.spatial_ref = proj_wkt
        var.GeoTransform = "{} {} {} {} {} {} ".format(*[geot[i] for i in range(6)])


def jp_superfunc(A, arr, mask):
    res = np.zeros((arr.shape[0], A.shape[1]), dtype=np.uint8)
    for i in range(arr.shape[0]):
        #if mask[i]:
        if mask[i] and (arr[i, :7] < 0).sum() == 0 and (arr[i, :7] > 1).sum() == 0:
            res[i, :] = (opt.nnls(A, arr[i, :])[0].clip(0, 2.54)*100).astype(np.uint8)
        else:
            res[i, :] = np.ones((A.shape[1]), dtype=np.uint8)*(255)
    return res

def input_mask(pq_mask_path):
    pq_ds = gdal.Open('HDF4_EOS:EOS_GRID:"{}":MOD_Grid_BRDF:BRDF_Albedo_Ancillary'.format(pq_mask_path))
    pq_raw = pq_ds.GetRasterBand(1).ReadAsArray()

    mask = 0b0000000000001111
    pq = pq_raw>>4 & mask

    snow_ds = gdal.Open('HDF4_EOS:EOS_GRID:"{}":MOD_Grid_BRDF:Snow_BRDF_Albedo'.format(pq_mask_path))
    snow_pq = np.equal(snow_ds.ReadAsArray(), 0)

    # True if pq is 1, 2 or 4
    pq = np.logical_and(snow_pq, np.logical_or(np.logical_or(np.equal(pq, np.ones(1)), np.equal(pq, np.ones(1)*2)), np.equal(pq, np.ones(1)*4)))
    #pq = np.logical_or(np.logical_or(np.equal(pq, np.ones(1)), np.equal(pq, np.ones(1)*2)), np.equal(pq, np.ones(1)*4))

    return pq.reshape(2400*2400)

def input_stack(tile_path):

    bands = [gdal.Open('HDF4_EOS:EOS_GRID:"{}":MOD_Grid_BRDF:Nadir_Reflectance_Band{}'.format(tile_path, b)) for b in range(1, 8)]

    #add band
    bands_stack = (bands[0].GetRasterBand(1).ReadAsArray()*.0001).astype(np.float32)
    for b in range(1, len(bands)):
        bands_stack = np.dstack((bands_stack, (bands[b].GetRasterBand(1).ReadAsArray()*.0001).astype(np.float32)))

    bands_stack = bands_stack.reshape(2400*2400, 7)

    #add log(band)
    for b in range(7):
        bands_stack = np.hstack((bands_stack, np.expand_dims(np.log(bands_stack[:, b]), axis=1)))

    #add band*log(band)
    for b in range(7):
        bands_stack = np.hstack((bands_stack, np.expand_dims(np.multiply(bands_stack[:, b], bands_stack[:, b+7]), axis=1)))

    #add band*next_band
    for b in range(7):
        for b2 in range(b+1, 7):
            bands_stack = np.hstack((bands_stack, np.expand_dims(np.multiply(bands_stack[:, b], bands_stack[:, b2]), axis=1)))

    #add log(band)*log(next_band)
    for b in range(7):
        for b2 in range(b+1, 7):
            bands_stack = np.hstack((bands_stack, np.expand_dims(np.multiply(bands_stack[:, b+7], bands_stack[:, b2+7]), axis=1)))

    #add (next_band-band)/(next_band+band)
    for b in range(7):
        for b2 in range(b+1, 7):
            bands_stack = np.hstack((bands_stack, np.expand_dims(np.divide(bands_stack[:, b2]-bands_stack[:, b], bands_stack[:, b2]+bands_stack[:, b]), axis=1)))

    bands_stack = np.nan_to_num(bands_stack)

    #--------------------------------
    # this bit new, by JPG
    ones = np.ones(bands_stack.shape[0])
    ones = ones.reshape(ones.shape[0], 1)
    bands_stack = np.concatenate((bands_stack, ones), axis=1)
    #--------------------------------

    return bands_stack


def members():

    A = np.load("endmembers.npy")
    #--------------------------------
    # this bit new, by JPG
    SumToOneWeight = 0.02
    ones = np.ones(A.shape[1]) * SumToOneWeight
    ones = ones.reshape(1, A.shape[1])
    A = np.concatenate((A, ones), axis=0).astype(np.float32)

    return A


if __name__ == "__main__":

    root_path = "/g/data2/u39/public/data/modis/lpdaac-tiles-c5/"

    parser = argparse.ArgumentParser(description="""Modis Vegetation Analysis argument parser""")
    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()

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

    tile_date = datetime.datetime(tile_year, 1, 1) + datetime.timedelta(days=tile_julday-1)

    tile_path = os.path.join(root_path, "MCD43A4.005", "{0:04d}.{1:02d}.{2:02d}".format(tile_date.year, tile_date.month, tile_date.day))
    tile = [tile for tile in os.listdir(tile_path) if (tile.endswith(".hdf") and tile.startswith("MCD43A4.A{0:04d}{1:03d}.h{2:02d}v{3:02d}.".format(tile_year, tile_julday, tile_h, tile_v)))]
    if len(tile) != 1:
        print("Tile not found", tile_year, tile_julday, tile_h, tile_v)
        sys.exit(1)

    modis_tile = os.path.join(tile_path, tile[0])

    file_name = get_nc_name(modis_tile, dest_path)

    if os.path.isfile(file_name):
        print("Tile {} already exists".format(file_name))
        sys.exit(0)

    bands_stack = input_stack(modis_tile)

    mask_path = os.path.join(root_path, "MCD43A2.005", "{0:04d}.{1:02d}.{2:02d}".format(tile_date.year, tile_date.month, tile_date.day))
    mask = [tile for tile in os.listdir(mask_path) if (tile.endswith(".hdf") and tile.startswith("MCD43A2.A{0:04d}{1:03d}.h{2:02d}v{3:02d}.".format(tile_year, tile_julday, tile_h, tile_v)))]
    if len(mask) != 1:
        print("Mask not found", tile_year, tile_julday, tile_h, tile_v)
        sys.exit(1)

    mask_tile = os.path.join(mask_path, mask[0])
    mask = input_mask(mask_tile)

    A = members()

    out = jp_superfunc(A, bands_stack, mask).reshape(2400, 2400, 3)
    pack_data(modis_tile, out, dest_path)
    pack_data_lowres(modis_tile, out, dest_path)