

import time
import os
import datetime as dt
import numpy as np
import gdal

from scipy import stats
import glob, tarfile
import psutil

from joblib import Parallel, delayed


from multiprocessing import Process, Pipe

import argparse

import gc
import utility_functions as utility_functions
import trends_calculation as trends_calculation

#for testing
#start_year = 1988
#end_year = 1989

#start_year = 1984
#end_year = 1992

#start_year = 1984
#end_year = 2018
#start_year = 1984
#end_year = 2019
#start_year = 1984
#end_year = 2018

#############################From: lcmap.py############################################
qamap = {'fill': 0,
         'clear': 1,
         'water': 2,
         'shadow': 3,
         'snow': 4,
         'cloud': 5,
         'cloud_conf1': 6,
         'cloud_conf2': 7,
         'cirrus1': 8,
         'cirrus2': 9,
         'occulsion': 10}

#############################From: Kelcy - SR-stuff.ipynb############################################
def bitmask(arr, bit):
    """
    Create a bool mask where the bit applies to the array
    """
    #reza explanation
    #take a 1 which is in bit position 0 and shift it left to the bit index with <<
    #then do an & to get the matching bits
    #then see if we have greater than 0 matching bits
    return arr & 1 << bit > 0

##############################################################################
def get_LSTclear_idx(LST, st_qa, px_qa):
    '''
    Context on scale values from ard dfcb
    If you want to un-scale the ST data that you read from ARD, you would have to multiply by 0.1
    similarly, if you want to un-scale the ST_QA read from ARD, you would have to multiply by 0.01
    However, in this example, we are not reading from ARD. We are adding a hardcoded threshold which
    is in units of Kelvin. The un-scaled value is 5 Kelvin. If you want to scale it using the ST_QA
    scale factor, then we must divide it by the scale factor (/0.01 or multiply it by 100)

    we could also just unscale both values if we wanted to
    '''
    five_kelvin_scaled = (5/0.01)
    return np.asarray((st_qa > five_kelvin_scaled)*(st_qa >= 0)\
           * (bitmask(px_qa, qamap['water']) + bitmask(px_qa, qamap['clear']))\
           * (LST > 2731.5))

def get_Thermalclear_idx(px_qa):
    return bitmask(px_qa, qamap['water']) + bitmask(px_qa, qamap['clear'])



def extract_ST(workpath, HV):
    inputdir = os.path.join(workpath, HV)
    outDir = os.path.join(workpath, HV, 'extract')

    if not os.path.exists(outDir):  ## if outfolder is not already existed creating one
        os.makedirs(outDir,exist_ok=True)
    print(inputdir)
    tarFiles = glob.glob(inputdir + os.sep + "*.tar")

    #--------------------------------------------------multi threaded version

    def extract(file_path):
        print(file_path)
        baz = tarfile.open(file_path, 'r:')
        for item in baz:
            # print(item)
            if '_ST.tif' in item.name:
                baz.extract(item, path=outDir)
            if '_STQA.tif' in item.name:
                baz.extract(item, path=outDir)
            if '_PIXELQA.tif' in item.name:
                baz.extract(item, path=outDir)
        baz.close()

    file_paths = []
    for count, tarF in enumerate(tarFiles):
        file_paths.append(tarF)

    num_cpu = psutil.cpu_count(logical=True)
    Parallel(n_jobs=num_cpu)(delayed(extract)(some_file_path) for some_file_path in file_paths)

    #--------------------------------------------------single threaded version
    '''
    for count, tarF in enumerate(tarFiles):
        print(tarF)
        baz = tarfile.open(tarF, 'r:')
        for item in baz:
            # print(item)
            if '_ST.tif' in item.name:
                baz.extract(item, path=outDir)
            if '_STQA.tif' in item.name:
                baz.extract(item, path=outDir)
            if '_PIXELQA.tif' in item.name:
                baz.extract(item, path=outDir)
        baz.close()
    '''
    #--------------------------------------------------



def ST_annual_stats(workpath, HV,start_year_inclusive,end_year_inclusive):

    start_year = start_year_inclusive
    end_year = end_year_inclusive

    inputdir = os.path.join(workpath, HV, 'extract')
    outDir = os.path.join(workpath, HV, 'stats')

    if not os.path.exists(outDir):  ## if outfolder is not already existed creating one
        os.makedirs(outDir,exist_ok=True)

    st_files = np.asarray(glob.glob(inputdir + os.sep + "*_ST.tif"))
    trans, prj = utility_functions.get_geo(st_files[0])
    dates_list = []
    for st_file in st_files:
        dates_list.append(dt.datetime.strptime(os.path.basename(st_file)[15:23], "%Y%m%d"))

    for iyear in range(start_year, end_year + 1):
        year_idx = np.asarray([iyear == dat.year for dat in dates_list])
        if sum(year_idx) == 0:
            continue
        
        
        print("total files for year "+str(iyear)+" "+str(len(st_files[year_idx])))
        data = None
        #get all the files that have that year
        for count, st_file in enumerate(st_files[year_idx]):
            print(iyear, st_file)
            # dates.append(dt.datetime.strptime(os.path.basename(st_file)[15:23], "%Y%m%d"))
            stQA_file = st_file.replace('_ST.tif', '_STQA.tif')
            pxQA_file = st_file.replace('_ST.tif', '_PIXELQA.tif')

            st_data = utility_functions.get_tifflayer(st_file)
            stQA_data = utility_functions.get_tifflayer(stQA_file,dtype=np.int)
            pxQA_data = utility_functions.get_tifflayer(pxQA_file,dtype=np.uint16)

            mask = get_LSTclear_idx(st_data, stQA_data, pxQA_data)
            st_data[~mask] = np.nan

            if data is None:
                data = np.zeros((st_data.shape[0], st_data.shape[1], sum(year_idx)))

            # print('dims: ', count, st_data.shape[0], st_data.shape[1], sum(year_idx), st_data.shape, data.shape)
            data[:, :, count] = st_data

        y_res, x_res, t_res = data.shape
        if 'data_annual' not in locals():
            data_annual = np.zeros((st_data.shape[0], st_data.shape[1], (end_year - start_year + 1)))

        ##########################Annual mean ST#########################################################

        #axis goes 0,1,2, so we are doing these stats per pixel through the time dimension
        # 10 - 273.15 is for converting the land surface temperature celsius
        data_annual[:, :, iyear - start_year] = np.nanmean(data, axis=2) / 10 - 273.15
        utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MeanT_{}.tif'.format(HV, iyear)),
                    data_annual[:, :, iyear - start_year],
                    [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        #########################Annual max ST#####################################################
        data_annual[:, :, iyear - start_year] = np.nanmax(data, axis=2) / 10 - 273.15
        utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MaxT_{}.tif'.format(HV, iyear)),
                    data_annual[:, :, iyear - start_year],
                    [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        #########################Annual min ST#####################################################
        data_annual[:, :, iyear - start_year] = np.nanmin(data, axis=2) / 10 - 273.15
        utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MinT_{}.tif'.format(HV, iyear)),
                    data_annual[:, :, iyear - start_year],
                    [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        #########################Annual median ST#####################################################
        data_annual[:, :, iyear - start_year] = np.nanmedian(data, axis=2) / 10 - 273.15
        utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MedianT_{}.tif'.format(HV, iyear)),
                    data_annual[:, :, iyear - start_year],
                    [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)
                    
        gc.collect()

def ST_seasonal_stats(workpath, HV,start_year_inclusive,end_year_inclusive):

    start_year = start_year_inclusive
    end_year = end_year_inclusive

    inputdir = os.path.join(workpath, HV, 'extract')
    outDir = os.path.join(workpath, HV, 'stats')

    if not os.path.exists(outDir):  ## if outfolder is not already existed creating one
        os.makedirs(outDir,exist_ok=True)

    st_files = np.asarray(glob.glob(inputdir + os.sep + "*_ST.tif"))
    trans, prj = utility_functions.get_geo(st_files[0])
    dates_list = []
    for st_file in st_files:
        dates_list.append(dt.datetime.strptime(os.path.basename(st_file)[15:23], "%Y%m%d"))
    dates_list = np.asarray(dates_list)


    #--------------------------------------------------|


    def getPixelStackForMonthsFromYear(chosen_months,chosen_year):
        year_idx = np.asarray([iyear == dat.year for dat in dates_list])
        if sum(year_idx) == 0:
            return None
        
        data = None
        for count, st_file in enumerate(st_files[year_idx]):
            print(iyear, st_file)
            # dates.append(dt.datetime.strptime(os.path.basename(st_file)[15:23], "%Y%m%d"))
            stQA_file = st_file.replace('_ST.tif', '_STQA.tif')
            pxQA_file = st_file.replace('_ST.tif', '_PIXELQA.tif')

            st_data = utility_functions.get_tifflayer(st_file)
            stQA_data = utility_functions.get_tifflayer(stQA_file, dtype=np.int)
            pxQA_data = utility_functions.get_tifflayer(pxQA_file, dtype=np.uint16)

            mask = get_LSTclear_idx(st_data, stQA_data, pxQA_data)
            st_data[~mask] = np.nan

            if data is None:
                data = np.zeros((st_data.shape[0], st_data.shape[1], sum(year_idx)))
            # print('dims: ', count, st_data.shape[0], st_data.shape[1], sum(year_idx), st_data.shape, data.shape)
            data[:, :, count] = st_data

        y_res, x_res, t_res = data.shape
    
    
    previous_year = None
    previous_year_data = None


    ABC = None

    #--------------------------------------------------|




    for iyear in range(start_year, end_year + 1):
        year_idx = np.asarray([iyear == dat.year for dat in dates_list])
        if sum(year_idx) == 0:
            continue

        data = None
        for count, st_file in enumerate(st_files[year_idx]):
            print(iyear, st_file)
            # dates.append(dt.datetime.strptime(os.path.basename(st_file)[15:23], "%Y%m%d"))
            stQA_file = st_file.replace('_ST.tif', '_STQA.tif')
            pxQA_file = st_file.replace('_ST.tif', '_PIXELQA.tif')

            st_data = utility_functions.get_tifflayer(st_file)
            stQA_data = utility_functions.get_tifflayer(stQA_file,dtype=np.int)
            pxQA_data = utility_functions.get_tifflayer(pxQA_file,dtype=np.uint16)

            mask = get_LSTclear_idx(st_data, stQA_data, pxQA_data)
            st_data[~mask] = np.nan

            if data is None:
                data = np.zeros((st_data.shape[0], st_data.shape[1], sum(year_idx)))
            # print('dims: ', count, st_data.shape[0], st_data.shape[1], sum(year_idx), st_data.shape, data.shape)
            data[:, :, count] = st_data

        y_res, x_res, t_res = data.shape


        #----------------------------------------|

        '''
            since there is no data before 1984, have to do a mod
            so for 1984, (jan feb),(mar april may),(june july aug),(sep oct nov)
            but then 1985+ (dec jan feb),(mar april may),(june july aug),(sep oct nov)
        '''

        '''
        if previous_year is None:
            dec_jan_feb_season = [1, 2]
        else:
            dec_jan_feb_season = [1, 2, 3]

        mar_april_may_season = [3, 4, 5]
        june_july_aug_season = [6, 7, 8]
        sept_oct_nov_season = [9, 10, 11]
        '''

        season1 = np.asarray([1,2,3])
        season_idx = np.asarray([(i.month in season1) for i in dates_list[year_idx]])
        images_to_use = data[:, :, season_idx]

        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MaxSeason1_{}.tif'.format(HV, iyear)),
                        np.nanmax(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        #-----------------------------------------|


        ###############################Seasonal max################################################
        season1 = np.asarray([1,2,3])
        season_idx = np.asarray([(i.month in season1) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MaxSeason1_{}.tif'.format(HV, iyear)),
                        np.nanmax(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season2 = np.asarray([4,5,6])
        season_idx = np.asarray([(i.month in season2) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MaxSeason2_{}.tif'.format(HV, iyear)),
                        np.nanmax(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season3 = np.asarray([7,8,9])
        season_idx = np.asarray([(i.month in season3) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MaxSeason3_{}.tif'.format(HV, iyear)),
                        np.nanmax(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season4 = np.asarray([10,11,12])
        season_idx = np.asarray([(i.month in season4) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MaxSeason4_{}.tif'.format(HV, iyear)),
                        np.nanmax(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        # del data
        ###############################Seasonal min################################################
        season1 = np.asarray([1,2,3])
        season_idx = np.asarray([(i.month in season1) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MinSeason1_{}.tif'.format(HV, iyear)),
                        np.nanmin(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season2 = np.asarray([4,5,6])
        season_idx = np.asarray([(i.month in season2) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MinSeason2_{}.tif'.format(HV, iyear)),
                        np.nanmin(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season3 = np.asarray([7,8,9])
        season_idx = np.asarray([(i.month in season3) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MinSeason3_{}.tif'.format(HV, iyear)),
                        np.nanmin(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season4 = np.asarray([10,11,12])
        season_idx = np.asarray([(i.month in season4) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MinSeason4_{}.tif'.format(HV, iyear)),
                        np.nanmin(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        # del data
        #############################Seasonal mean###########################################
        season1 = np.asarray([1,2,3])
        season_idx = np.asarray([(i.month in season1) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MeanSeason1_{}.tif'.format(HV, iyear)),
                        np.nanmean(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season2 = np.asarray([4,5,6])
        season_idx = np.asarray([(i.month in season2) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MeanSeason2_{}.tif'.format(HV, iyear)),
                        np.nanmean(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season3 = np.asarray([7,8,9])
        season_idx = np.asarray([(i.month in season3) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MeanSeason3_{}.tif'.format(HV, iyear)),
                        np.nanmean(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season4 = np.asarray([10,11,12])
        season_idx = np.asarray([(i.month in season4) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MeanSeason4_{}.tif'.format(HV, iyear)),
                        np.nanmean(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        # del data
        #############################Seasonal median###########################################
        season1 = np.asarray([1,2,3])
        season_idx = np.asarray([(i.month in season1) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MedianSeason1_{}.tif'.format(HV, iyear)),
                        np.nanmedian(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season2 = np.asarray([4,5,6])
        season_idx = np.asarray([(i.month in season2) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MedianSeason2_{}.tif'.format(HV, iyear)),
                        np.nanmedian(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season3 = np.asarray([7,8,9])
        season_idx = np.asarray([(i.month in season3) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MedianSeason3_{}.tif'.format(HV, iyear)),
                        np.nanmedian(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        season4 = np.asarray([10,11,12])
        season_idx = np.asarray([(i.month in season4) for i in dates_list[year_idx]])
        if np.sum(season_idx) > 0:
            utility_functions.saveGeoTiff(os.path.join(outDir, '{}_MedianSeason4_{}.tif'.format(HV, iyear)),
                        np.nanmedian(data[:, :, season_idx], axis=2) / 10 - 273.15,
                        [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_Float32)

        del data




def num_clear(workpath, HV, start_year_inclusive, end_year_inclusive):

    start_year = start_year_inclusive
    end_year = end_year_inclusive

    inputdir = os.path.join(workpath, HV, 'extract')
    outDir = os.path.join(workpath, HV, 'stats')

    if not os.path.exists(outDir):  ## if outfolder is not already existed creating one
        os.makedirs(outDir,exist_ok=True)

    st_files = np.asarray(glob.glob(inputdir + os.sep + "*_ST.tif"))
    trans, prj = utility_functions.get_geo(st_files[0])
    dates_list = []
    for st_file in st_files:
        dates_list.append(dt.datetime.strptime(os.path.basename(st_file)[15:23], "%Y%m%d"))
    dates_list = np.asarray(dates_list)

    for iyear in range(start_year, end_year + 1):
        year_idx = np.asarray([iyear == dat.year for dat in dates_list])
        clear_img = np.zeros((5000, 5000), dtype=np.uint16)
        for st_file in st_files[year_idx]:
            stQA_file = st_file.replace('_ST.tif', '_STQA.tif')
            pxQA_file = st_file.replace('_ST.tif', '_PIXELQA.tif')

            st_data = utility_functions.get_tifflayer(st_file)
            stQA_data = utility_functions.get_tifflayer(stQA_file, dtype=np.int)
            pxQA_data = utility_functions.get_tifflayer(pxQA_file, dtype=np.uint16)

            clear_img += get_LSTclear_idx(st_data, stQA_data, pxQA_data)

        print('Saving', HV, iyear)
        y_res, x_res = clear_img.shape
        utility_functions.saveGeoTiff(os.path.join(outDir, '{}_numClear_{}.tif'.format(HV, iyear)),
                    clear_img, [x_res, y_res, 1], trans, prj, dtype=gdal.GDT_UInt16)


def determine_hv(x, y, aff=(-2565585, 150000, 0, 3314805, 0, -150000)):
    col = (x - aff[0] - aff[3] * aff[2]) / aff[1]
    row = (y - aff[3] - aff[0] * aff[4]) / aff[5]

    return int(col), int(row)



def st_stats_forHua(workpath,HV_List,start_year_inclusive,end_year_inclusive):

    # workpath = r'C:\QZhou\Therm_ARD'
    #workpath = r'S:\Temp\Shi\Impervious_surface\UHI\data\siouxfalls'
    #workpath = "C:\\Users\\rhussain\\Desktop\\LSRD\\data\\phoenix"
    #workpath = "C:\\Users\\rhussain\\Desktop\\LSRD\\data\\phoenix_dev"
    
    
    #official
    #workpath = "C:\\Users\\rhussain\\Desktop\\LSRD\\data\\atlanta"
    #HV_List = ["H23V13", "H23V14", "H24V13", "H24V14"]

    #test
    #workpath = "/caldera/projects/usgs/eros/urban_heat_islands/atlanta"
    #HV_List = ["H23V13"]
    
    #'H06V13' #['H06V13','H07V13']
    # workpath = r'/lcmap_data/bulk/Qiang/ST_ARD_mn'
    # HV = 'H16V06'

    start_year_inclusive = int(start_year_inclusive)
    end_year_inclusive = int(end_year_inclusive)

    for HV in HV_List:

        # unzips the tar files to the tile folder
        extract_ST(workpath, HV)

        #converts the files into {}_MeanT_{}.tif, MedianT, maxT, minT files
        #for each year

        ST_annual_stats(workpath, HV,start_year_inclusive, end_year_inclusive)

        #converts the files into MeanT,MedianT,maxT,minT for
        #these groups of months [1,2,3][4,5,6][7,8,9][10,11,12]
        ST_seasonal_stats(workpath, HV,start_year_inclusive, end_year_inclusive)

        print("starting trend_stats")

        start_time = time.time()
        trends_calculation.trend_stats(workpath, HV, 'MeanT', start_year_inclusive, end_year_inclusive)
        trends_calculation.trend_stats(workpath, HV, 'MedianT', start_year_inclusive, end_year_inclusive)
        trends_calculation.trend_stats(workpath, HV, 'MinT',start_year_inclusive, end_year_inclusive)
        trends_calculation.trend_stats(workpath, HV, 'MaxT',start_year_inclusive, end_year_inclusive)
        end_time = time.time()
        print("time taken " + str(end_time - start_time))
        num_clear(workpath, HV, start_year_inclusive, end_year_inclusive)


if __name__ == '__main__':
    #https://stackabuse.com/command-line-arguments-in-python/
    parser = argparse.ArgumentParser()
    #(workpath,HV_List,start_year_inclusive,end_year_inclusive)
    parser.add_argument("--workpath", "-w", help="specify workpath directory")
    parser.add_argument("--hv_tilename", "-t", help="specify tilename singular")
    parser.add_argument("--start_year_inclusive", "-s", help="specify start_year_inclusive")
    parser.add_argument("--end_year_inclusive", "-e", help="specify end_year_inclusive")

    # Read arguments from the command line
    args = parser.parse_args()

    if args.workpath and args.hv_tilename and args.start_year_inclusive and args.end_year_inclusive:
        print(str(args))
        st_stats_forHua(args.workpath,[args.hv_tilename],args.start_year_inclusive,args.end_year_inclusive)
        

        





