'''
/***************************************************************************
Name	   :  geowrsi_utilities.py
Description:  GeoWRSI Utility Functions for FEWSTools plugin
copyright  :  (C) 2022-2023 by FEWS
email      :  dhackman@contractor.usgs.gov
Created    :  07/29/2022 - dhackman
Modified   :


***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
'''
import numpy as np
import os

from PyQt5.QtWidgets import QMessageBox

from fews_tools import fews_tools_config as config
from fews_tools.utilities import geoclim_gdal_utilities as gdal_util


def fill_data_cube_from_files(reg_dic, file_list, resampled_output_path):
    '''
    This function will fill a 3-D cube with the data from a list of file(s)
    params(dic) - reg_dic - region dictionary
    params(list) - file_list - List of the files to be converted to data cube
    params(string) - resampled_output_path - path to resampled files output
    return(3-D cube) - data_cube
    '''
    # call gdal_util resample_input_files
    _, resampled_file_list = gdal_util.resample_input_files(
        resampled_output_path, None, file_list, reg_dic)
    _, _, _, _, txt_data_type = gdal_util.get_geotiff_info(
        resampled_file_list[0])
    np_data_type = gdal_util.TYPE_DIC[txt_data_type]["NP"]
    # call gdal_util get_data_cube
    data_cube = gdal_util.get_data_cube(resampled_file_list, np_data_type)
    return data_cube


def fill_data_cube_from_value(value, resampled_mask_file, np_data_type):
    '''
    This will take in a single specified value to be turned into a data_cube
    params(int or float) - value - either a integer LGP or float WHC value
    params(string) - resampled_mask_file - path ot resampled mask for region
    params(datatype) - np_data_type = this is a value from the TYPE_DIC
                                    in geoclim_gdal_utilities.py
    '''
    _, row_ct, col_ct, _, _ = gdal_util.get_geotiff_info(resampled_mask_file)
    data_array = np.full((row_ct, col_ct), value, np_data_type)
    data_cube = data_array.reshape(1, row_ct, col_ct)
    return data_cube


def fill_lgp(set_dic, reg_dic, resampled_mask_file, resampled_output_path):
    '''
    Create an integer 3-D cube of LGP values taken from the settings dictionary
    params(dic) - set_dic - settings dictionary
    params(dic) - reg_dic - region dictionary
    params(string) - resampled_mask_file - path to resampled mask for region
    params(string) - resampled_output_path - path to resampled files output
    return(3-D cube) - lgp_data_cube
    '''
    if set_dic['lgp_type'] != 'Specified':  # means it is default or selected
        # use method from fill gdas/GeoRFE
        file_list = []
        file_list.append(set_dic['lgp_file'])
        lgp_data_cube = fill_data_cube_from_files(reg_dic,
                                                  file_list,
                                                  resampled_output_path)
    else:
        np_data_type = gdal_util.TYPE_DIC["Byte"]["NP"]
        lgp_data_cube = fill_data_cube_from_value(
            set_dic['lgp_period'], resampled_mask_file, np_data_type)
    return lgp_data_cube


def fill_whc(set_dic, reg_dic, resampled_mask_file, resampled_output_path):
    '''
    Create a float 3-D cube of WHC values taken from the settings dictionary
    params(dic) - set_dic - settings dictionary
    params(dic) - reg_dic - region dictionary
    params(string) - resampled_mask_file - path to resampled mask for region
    params(string) - resampled_output_path - path to resampled files output
    return(3-D cube) - whc_data_cube
    '''
    if set_dic['whc_type'] != 'Specified':  # means it is default or selected
        # use method from fill gdas/GeoRFE
        file_list = []
        file_list.append(set_dic['whc_file'])
        whc_data_cube = fill_data_cube_from_files(reg_dic,
                                                  file_list,
                                                  resampled_output_path)
    else:
        np_data_type = gdal_util.TYPE_DIC["Float32"]["NP"]
        whc_data_cube = fill_data_cube_from_value(
            set_dic['whc_mm'], resampled_mask_file, np_data_type)
    return whc_data_cube


def fill_sos(set_dic, reg_dic, resampled_mask_file, resampled_output_path):
    '''
    Create a 3-D cube of WHC values taken from the settings dictionary
    params(dic) - set_dic - settings dictionary
    params(dic) - reg_dic - region dictionary
    params(string) - resampled_mask_file - path to resampled mask for region
    params(string) - resampled_output_path - path to resampled files output
    return(3-D cube) - sos_data_cube
    '''
    sos_data_cube = None
    if set_dic['sos_type'] == 'Calculated' or\
            set_dic['sos_type'] == 'Climatological':
        file_list = []
        file_list.append(reg_dic['SOS'])
        sos_data_cube = fill_data_cube_from_files(reg_dic,
                                                  file_list,
                                                  resampled_output_path)
    elif set_dic['sos_type'] == 'Selected':
        file_list = []
        file_list.append(set_dic['sos_file'])
        sos_data_cube = fill_data_cube_from_files(reg_dic,
                                                  file_list,
                                                  resampled_output_path)
    elif set_dic['sos_type'] == 'Specified':
        np_data_type = gdal_util.TYPE_DIC["Byte"]["NP"]
        sos_data_cube = fill_data_cube_from_value(set_dic['sos_period'],
                                                  resampled_mask_file,
                                                  np_data_type)
    return sos_data_cube


def fill_mask(resampled_mask_file):
    '''
    Create a 3-D cube of mask values taken from the resampled mask
    params(string) - resampled_mask_file - path to resampled mask for region
    returns(3-D cube) - mask_data_cube
    '''
    file_list = []
    file_list.append(resampled_mask_file)
    np_data_type = gdal_util.TYPE_DIC["Byte"]["NP"]
    mask_data_cube = gdal_util.get_data_cube(file_list, np_data_type)
    return mask_data_cube


def get_kc_weights(crop_dic, lgp):
    '''
    This function is designed to take a number of dekads specifying the
    length of the growing season (LGP) and return a 1-D array of weights
    for each period
    params(dic) - crop_dic - crop dictionary
    params(int) - lgp - integer of the number of dekads in lgp
    return(float array) - kc_wt - kc weight value for each dekad in LGP
    '''
    # Fraction of the season at each dekad
    dek_frac = [float((i+1)/lgp) for i in range(lgp)]
    # Set time periods of growing season when kc changes
    # Initialize them all to the initial stage value
    kc_wt = [float(crop_dic['K1']) for i in range(lgp)]
    # Set kc values for each dekad of the growing season
    for i in range(lgp):
        # Vegatative stage
        if (dek_frac[i] > crop_dic['F1']) and (dek_frac[i] <= crop_dic['F2']):
            kc_wt[i] = (crop_dic['K1'] + (((crop_dic['K2'] - crop_dic['K1'])
                                           / (crop_dic['F2'] - crop_dic['F1']))
                                          * (dek_frac[i] - crop_dic['F1'])))
        # Reproductive Stage
        elif (dek_frac[i] > crop_dic['F2']) and\
                (dek_frac[i] <= crop_dic['F3']):
            kc_wt[i] = crop_dic['K2']
        # Ripening Stage
        elif dek_frac[i] > crop_dic['F3']:
            kc_wt[i] = (crop_dic['K2'] + (((crop_dic['K3'] - crop_dic['K2'])
                                           / (1.0 - crop_dic['F3']))
                                          * (dek_frac[i] - crop_dic['F3'])))
    return kc_wt


def get_rd_fraction(crop_dic, lgp):
    '''
    This function is designed to take a number of periods specifying the
    length of the growing season (LGP) and return a 1-D array of weights
    for each period
    params(dic) - crop_dic - crop dictionary
    params(int) - lgp - integer of the number of dekads in lgp
    return(float array) - rdf - root depth fraction value for each dekad in LGP
    '''
    # Fraction of the season at each dekad
    dek_frac = [float((i+1)/lgp) for i in range(lgp)]
    # Effective Root
    rdf = [float(1.0) for i in range(lgp)]
    # Set RDF values for each dekad of the growing season
    for i in range(lgp):
        if dek_frac[i] <= crop_dic['F2']:
            rdf[i] = crop_dic['RootIni'] + ((1-crop_dic['RootIni'])
                                            * (dek_frac[i] / crop_dic['F2']))
    return rdf


def get_critical_soil_water(rdf, whc, swf):
    '''
    Function gets the Critical Soil Water (SWC) for a pixel value
    params(float array) - rdf - root depth fraction value for each dekad in LGP
    params(float) - whc - the water holding capacity of that pixel
    params(float) - swf - soil water fraction
    '''
    swc = [float(entry * whc * swf) for entry in rdf]
    return swc


def stack_cubes(cubes_arr):
    '''
    This function will stack the needed ppt, pet, lgp, whc cubes in order
    to slice the same pixel for all of them at the same time to pass int
    the run_geowrsi functions
    params(array) - cubes_arr - cubes to stack
    returns(3-D cube) - stacked_cube - the overall cubes stacked together
    '''
    stacked_cube = None
    for cube in cubes_arr:
        if stacked_cube is None:
            stacked_cube = cube.astype(gdal_util.TYPE_DIC['Float32']['NP'])
        else:
            stacked_cube = np.concatenate((
                stacked_cube,
                cube.astype(gdal_util.TYPE_DIC['Float32']['NP'])), axis=0)
    return stacked_cube


def get_pet_c(pet_arr, kc_arr):
    '''
    PETc is the crop specific potential evapotranspiration after an adjustment
    is made to the reference crop potential evapotranspiration (PET) by the use
    of appropriate dekadal crop coefficients (Kc), make sure arrays are same
    shape and size
    params(float array) - pet_arr - array of pet data
    params(float array) - kc_arr - kc weights array
    returns(float array) - pet_c - crop specific evapotranspiration
    '''
    if len(pet_arr) == len(kc_arr):
        pet_c = np.multiply(pet_arr, kc_arr)
    else:
        # If these arent equal lengths then it is likely a current/forecast/
        # extended run so we base it off of the pet_arr length to get pet_c
        # array of the same length
        pet_len = len(pet_arr)
        pet_c = np.multiply(pet_arr[0:pet_len], kc_arr[0:pet_len])
    return pet_c


def get_swb_a(ppt_arr, pet_arr, whc):
    '''
    Establish the antecedent soil water balance to carry over into the first
    dekad of the season, otherwise known as the plant available water, make
    sure ppt and pet arrays are same size and shape, only pass in the 6 dekads
    before SOS
    params(float array) - ppt_arr - array of ppt data (6 dekads)
    params(float array) - pet_arr - array of pet data (6 dekads)
    params(float) - whc - value for the water holding capacity
    returns(float) - swb_a - preseason soil water balance acumulation
    '''
    kc_arr = [float(0.15) for pet in pet_arr]
    pet_c = get_pet_c(pet_arr, kc_arr)
    tmp_paw = np.subtract(ppt_arr, pet_c)
    swb_a = 0
    for x_l in tmp_paw:
        swb_a += x_l
        if swb_a < 0:
            swb_a = 0
        elif swb_a > whc:
            swb_a = whc
    return swb_a


def get_wrsi_input_type(radio_button):
    '''
    Get the precipitation or evapotranspiration type
    params(gui_obj) - radio_buttion - The radio button to query
    returns(string) - ret_val - Climatological or Actual
    '''
    ret_val = 'Climatological'
    if radio_button.isChecked():
        ret_val = 'Actual'
    return ret_val


def check_empty_raster_arr(arr):
    '''
    This function is used to check if an array is full of only nan values,
    used so that only the dekads with WRSI values are output for the user
    to look at and view the information
    params(2d array) - arr - the raster file being passed in
    returns(bool) - true if empty, false if not empty
    '''
    return np.isnan(arr).all()


def run_geowrsi_historical(stacked_arr, crop_dic, periods, init_per_of_season,
                           final_per_of_season, cross_year, periodicity):
    '''
    This function will take an array, that contains the data for a single
    pixel and will then calculate to wrsi for that pixel, assumes data has
    been checked before being run, so no data checks will be used in this
    function
    params(float array) - stacked_arr - stacked ppt, pet, lgp, whc, mask, sos
    params(dic) - crop_dic - crop dictionary
    params(int) - periods - number of periods for the ppt and pet parts of cube
    params(int) - init_per_of_season - initial period from region
    params(int) - final_per_of_season - final period from region
    params(bool) - cross_year - if cross year or not
    params(string) - periodicity - tells whether Dekadal or Pentadal
    returns(arr) - wrsi_arr - the array with wrsi value for each dekad
    returns(arr) - swi_arr - the array with swi value for each dekad
    returns(arr) - lgp_phen_arr - the array with % of LGP phenology completed
    returns(arr) - eos_output_arr - the array with all the values for EOS files
    returns(arr) - crop_stages_arr - the array for all the crop stage totals
    '''
    # Used for the special case that user chooses an extremely short lgp
    # and causes the perods var to be to small to fit all of the output
    # in the output arrays
    if periods < 16:
        output_arr_len = 16
    else:
        output_arr_len = periods
    # initialize wrsi, swi, lgp_phen array with None values
    wrsi_arr = np.empty(output_arr_len) * np.nan
    swi_arr = np.empty(output_arr_len) * np.nan
    lgp_phen_arr = np.empty(output_arr_len) * np.nan
    # Used for EOS output values
    eos_output_arr = np.empty(output_arr_len) * np.nan
    # Used for crop_stages output values
    crop_stages_arr = np.empty(output_arr_len) * np.nan

    if periodicity == 'Dekadal':
        shift = config.DEKAD_SHIFT
        yr_shift = config.DEKADS_PER_YEAR
        no_start = config.DEKAD_NOSTART
    elif periodicity == 'Pentadal':
        shift = config.PENTAD_SHIFT
        yr_shift = config.PENTADS_PER_YEAR
        no_start = config.PENTAD_NOSTART
    # Check if within the mask region otherwise don't need to run the WRSI
    if int(stacked_arr[-1]) == 1:
        sos = int(stacked_arr[-2])
        whc = stacked_arr[-3]
        lgp = int(stacked_arr[-4])
        ppt_arr = stacked_arr[0:periods].astype(
            gdal_util.TYPE_DIC['Int16']['NP'])
        pet_arr = stacked_arr[periods:-
                              4].astype(gdal_util.TYPE_DIC['Float32']['NP'])

        if sos != no_start and sos != 0 and lgp != 0 and whc != 0:
            sos_temp = sos
            # Logic to handle seasons that extend over a cross year
            if cross_year:
                final_per_of_season = final_per_of_season + yr_shift
                if sos_temp < init_per_of_season:
                    sos_temp = sos_temp + yr_shift
            # Check that sos is within region initial and final period
            # if not set to no start/fail for pixel and return
            if not init_per_of_season <= sos_temp <= final_per_of_season:
                wrsi_arr[:] = config.FAIL_NO_START
                swi_arr[:] = config.FAIL_NO_START
                lgp_phen_arr[:] = config.FAIL_NO_START
                return wrsi_arr, swi_arr, lgp_phen_arr,\
                    eos_output_arr, crop_stages_arr
            # Logic to line up the sos and init/final_per_of_season with the
            # index values in the ppt_arr and pet_arr
            sos_idx = sos_temp - init_per_of_season + shift
            final_per_idx = final_per_of_season - init_per_of_season + shift
            kc = get_kc_weights(crop_dic, lgp)
            rdf = get_rd_fraction(crop_dic, lgp)
            swc = get_critical_soil_water(rdf, whc, crop_dic['KP'])
            pet_c = get_pet_c(pet_arr[sos_idx:sos_idx + lgp], kc)
            swb_a = get_swb_a(
                ppt_arr[sos_idx-shift:sos_idx],
                pet_arr[sos_idx-shift:sos_idx], whc)

            np_data_type = gdal_util.TYPE_DIC["Float32"]["NP"]
            aetout = np.full(lgp, 0.0, np_data_type)
            swi_processing_arr = np.full(lgp, 0.0, np_data_type)
            paw1 = 0.0
            paw2 = 0.0
            aet1 = 0.0
            aet2 = 0.0
            aets = 0.0
            pets = 0.0
            paw1_arr = np.full(lgp, 0.0, np_data_type)
            paw2_arr = np.full(lgp, 0.0, np_data_type)
            aet1_arr = np.full(lgp, 0.0, np_data_type)
            aet2_arr = np.full(lgp, 0.0, np_data_type)
            swb_a_arr = np.full(lgp, 0.0, np_data_type)
            wrsi_processing_arr = np.full(lgp, 0.0, np_data_type)
            lgp_phen_processing_arr = np.full(lgp, 0.0, np_data_type)
            surplus_water_total = 0.0
            max_surplus_water = 0.0
            water_deficit_total = 0.0
            max_water_deficit = 0.0
            aet_init = 0.0
            aet_veg = 0.0
            aet_flow = 0.0
            aet_ripe = 0.0
            wreq_init = 0.0
            wreq_veg = 0.0
            wreq_flow = 0.0
            wreq_ripe = 0.0
            wsurp_init = 0.0
            wsurp_veg = 0.0
            wsurp_flow = 0.0
            wsurp_ripe = 0.0
            wdef_init = 0.0
            wdef_veg = 0.0
            wdef_flow = 0.0
            wdef_ripe = 0.0

            for i in range(sos_idx, sos_idx+lgp):
                loop_index = i-sos_idx
                mad = swc[loop_index]
                paw1 = swb_a + ppt_arr[i]
                if paw1 >= mad:
                    aet1 = pet_c[loop_index]
                    if aet1 > paw1:
                        aet1 = paw1
                    paw2 = paw1-aet1
                    if paw2 >= mad:
                        aet2 = pet_c[loop_index]
                        if aet2 > paw2:
                            aet2 = paw2
                    else:
                        if (0.0 < paw2 < mad):
                            aet2 = (paw2/mad)*pet_c[loop_index]
                            if aet2 > paw2:
                                aet2 = paw2
                        else:
                            aet2 = 0.0
                else:
                    aet1 = (paw1/mad)*pet_c[loop_index]
                    if aet1 > paw1:
                        aet1 = paw1
                    paw2 = paw1 - aet1
                    if paw2 > 0.0:
                        aet2 = (paw2/mad)*pet_c[loop_index]
                        if aet2 > paw2:
                            aet2 = paw2
                        else:
                            aet2 = 0.0
                    else:
                        aet2 = 0.0

                paw1_arr[loop_index] = paw1
                paw2_arr[loop_index] = paw2
                aet1_arr[loop_index] = aet1
                aet2_arr[loop_index] = aet2

                aet = (crop_dic['C1'] * aet1) + (crop_dic['C2'] * aet2)

                swb_a = 0
                tmp_paw = paw1 - aet
                surplus_water = 0.0
                if tmp_paw > whc:
                    surplus_water = tmp_paw - whc
                    swb_a = whc
                else:
                    swb_a = tmp_paw
                swb_a_arr[loop_index] = swb_a
                swi_processing_arr[loop_index] = round((swb_a/whc)*100)
                # Calculate per period wrsi
                # Calculate the sum of AET and sum of PET for seasonal WRSI
                aets = aets + aet
                pets = pets + pet_c[loop_index]
                aetout[loop_index] = aet
                if aets > pets:
                    wrsi_processing_arr[loop_index] = 100
                else:
                    wrsi_processing_arr[loop_index] = round((aets/pets)*100)

                lgp_phen_processing_arr[loop_index] = round(
                    ((loop_index+1)/lgp)*100)

                # Surplus water processing
                if surplus_water > 0:
                    surplus_water_total += surplus_water
                    if surplus_water > max_surplus_water:
                        max_surplus_water = surplus_water
                # Water Deficit processing
                water_deficit = aet - pet_c[loop_index]
                if water_deficit < 0:
                    water_deficit_total -= water_deficit
                    if (-1 * water_deficit) > max_water_deficit:
                        max_water_deficit = -1 * water_deficit

                percent_season = ((loop_index+1)/lgp)
                # Crop stage totals calculated here
                if percent_season < crop_dic['F1']:
                    # Initial Stage
                    aet_init += aet
                    wreq_init += pet_c[loop_index]
                    if surplus_water > 0:
                        wsurp_init += surplus_water
                    if water_deficit < 0:
                        wdef_init -= water_deficit
                elif percent_season >= crop_dic['F1'] and\
                        percent_season < crop_dic['F2']:
                    # Vegetative stage
                    aet_veg += aet
                    wreq_veg += pet_c[loop_index]
                    if surplus_water > 0:
                        wsurp_veg += surplus_water
                    if water_deficit < 0:
                        wdef_veg -= water_deficit
                elif percent_season >= crop_dic['F2'] and\
                        percent_season < crop_dic['F3']:
                    # Flowering Stage
                    aet_flow += aet
                    wreq_flow += pet_c[loop_index]
                    if surplus_water > 0:
                        wsurp_flow += surplus_water
                    if water_deficit < 0:
                        wdef_flow -= water_deficit
                elif percent_season >= crop_dic['F3']:
                    # Ripening stage
                    aet_ripe += aet
                    wreq_ripe += pet_c[loop_index]
                    if surplus_water > 0:
                        wsurp_ripe += surplus_water
                    if water_deficit < 0:
                        wdef_ripe -= water_deficit

            # Outputs at every period
            # If season has not started then fill that part of the array with
            # YET_TO_START value (254)
            wrsi_arr[0:sos_idx] = config.YET_TO_START
            wrsi_arr[sos_idx:sos_idx + lgp] = wrsi_processing_arr
            wrsi_arr[sos_idx + lgp:] = wrsi_processing_arr[-1]
            # Stop WRSI calculation after Last Dek of Season + 6
            wrsi_arr[final_per_idx + shift:] = wrsi_arr[final_per_idx + shift]

            swi_arr[0:sos_idx] = config.YET_TO_START
            swi_arr[sos_idx:sos_idx + lgp] = swi_processing_arr
            swi_arr[sos_idx + lgp:] = swi_processing_arr[-1]
            # Stop SWI calculation after Last Dek of Season + 6
            swi_arr[final_per_idx + shift:] = swi_arr[final_per_idx + shift]

            lgp_phen_arr[0:sos_idx] = config.YET_TO_START
            lgp_phen_arr[sos_idx:sos_idx + lgp] = lgp_phen_processing_arr
            lgp_phen_arr[sos_idx + lgp:] = 200
            # Stop lgp_phen calculation after Last Dek of Season + 6
            lgp_phen_arr[final_per_idx +
                         shift:] = lgp_phen_arr[final_per_idx + shift]

            # EOS outputs section
            eos_output_arr[0] = round(wrsi_arr[final_per_idx + shift])
            eos_output_arr[2] = round(lgp_phen_arr[final_per_idx + shift])
            eos_output_arr[3] = round(swi_arr[final_per_idx + shift])
            eos_output_arr[4] = round(pets)
            eos_output_arr[5] = round(max_water_deficit)
            eos_output_arr[6] = round(water_deficit_total)
            eos_output_arr[7] = round(max_surplus_water)
            eos_output_arr[8] = round(surplus_water_total)
            eos_output_arr[9] = round(aets)

            # Crop stages outputs secton
            crop_stages_arr[0] = round(aet_init)
            crop_stages_arr[1] = round(aet_veg)
            crop_stages_arr[2] = round(aet_flow)
            crop_stages_arr[3] = round(aet_ripe)
            crop_stages_arr[4] = round(wreq_init)
            crop_stages_arr[5] = round(wreq_veg)
            crop_stages_arr[6] = round(wreq_flow)
            crop_stages_arr[7] = round(wreq_ripe)
            crop_stages_arr[8] = round(wsurp_init)
            crop_stages_arr[9] = round(wsurp_veg)
            crop_stages_arr[10] = round(wsurp_flow)
            crop_stages_arr[11] = round(wsurp_ripe)
            crop_stages_arr[12] = round(wdef_init)
            crop_stages_arr[13] = round(wdef_veg)
            crop_stages_arr[14] = round(wdef_flow)
            crop_stages_arr[15] = round(wdef_ripe)
        else:
            if sos == 0 or lgp == 0 or whc == 0:
                wrsi_arr[:] = config.NO_DATA
                swi_arr[:] = config.NO_DATA
                lgp_phen_arr[:] = config.NO_DATA

                # EOS outputs section
                eos_output_arr[0] = config.NO_DATA
                eos_output_arr[2] = config.NO_DATA
                eos_output_arr[3] = config.NO_DATA
            else:
                wrsi_arr[:] = config.FAIL_NO_START
                swi_arr[:] = config.FAIL_NO_START
                lgp_phen_arr[:] = config.FAIL_NO_START
                # EOS outputs section
                eos_output_arr[0] = config.FAIL_NO_START
                eos_output_arr[2] = config.FAIL_NO_START
                eos_output_arr[3] = config.FAIL_NO_START

    return wrsi_arr, swi_arr, lgp_phen_arr, eos_output_arr, crop_stages_arr


def calc_sos(stacked_cube, reg_dic, sos_dic, periods_early, ppt_cube_periods):
    '''
    Calculates the sos per pixel for the given data

    :params(float array) - Stacked_cube - stacked ppt, pet, lgp, whc, climsos,
                                          and mask data
    :params(dic) - reg_dic - Dictionary of region values
    :params(dic) - sos_dic - Dictionary of geowrsi settings
    :params(int) - periods_early - Number of periods before the inital period
                                   the ppt data starts
    :params(int) - ppt_cube_periods - number of periods in the ppt cube

    :returns(int array) - stacked_cube - Returns the origional stacked_cube
                                         with the sos slice replaced with the
                                         newly calculated sos
    '''

    if reg_dic["PeriodType"] == "Dekadal":
        no_start = config.DEKAD_NOSTART
        max_periods = config.DEKADS_PER_YEAR

    elif reg_dic["PeriodType"] == "Pentadal":
        no_start = config.PENTAD_NOSTART
        max_periods = config.PENTADS_PER_YEAR

    else:
        QMessageBox.warning(
            None,
            'Period Type Missing!!',
            'No period type selected, choose a differnt region',
            QMessageBox.Ok)

    # Calculate the 0-Based SOS Period for each good pixel
    sos = np.apply_along_axis(sos_1d, 0, stacked_cube,
                              reg_dic["InitialPeriod"], reg_dic["FinalPeriod"],
                              no_start, max_periods, sos_dic, periods_early,
                              ppt_cube_periods)

    stacked_cube[-2] = sos

    return stacked_cube


def sos_1d(stacked_arr, region_sos, region_eos, no_start, max_periods, sos_dic,
           periods_early, ppt_cube_periods):
    '''
    Calculates an integer representing the period in which a pixel's growing
    season begins based on when/if a rainfall threashold is met.

    :params(float array) - stacked_arr - A slice along the time axis of the
                                         stacked ppt, pet, lgp, whc, mask data
    :params(int) - region_sos - Earliest period of the season for the given
    :params(int) - region_eos - Last period of the season for the given region
    :params(int) - no_start - No start value for the given period type
    :params(int) - max_periods - Max number of periods in a season
    :params(dic) - sos_dic - Dictionary of sos settings values
    :params(int) - periods_early - Number of periods before the inital period
                                   the ppt data starts
    :params(int) - ppt_cube_periods - number of periods in the ppt cube

    :returns(int) - sos - Calculated start of season value for the given pixel
    '''

    mediansos = int(stacked_arr[-2])
    lgp = int(stacked_arr[-4])

    # Check if within the mask region, mediansos < no_start, and lgp != 0
    if stacked_arr[-1] == 1 and mediansos < no_start and lgp != 0:

        rain_arr = stacked_arr[0:ppt_cube_periods].astype(
            gdal_util.TYPE_DIC['Int16']['NP'])

        sos = no_start  # default is no start
        earliest_sos_ind = region_sos - 1
        region_eos_ind = region_eos - 1
        medsosdek_ind = mediansos - 1

        # Account for the ppt files starting periods_early
        # before the region sos
        ppt_file_shift = region_sos - periods_early

        # Check if the region is a cross year and adjust the eos and medsos
        if region_eos < region_sos:
            region_eos_ind += max_periods

            if medsosdek_ind < region_sos:
                medsosdek_ind += max_periods

        if sos_dic['ignore_sos_clim'] == "True":
            start_looking = earliest_sos_ind
            stop_looking = region_eos_ind

        else:
            # sos cannot occur before the region's Earliest sos period OR
            # sooner than max_periods_early periods before Clim sos
            if ((medsosdek_ind - int(sos_dic['max_periods_early']))
                    >= earliest_sos_ind and medsosdek_ind < no_start):
                start_looking = medsosdek_ind - \
                    int(sos_dic['max_periods_early'])
            else:
                start_looking = earliest_sos_ind

            # sos is considered 'No Start' if not found within max_periods_late
            # periods of the Climatological sos
            if ((medsosdek_ind + int(sos_dic['max_periods_late']))
                    < region_eos_ind and medsosdek_ind < no_start):
                stop_looking = medsosdek_ind + int(sos_dic['max_periods_late'])
            else:
                stop_looking = region_eos_ind

        periods_to_search = stop_looking - start_looking

        # The starting period to search taking into account that the ppt data
        # starts at earliest_sos_ind - periods_early
        start_sos_search = periods_early + start_looking - earliest_sos_ind
        end_sos_search = start_sos_search + periods_to_search + 1

        for i in range(start_sos_search, end_sos_search):
            if ((rain_arr[i] >= int(sos_dic['period_1_rainfall'])) and
                ((rain_arr[i + 1] + rain_arr[i + 2]) >=
                    int(sos_dic['period_2_3_rainfall']))):

                sos = i + ppt_file_shift

                if sos != no_start and sos > max_periods:
                    sos -= max_periods

                break

        # Exclude areas that will be less than X% complete by EoS
        # (max_growing_periods_after_eos after region's final dek)
        if sos_dic['exclude_areas_less_than'] == "True":
            if sos not in [0, no_start]:
                if 0 < ((region_eos + sos_dic['max_growing_periods_after_eos']
                        - sos) / lgp * 100.0) < \
                        float(sos_dic['exclude_areas_less_than_percent']):
                    sos = no_start

    else:
        sos = 0

    return sos


def anomaly_calc(array_orig, reg_dic, resampled_mask,
                 resample_clim_array_orig, clim_type, cross_year):
    '''
    Function for calculating anomaly for WRSI or SOS depending on the
    clim_type input
    params(2-D array) - array_orig - input array to be compared against the
                                    historic median
    params(dic) - reg_dic - region dictionary
    params(2-D array) - resampled_mask - resampled mask for the region
    params(2-D array) - resampled_clim_array_orig - climatological WRSI or SOS
                                                array resampled for the region
                                                being worked with
    params(string) - clim_type - "WRSI" or "SOS" depending on the calculation
                                desired
    params(bool) - cross_year - if cross year or not
    return(2-D array) - anom - the calculated anomaly
    '''
    anom = None
    # Create a copy of array so that the original array is not changed
    array = array_orig.copy()
    resample_clim_array = resample_clim_array_orig.copy()
    if reg_dic["PeriodType"] == 'Dekadal':
        yr_shift = config.DEKADS_PER_YEAR
        mid = yr_shift/2
        no_start = config.DEKAD_NOSTART
    elif reg_dic["PeriodType"] == 'Pentadal':
        yr_shift = config.PENTADS_PER_YEAR
        mid = yr_shift/2
        no_start = config.PENTAD_NOSTART

    if clim_type == "WRSI":
        anom = calc_masked_pct_of_avg_wrsi(array, resample_clim_array,
                                                resampled_mask)
        anom = anom_fix(anom, array, [config.FAIL_NO_START,
                                      config.YET_TO_START,
                                      config.NO_DATA])
    elif clim_type == "SOS":
        if cross_year:
            # add yr_shift to values less than mid when it is a cross year
            array[array <= mid] += yr_shift
            resample_clim_array[resample_clim_array <= mid] += yr_shift

        anom = gdal_util.calc_masked_array_diff(array, resample_clim_array,
                                                resampled_mask)
        anom = anom_fix(anom, array, [no_start])

    return anom


def anom_fix(anom, array, carryover_values):
    '''
    Function to carry over specified values to the anomaly

    params(2-D array) - anom - the anomaly array
    params(2-D array) - array - the WRSI or SOS array
    params(array) - carryover_values - The values that should replace the
                                       result in the anom array if they are
                                       present for a given pixel in the input
                                       array

    return (2_D array) - anom - the resulting anomaly array
    '''
    for j in range(anom.shape[1]):
        for i in range(anom.shape[0]):
            if array[i, j] in carryover_values:
                anom[i, j] = array[i, j]

    return anom


def get_default_output_dic():
    '''
    Function to setup the default output settings dictionary for use in the
    output_settings_controller and the geowrsi_controller
    returns(dic) - output_dic - the output settings dictionary
    '''
    output_dic = {  # Fill with defaults to start
        'current_wrsi': True,
        'current_wrsi_anomaly': False,
        'current_lgp_phenology': True,
        'current_soil_water_index': True,
        'current_water_req_totals': False,
        'current_max_water_deficits': False,
        'current_water_deficit_totals': False,
        'current_max_water_surplus': False,
        'current_water_surplus_totals': False,
        'current_aet_totals': False,
        'eos_wrsi': True,
        'eos_wrsi_anomaly': True,
        'eos_lgp_phenology': False,
        'eos_soil_water_index': False,
        'eos_water_req_totals': False,
        'eos_max_water_deficits': False,
        'eos_water_deficit_totals': False,
        'eos_max_water_surplus': False,
        'eos_water_surplus_totals': False,
        'eos_aet_totals': False,
        'keep_temp_wb_files': False,
        'keep_until': 0.0,
        'output_directory': '',
        'save_every_wrsi': False,
        'save_every_lgp_phenology': False,
        'save_every_soil_water_index': False,
        'save_all_files_at_given': False,
        'given_period': 0,
        'save_crop_stage_totals_eos': True,
        'format_outputs_for_arcgis': False,
        'show_title': True,
        'show_full_title': False,
        'font_name': '',
        'font_size': 0,
        'font_style': '',
        'output_prefix': ''
    }
    return output_dic


def set_output_dic_from_file(file_path, output_dic):
    '''
    If the config.GEOWRSI_OUT_SETTINGS_FILE file exists, then fill the
    output_dic with the values from that file in order to initialize settings
    params(string) - file_path - file path to geowrsi_output_settings.txt
    params(dic) - output_dic - this is the current output_dic, it will be
                            overwritten and then returned if file path exists
    returns(dic) - output_dic - returns either the original dic that was passed
                                in or the new updated dic
    '''
    if os.path.exists(file_path):
        lines = []
        with open(file_path, "r") as f_obj:
            lines = f_obj.readlines()
        for line in lines:
            temp_val = line.split(' ', 1)
            if temp_val[1].strip() == "True":
                output_dic[temp_val[0]] = True
            elif temp_val[1].strip() == "False":
                output_dic[temp_val[0]] = False
            else:
                output_dic[temp_val[0]] = temp_val[1].strip()
    return output_dic


def get_default_settings_dic():
    '''
    Function to setup the default settings dictionary for use in the
    geowrsi_settings_controller
    returns(dic) - set_dic - the settings dictionary
    '''
    set_dic = {
        # Region/crop tab
        'analysis_region': '',
        'crop_type': '',
        # sos tab
        'sos_type': 'Calculated',
        'sos_offset': 0,
        'sos_period': 1,
        'sos_file': '',
        # lgp tab
        'lgp_type': 'Default',
        'lgp_period': 1,
        'lgp_file': '',
        # WHC Tab
        'whc_type': 'Default',
        'whc_mm': 0.0,
        'whc_file': '',
        # Mask Tab
        'mask_type': 'Default',
        'mask_file': '',
        # Precip/PET data
        'precip_type': 'Actual',
        'pet_type': 'Actual',
        'use_simulated_data': False,
        'ppt_dataset': '',
        'pet_dataset': ''
    }
    return set_dic


def set_settings_dic_from_file(file_path, settings_dic):
    '''
    If the config.GEOWRSI_SETTINGS_FILE file exists, then fill the
    settings_dic with the values from that file in order to initialize settings
    params(string) - file_path - file path to geowrsi_settings.txt
    params(dic) - settings_dic - this is the current settings_dic, it will be
                            overwritten and then returned if file path exists
    returns(dic) - settings_dic - returns either the original dic that was passed
                                in or the new updated dic
    '''
    if os.path.exists(file_path):
        lines = []
        with open(file_path, "r") as f_obj:
            lines = f_obj.readlines()
        for line in lines:
            temp_val = line.split(' ', 1)
            if temp_val[0] in settings_dic.keys():
                if temp_val[1].strip() == "True":
                    settings_dic[temp_val[0]] = True
                elif temp_val[1].strip() == "False":
                    settings_dic[temp_val[0]] = False
                else:
                    settings_dic[temp_val[0]] = temp_val[1].strip()
    return settings_dic


def calculate_max_sos_plus_lgp(cross_year, periodicity, output_file_path,
                               set_dic, reg_dic):
    '''
    Function to calculate the number of files needed based on SOS and LGP
    params(bool) - cross_year - if it is a cross year
    params(string) - periodicity - Dekadal or Pentadal
    params(string) - output_file_path - output file path from wrksp_setup
    params(dict) - set_dic - the geowrsi settings dictionary
    params(dict) - reg_dic - the region dictionary for selected region
    returns(int) - max_value - integer of the max combined pixel for SOS = LGP
    '''
    # Set  values based on periodicity
    if periodicity == 'Dekadal':
        shift = config.DEKADS_PER_YEAR
        half_year_val = shift/2
    elif periodicity == 'Pentadal':
        shift = config.PENTADS_PER_YEAR
        half_year_val = shift/2
    # Resample mask to reg_dic to pass into the get sos and lgp cubes
    resampled_mask, _ = \
        gdal_util.resample_input_files(output_file_path,
                                       set_dic['mask_file'],
                                       [], reg_dic)
    # Get the SOS cube
    sos_data_cube = fill_sos(set_dic, reg_dic, resampled_mask,
                             output_file_path)
    # If cross year add shift (36 or 72) for SOS values in first half of year
    # but exclude 0 since that is a no data spot
    if cross_year:
        sos_data_cube = np.where(np.logical_and(sos_data_cube <= half_year_val,
                                                0 < sos_data_cube),
                                 sos_data_cube + shift,
                                 sos_data_cube)
    # If the sos type is calculated, then to be safe we need to add on the
    # max_periods_late value to this cube since it is using the clim avg sos
    # file for the sos cube before we have done the entire calc sos function
    if set_dic['sos_type'] == 'Calculated':
        sos_data_cube = sos_data_cube + int(set_dic['max_periods_late'])

    # Get the LGP cube
    lgp_data_cube = fill_lgp(set_dic, reg_dic, resampled_mask,
                             output_file_path)
    # Add SOS and LGP cubes together
    sos_plus_lgp = np.add(sos_data_cube, lgp_data_cube)
    # Get the max of SOS + LGP
    max_value = np.max(sos_plus_lgp)
    # Return max value as an int
    return max_value


def run_geowrsi_current(stacked_arr, crop_dic, periods, init_per_of_season,
                        final_per_of_season, cross_year, periodicity):
    '''
    This function will take an array, that contains the data for a single
    pixel and will then calculate to wrsi for that pixel, assumes data has
    been checked before being run, so no data checks will be used in this
    function
    params(float array) - stacked_arr - stacked ppt, pet, lgp, whc, mask, sos
    params(dic) - crop_dic - crop dictionary
    params(int) - periods - number of periods for the ppt and pet parts of cube
    params(int) - init_per_of_season - initial period from region
    params(int) - final_per_of_season - final period from region
    params(bool) - cross_year - if cross year or not
    params(string) - periodicity - tells whether Dekadal or Pentadal
    returns(arr) - wrsi_arr - the array with wrsi value for each dekad
    returns(arr) - swi_arr - the array with swi value for each dekad
    returns(arr) - lgp_phen_arr - the array with % of LGP phenology completed
    returns(arr) - cur_output_arr - the array with all the Current period 
                                    outputs
    '''
    # Used for the special case that user chooses an extremely short lgp or
    # since its current the PET list length which sets the periods is short
    # and causes the perods var to be to small to fit all of the output
    # in the output arrays
    if periods < 10:
        output_arr_len = 10
    else:
        output_arr_len = periods
    # initialize wrsi, swi, lgp_phen array with None values
    wrsi_arr = np.empty(output_arr_len) * np.nan
    swi_arr = np.empty(output_arr_len) * np.nan
    lgp_phen_arr = np.empty(output_arr_len) * np.nan
    # Used for current output values
    cur_output_arr = np.empty(output_arr_len) * np.nan

    if periodicity == 'Dekadal':
        shift = config.DEKAD_SHIFT
        yr_shift = config.DEKADS_PER_YEAR
        no_start = config.DEKAD_NOSTART
    elif periodicity == 'Pentadal':
        shift = config.PENTAD_SHIFT
        yr_shift = config.PENTADS_PER_YEAR
        no_start = config.PENTAD_NOSTART
    # Check if within the mask region otherwise don't need to run the WRSI
    if int(stacked_arr[-1]) == 1:
        sos = int(stacked_arr[-2])
        whc = stacked_arr[-3]
        lgp = int(stacked_arr[-4])
        ppt_arr = stacked_arr[0:periods].astype(
            gdal_util.TYPE_DIC['Int16']['NP'])
        pet_arr = stacked_arr[periods:-
                              4].astype(gdal_util.TYPE_DIC['Float32']['NP'])

        if sos != no_start and sos != config.SOS_YET_TO_START and\
                sos != 0 and lgp != 0 and whc != 0:
            sos_temp = sos
            # Logic to handle seasons that extend over a cross year
            if cross_year:
                final_per_of_season = final_per_of_season + yr_shift
                if sos_temp < init_per_of_season:
                    sos_temp = sos_temp + yr_shift
            # Check that sos is within region initial and final period
            # if not set to no start/fail for pixel and return
            if not init_per_of_season <= sos_temp <= final_per_of_season:
                wrsi_arr[:] = config.FAIL_NO_START
                swi_arr[:] = config.FAIL_NO_START
                lgp_phen_arr[:] = config.FAIL_NO_START
                return wrsi_arr, swi_arr, lgp_phen_arr,\
                    cur_output_arr
            # Logic to line up the sos and init/final_per_of_season with the
            # index values in the ppt_arr and pet_arr
            sos_idx = sos_temp - init_per_of_season + shift
            # If the sos_idx is >= periods then the SOS is beyond the current
            # dekad so set the outputs to yet to start
            if sos_idx >= periods:
                wrsi_arr[:] = config.YET_TO_START
                swi_arr[:] = config.YET_TO_START
                lgp_phen_arr[:] = config.YET_TO_START
                cur_output_arr[0] = config.YET_TO_START
                cur_output_arr[2] = config.YET_TO_START
                cur_output_arr[3] = config.YET_TO_START
                return wrsi_arr, swi_arr, lgp_phen_arr,\
                    cur_output_arr

            final_per_idx = final_per_of_season - init_per_of_season + shift
            kc = get_kc_weights(crop_dic, lgp)
            rdf = get_rd_fraction(crop_dic, lgp)
            swc = get_critical_soil_water(rdf, whc, crop_dic['KP'])
            # if sos_idx + lgp is >= periods then we need to only run up
            # to the current period to avoid index out of bounds errors
            # else run from sos_idx to sos_idx + lgp
            if (sos_idx + lgp) >= periods:
                cur_per_idx = periods
            else:
                cur_per_idx = sos_idx + lgp
            pet_c = get_pet_c(pet_arr[sos_idx:cur_per_idx], kc)
            swb_a = get_swb_a(
                ppt_arr[sos_idx-shift:sos_idx],
                pet_arr[sos_idx-shift:sos_idx], whc)

            proc_arr_len = cur_per_idx - sos_idx
            np_data_type = gdal_util.TYPE_DIC["Float32"]["NP"]
            aetout = np.full(proc_arr_len, 0.0, np_data_type)
            swi_processing_arr = np.full(proc_arr_len, 0.0, np_data_type)
            paw1 = 0.0
            paw2 = 0.0
            aet1 = 0.0
            aet2 = 0.0
            aets = 0.0
            pets = 0.0
            paw1_arr = np.full(proc_arr_len, 0.0, np_data_type)
            paw2_arr = np.full(proc_arr_len, 0.0, np_data_type)
            aet1_arr = np.full(proc_arr_len, 0.0, np_data_type)
            aet2_arr = np.full(proc_arr_len, 0.0, np_data_type)
            swb_a_arr = np.full(proc_arr_len, 0.0, np_data_type)
            wrsi_processing_arr = np.full(proc_arr_len, 0.0, np_data_type)
            lgp_phen_processing_arr = np.full(proc_arr_len, 0.0, np_data_type)
            surplus_water_total = 0.0
            max_surplus_water = 0.0
            water_deficit_total = 0.0
            max_water_deficit = 0.0

            for i in range(sos_idx, cur_per_idx):
                loop_index = i-sos_idx
                mad = swc[loop_index]
                paw1 = swb_a + ppt_arr[i]
                if paw1 >= mad:
                    aet1 = pet_c[loop_index]
                    if aet1 > paw1:
                        aet1 = paw1
                    paw2 = paw1-aet1
                    if paw2 >= mad:
                        aet2 = pet_c[loop_index]
                        if aet2 > paw2:
                            aet2 = paw2
                    else:
                        if (0.0 < paw2 < mad):
                            aet2 = (paw2/mad)*pet_c[loop_index]
                            if aet2 > paw2:
                                aet2 = paw2
                        else:
                            aet2 = 0.0
                else:
                    aet1 = (paw1/mad)*pet_c[loop_index]
                    if aet1 > paw1:
                        aet1 = paw1
                    paw2 = paw1 - aet1
                    if paw2 > 0.0:
                        aet2 = (paw2/mad)*pet_c[loop_index]
                        if aet2 > paw2:
                            aet2 = paw2
                        else:
                            aet2 = 0.0
                    else:
                        aet2 = 0.0

                paw1_arr[loop_index] = paw1
                paw2_arr[loop_index] = paw2
                aet1_arr[loop_index] = aet1
                aet2_arr[loop_index] = aet2

                aet = (crop_dic['C1'] * aet1) + (crop_dic['C2'] * aet2)

                swb_a = 0
                tmp_paw = paw1 - aet
                surplus_water = 0.0
                if tmp_paw > whc:
                    surplus_water = tmp_paw - whc
                    swb_a = whc
                else:
                    swb_a = tmp_paw
                swb_a_arr[loop_index] = swb_a
                swi_processing_arr[loop_index] = round((swb_a/whc)*100)
                # Calculate per period wrsi
                # Calculate the sum of AET and sum of PET for seasonal WRSI
                aets = aets + aet
                pets = pets + pet_c[loop_index]
                aetout[loop_index] = aet
                if aets > pets:
                    wrsi_processing_arr[loop_index] = 100
                else:
                    wrsi_processing_arr[loop_index] = round((aets/pets)*100)

                lgp_phen_processing_arr[loop_index] = round(
                    ((loop_index+1)/lgp)*100)

                # Surplus water processing
                if surplus_water > 0:
                    surplus_water_total += surplus_water
                    if surplus_water > max_surplus_water:
                        max_surplus_water = surplus_water
                # Water Deficit processing
                water_deficit = aet - pet_c[loop_index]
                if water_deficit < 0:
                    water_deficit_total -= water_deficit
                    if (-1 * water_deficit) > max_water_deficit:
                        max_water_deficit = -1 * water_deficit

            # Outputs at every period
            # If season has not started then fill that part of the array with
            # YET_TO_START value (254)
            wrsi_arr[0:sos_idx] = config.YET_TO_START
            wrsi_arr[sos_idx:cur_per_idx] = wrsi_processing_arr
            wrsi_arr[cur_per_idx:] = wrsi_processing_arr[-1]
            # Stop WRSI calculation after Last Dek of Season + 6
            if (final_per_idx + shift) < output_arr_len:
                wrsi_arr[final_per_idx + shift:] = \
                    wrsi_arr[final_per_idx + shift]

            swi_arr[0:sos_idx] = config.YET_TO_START
            swi_arr[sos_idx:cur_per_idx] = swi_processing_arr
            swi_arr[cur_per_idx:] = swi_processing_arr[-1]
            # Stop SWI calculation after Last Dek of Season + 6
            if (final_per_idx + shift) < output_arr_len:
                swi_arr[final_per_idx + shift:] = \
                    swi_arr[final_per_idx + shift]

            lgp_phen_arr[0:sos_idx] = config.YET_TO_START
            lgp_phen_arr[sos_idx:cur_per_idx] = lgp_phen_processing_arr
            lgp_phen_arr[cur_per_idx:] = 200
            # Stop lgp_phen calculation after Last Dek of Season + 6
            if (final_per_idx + shift) < output_arr_len:
                lgp_phen_arr[final_per_idx +
                             shift:] = lgp_phen_arr[final_per_idx + shift]

            # Cur outputs section
            cur_output_arr[0] = round(wrsi_arr[-1])
            cur_output_arr[2] = round(lgp_phen_arr[-1])
            cur_output_arr[3] = round(swi_arr[-1])
            cur_output_arr[4] = round(pets)
            cur_output_arr[5] = round(max_water_deficit)
            cur_output_arr[6] = round(water_deficit_total)
            cur_output_arr[7] = round(max_surplus_water)
            cur_output_arr[8] = round(surplus_water_total)
            cur_output_arr[9] = round(aets)

        else:
            if sos == config.SOS_YET_TO_START and lgp != 0 and whc != 0:
                wrsi_arr[:] = config.YET_TO_START
                swi_arr[:] = config.YET_TO_START
                lgp_phen_arr[:] = config.YET_TO_START

                # CUR outputs section
                cur_output_arr[0] = config.YET_TO_START
                cur_output_arr[2] = config.YET_TO_START
                cur_output_arr[3] = config.YET_TO_START

            elif sos == 0 or lgp == 0 or whc == 0:
                wrsi_arr[:] = config.NO_DATA
                swi_arr[:] = config.NO_DATA
                lgp_phen_arr[:] = config.NO_DATA

                # CUR outputs section
                cur_output_arr[0] = config.NO_DATA
                cur_output_arr[2] = config.NO_DATA
                cur_output_arr[3] = config.NO_DATA
            else:
                wrsi_arr[:] = config.FAIL_NO_START
                swi_arr[:] = config.FAIL_NO_START
                lgp_phen_arr[:] = config.FAIL_NO_START
                # CUR outputs section
                cur_output_arr[0] = config.FAIL_NO_START
                cur_output_arr[2] = config.FAIL_NO_START
                cur_output_arr[3] = config.FAIL_NO_START

    return wrsi_arr, swi_arr, lgp_phen_arr, cur_output_arr


def calc_current_sos(stacked_cube, reg_dic, sos_dic, periods_early, ppt_cube_periods):
    '''
    Calculates the sos per pixel for the given data

    :params(float array) - Stacked_cube - stacked ppt, pet, lgp, whc, climsos,
                                          and mask data
    :params(dic) - reg_dic - Dictionary of region values
    :params(dic) - sos_dic - Dictionary of geowrsi settings
    :params(int) - periods_early - Number of periods before the inital period
                                   the ppt data starts
    :params(int) - ppt_cube_periods - number of periods in the ppt cube

    :returns(int array) - stacked_cube - Returns the origional stacked_cube
                                         with the sos slice replaced with the
                                         newly calculated sos
    '''

    if reg_dic["PeriodType"] == "Dekadal":
        no_start = config.DEKAD_NOSTART
        max_periods = config.DEKADS_PER_YEAR

    elif reg_dic["PeriodType"] == "Pentadal":
        no_start = config.PENTAD_NOSTART
        max_periods = config.PENTADS_PER_YEAR

    else:
        QMessageBox.warning(
            None,
            'Period Type Missing!!',
            'No period type selected, choose a differnt region',
            QMessageBox.Ok)

    # Calculate the 0-Based SOS Period for each good pixel
    sos = np.apply_along_axis(sos_1d_current, 0, stacked_cube,
                              reg_dic["InitialPeriod"], reg_dic["FinalPeriod"],
                              no_start, max_periods, sos_dic, periods_early,
                              ppt_cube_periods)

    stacked_cube[-2] = sos

    return stacked_cube


def sos_1d_current(stacked_arr, region_sos, region_eos, no_start, max_periods,
                   sos_dic, periods_early, ppt_cube_periods):
    '''
    Calculates an integer representing the period in which a pixel's growing
    season begins based on when/if a rainfall threashold is met.

    :params(float array) - stacked_arr - A slice along the time axis of the
                                         stacked ppt, pet, lgp, whc, mask data
    :params(int) - region_sos - Earliest period of the season for the given
    :params(int) - region_eos - Last period of the season for the given region
    :params(int) - no_start - No start value for the given period type
    :params(int) - max_periods - Max number of periods in a season
    :params(dic) - sos_dic - Dictionary of sos settings values
    :params(int) - periods_early - Number of periods before the inital period
                                   the ppt data starts
    :params(int) - ppt_cube_periods - number of periods in the ppt cube

    :returns(int) - sos - Calculated start of season value for the given pixel
    '''

    mediansos = int(stacked_arr[-2])
    lgp = int(stacked_arr[-4])

    # Check if within the mask region, mediansos < no_start, and lgp != 0
    if stacked_arr[-1] == 1 and mediansos < no_start and lgp != 0:

        rain_arr = stacked_arr[0:ppt_cube_periods].astype(
            gdal_util.TYPE_DIC['Int16']['NP'])

        sos = config.SOS_YET_TO_START  # default for current SOS_YET_TO_START
        earliest_sos_ind = region_sos - 1
        region_eos_ind = region_eos - 1
        medsosdek_ind = mediansos - 1

        # Account for the ppt files starting periods_early
        # before the region sos
        ppt_file_shift = region_sos - periods_early

        # Check if the region is a cross year and adjust the eos and medsos
        if region_eos < region_sos:
            region_eos_ind += max_periods

            if medsosdek_ind < region_sos:
                medsosdek_ind += max_periods

        if sos_dic['ignore_sos_clim'] == "True":
            start_looking = earliest_sos_ind
            stop_looking = region_eos_ind

        else:
            # sos cannot occur before the region's Earliest sos period OR
            # sooner than max_periods_early periods before Clim sos
            if ((medsosdek_ind - int(sos_dic['max_periods_early']))
                    >= earliest_sos_ind and medsosdek_ind < no_start):
                start_looking = medsosdek_ind - \
                    int(sos_dic['max_periods_early'])
            else:
                start_looking = earliest_sos_ind

            # sos is considered 'No Start' if not found within max_periods_late
            # periods of the Climatological sos
            if ((medsosdek_ind + int(sos_dic['max_periods_late']))
                    < region_eos_ind and medsosdek_ind < no_start):
                stop_looking = medsosdek_ind + int(sos_dic['max_periods_late'])
            else:
                stop_looking = region_eos_ind

        periods_to_search = stop_looking - start_looking

        # The starting period to search taking into account that the ppt data
        # starts at earliest_sos_ind - periods_early
        start_sos_search = periods_early + start_looking - earliest_sos_ind
        end_sos_search = start_sos_search + periods_to_search + 1

        # Add check for if our end_sos_search and the 2 periods needed for a
        # look ahead is beyond the length of rain_arr
        if (end_sos_search + 2) >= len(rain_arr):
            # set new end_sos_search to length of rain_arr - 2 for look ahead
            end_sos_search = len(rain_arr) - 2

        # Add check for if the end_sos_search ends up being less than the
        # start_sos_search
        if end_sos_search > start_sos_search:
            for i in range(start_sos_search, end_sos_search):
                if ((rain_arr[i] >= int(sos_dic['period_1_rainfall'])) and
                    ((rain_arr[i + 1] + rain_arr[i + 2]) >=
                        int(sos_dic['period_2_3_rainfall']))):

                    sos = i + ppt_file_shift

                    if sos != config.SOS_YET_TO_START and sos > max_periods:
                        sos -= max_periods

                    break

        # Exclude areas that will be less than X% complete by EoS
        # (max_growing_periods_after_eos after region's final dek)
        if sos_dic['exclude_areas_less_than'] == "True":
            if sos not in [0, no_start, config.SOS_YET_TO_START]:
                if 0 < ((region_eos + sos_dic['max_growing_periods_after_eos']
                        - sos) / lgp * 100.0) < \
                        float(sos_dic['exclude_areas_less_than_percent']):
                    sos = no_start

    else:
        sos = 0

    return sos


def calc_sos_mode_1d_for_extended(stacked_arr, periodicity, cross_year):
    """
    Calculates the mode for the extended SOS
    """
    if periodicity == 'Dekadal':
        no_start = config.DEKAD_NOSTART
        shift = config.DEKADS_PER_YEAR
        half_year_val = shift/2
    else:
        no_start = config.PENTAD_NOSTART
        shift = config.PENTADS_PER_YEAR
        half_year_val = shift/2

    valid_idx = np.where(np.logical_and(
            np.logical_and(
                stacked_arr != no_start,
                stacked_arr != config.SOS_YET_TO_START),
            stacked_arr != 0))
    if valid_idx[0].size > 0:
        vals, counts = np.unique(stacked_arr[valid_idx], return_counts=True)
        mode_value = np.argwhere(counts == np.max(counts))
        modes_arr = np.array(vals[mode_value].flatten().tolist())
        # If cross_year then add shift to values less than half_year_value
        if cross_year:
            modes_arr = np.where(np.logical_and(modes_arr <= half_year_val,
                                                0 < modes_arr),
                                 modes_arr + shift,
                                 modes_arr)
        mode = np.min(modes_arr)
        if mode > shift:
            mode -= shift
    else:
        # No indexes that do not contain 253, 254, or 255
        yet_to_start = np.where(stacked_arr == config.SOS_YET_TO_START)
        if yet_to_start[0].size > 0:
            # Set as yet to start
            mode = config.SOS_YET_TO_START
        else:
            # If there are 0 yet to starts then see if there are no starts
            no_start_arr = np.where(stacked_arr == no_start)
            if no_start_arr[0].size > 0:
                mode = no_start_arr
            else:
                # only 0 values in array
                mode = 0

    return mode


def calc_avg_med_1d_for_extended(stacked_arr, stat_type, ignore_vals=False):
    """
    Calculates the average or median of the array excluding the
    253, 254, 255 values unless that is all there is
    params(numpy arr) - stacked_arr - the 1-D numpy arr
    params(string) - stat_type - 'Average' or 'Median'
    ignore_vals(bool) - ignore_vals - defaults to false, true ignores
        253, 254, 255 values
    """
    if ignore_vals:
        # Then we ignore the 253, 254, 255 vals in the calculations
        valid_idx = np.where(np.logical_and(
            np.logical_and(
                stacked_arr != config.FAIL_NO_START,
                stacked_arr != config.YET_TO_START),
            stacked_arr != config.NO_DATA))
        if valid_idx[0].size > 0:
            if stat_type == "Average":
                stat_val = np.nanmean(stacked_arr[valid_idx])
            elif stat_type == "Median":
                stat_val = np.nanmedian(stacked_arr[valid_idx])
        else:
            # No indexes that do not contain 253, 254, or 255
            yet_to_start = np.where(stacked_arr == config.YET_TO_START)
            if yet_to_start[0].size > 0:
                # Set as yet to start
                stat_val = config.YET_TO_START
            else:
                # If there are 0 yet to starts then see if there are no starts
                no_start = np.where(stacked_arr == config.FAIL_NO_START)
                if no_start[0].size > 0:
                    stat_val = config.FAIL_NO_START
                else:
                    # only 255 values in array
                    stat_val = config.NO_DATA
    else:
        if stat_type == "Average":
            stat_val = np.nanmean(stacked_arr)
        elif stat_type == "Median":
            stat_val = np.nanmedian(stacked_arr)

    if np.isnan(stat_val):
        return stat_val
    else:
        return round(stat_val)


def calc_masked_pct_of_avg_wrsi(data_array, avg_array, mask_array, nd_val=-9999):
    """
    Calculate the numpy percent of average array, then
    apply the mask.  Formula is result = data * 100 / avg

    Arguments:
    data_array -- 2d numpy geotiff array.
    avg_array -- 2d numpy geotiff array.
    mask_array -- 2d numpy mask array.
    nd_val -- Integer no data value.
    Returns:
    masked_pct_avg_array -- 2d numpy masked percent of average array.
    """
    masked_pct_avg_array = None
    
    if mask_array.shape != data_array.shape != avg_array.shape:
        raise IOError("Mask and/or array shape doesn't match")

    # If the average array is 0 or the data array is nan, set the output to 255
    # otherwise compute the percent of average
    pct_avg_array = np.where(
        np.logical_or(avg_array == 0, np.isnan(data_array)), 255,
        ((data_array * 100.0) / avg_array).round())
    
    masked_pct_avg_array = np.ma.masked_array(
        pct_avg_array, mask=mask_array, fill_value = nd_val).filled()

    return masked_pct_avg_array
