# from os.path import expanduser
import pandas as pd, numpy as np,os
import matplotlib.pyplot as plt
import glob
import csv,gdal,tarfile
from scipy.ndimage import gaussian_filter

# home = expanduser("~")
home = r'Z:\townsenduser-rw\HyspexPro\Output\Cheesehead_V2\CHEESEHEAD_20190830\BRDF\Tiff'
images = glob.glob("%s/CHEESE*_rgbim.tif" % home)
images.sort()
for image in images:
    img = gdal.Open(image)
    green = img.GetRasterBand(2).ReadAsArray()
    # green = img.GetRasterBand(49).ReadAsArray() # for hyspex img, use band 49: 558.98nm as green band

    V=green.copy()
    V[green <0]=0
    VV=gaussian_filter(V,sigma=.5,truncate=4)
    
    W=0*green.copy()+1
    W[green<0]=0
    WW=gaussian_filter(W,sigma=.5,truncate=4)

    # Generate np.nan for the surrounding background
    Z=VV/WW
    
    green[green<0]= Z[green<0]
    green[np.isnan(green)] = -1
    
    # output_name = "/home/adam/Dropbox/rs/cheese/data/%s" % os.path.basename(image)
    # output_name =output_name.replace('rgbim','g')
    output_name = image.replace('_rgbim', '_g')

    
    datatype = gdal.GDT_Int16
    
    # Set the output raster transform and projection properties
    driver = gdal.GetDriverByName("GTIFF")
    tiff = driver.Create(output_name,img.RasterXSize,img.RasterYSize,1,datatype)
    tiff.SetGeoTransform(img.GetGeoTransform())
    tiff.SetProjection(img.GetProjection())
        

    tiff.GetRasterBand(1).WriteArray(green)
    tiff.GetRasterBand(1).SetNoDataValue(-1)

    del tiff, driver        

# Extract green band from naip image. (skip if have green band file ready)
# img = gdal.Open("%s/Downloads/chs_naip_2018_1m.tif" % home)
# green = img.GetRasterBand(2).ReadAsArray().astype(int)
#
# output_name = "/home/adam/Dropbox/rs/cheese/data/chs_naip_2018_1m_g.tif"
# output_name =output_name.replace('rgbim','g')
#
# datatype = gdal.GDT_Int16
#
# # Set the output raster transform and projection properties
# driver = gdal.GetDriverByName("GTIFF")
# tiff = driver.Create(output_name,img.RasterXSize,img.RasterYSize,1,datatype)
# tiff.SetGeoTransform(img.GetGeoTransform())
# tiff.SetProjection(img.GetProjection())
#
#
# tiff.GetRasterBand(1).WriteArray(green)
#
# del tiff, driver