'''
Module Name:  ImageFusion.py

Created on Apr 22, 2014

@author: Bo Yang

All Rights Reserved
'''


#===================================================================================
#              Import Lib & Set Module Parameter
#===================================================================================
import numpy as np
import math
import time
from numpy import matrix, linalg, squeeze, asarray
from os import chdir

#Set the working directory
chdir(r'E:\GoogleDrive\Pyworkspace_Lab\Image_Fusion_v16_GSoc\MODIS_NDVI_Clip')

#Set the raster nodata value
no_data = -3.40282346639e+038

#Set the Spatial t ratio
spatial_temporal_factor=1.0

#set the data range, bosolute values larger than it will be treat as no data
data_valid_range = -5000

#read Exp/gaussian_model fitting parameters from txt file.
parameter_text_file= r'Fitted1.txt'
fit_model_type= 'exponential_model'

#use smaller moving window for the quick fusion mode
bool_quick_mode= 'true'

#If output the acciociate uncertainty
bool_uncertainty= 'false'

#set output raster name
output_fusion_raster= 'output0724'

#===================================================================================
#                   load the parameter from text file
#===================================================================================
parameter_text=open(parameter_text_file,'r')
allLines=parameter_text.readlines()
parameter_text.close()

spatial_parameter_text=str(allLines[1])
temporal_parameter_text=str(allLines[3])

spatial_parameter_array = [float(s) for s in spatial_parameter_text.split(';') if s[0].isdigit()]
temporal_parameter_array = [float(s) for s in temporal_parameter_text.split(';') if s[0].isdigit()]

print "Reading the ", fit_model_type, " fitting model parameters"
print "Spatial fitting parameters are: ", spatial_parameter_array
print "Temporal fitting parameters are: ", temporal_parameter_array
#===================================================================================
#                   input fine Raster
#===================================================================================
#load fine raster to a matrix
fine_raster_text = open(r'lt5ndvi_0716.txt','r')
matrix_fine_raster = np.loadtxt(fine_raster_text.readlines(), skiprows=6)#First six lines are assumed to be image head
fine_raster_text.close()

#===================================================================================
#                    input Coarse Raster  
#===================================================================================
#load all coarse rasters
Coase1 = open(r'a2002203era3_07232015.txt','r')
Coase2 = open(r'A2002204Era2_07232015.txt','r')
Coase3 = open(r'A2002205Era1_07242015.txt','r')
CoaRaster1 = np.loadtxt(Coase1.readlines(), skiprows=6)#First six lines are assumed to be image head
CoaRaster2 = np.loadtxt(Coase2.readlines(), skiprows=6)#First six lines are assumed to be image head
CoaRaster3 = np.loadtxt(Coase3.readlines(), skiprows=6)#First six lines are assumed to be image head
Coase1.close()
Coase2.close()
Coase3.close()

#construct an array of all coarse resolution images
matrix_coarse_raster = [CoaRaster1,CoaRaster2,CoaRaster3]

#===================================================================================
#        calculate the spatio-temporal trend at the last predict day
#===================================================================================
trend_sum = 0.0
trend_count = 0
for j in range(matrix_coarse_raster[-1].shape[0]):
    for i in range(matrix_coarse_raster[-1].shape[1]):
        StackCoarse = []
        if matrix_coarse_raster[-1][j][i] > data_valid_range:
            trend_sum = trend_sum + matrix_coarse_raster[-1][j][i]
            trend_count = trend_count + 1
spatial_temporal_trend = trend_sum/trend_count

print "The ST trend of the last day is: " + str(spatial_temporal_trend)
#===================================================================================
#                       Create intermediate Rasters
#===================================================================================
matrix_new_fine = np.asarray([[0.0 for i in range(matrix_fine_raster.shape[1])] for j in range(matrix_fine_raster.shape[0])])
matrix_uncertainty_raster = np.asarray([[0.0 for i in range(matrix_fine_raster.shape[1])] for j in range(matrix_fine_raster.shape[0])])
matrix_new_coarse = np.asarray([[0.0 for i in range(matrix_fine_raster.shape[1])] for j in range(matrix_fine_raster.shape[0])])

#===================================================================================
#                    construct the moving window
#===================================================================================
moving_window_pixels = []
if str(bool_quick_mode) == 'false':
    windowBegin,windowEnd = 0,27
else:
    windowBegin,windowEnd = 9,18
for m in range(windowBegin,windowEnd):
    for n in range(windowBegin,windowEnd):
        moving_window_pixels.append(m*27+n)

#====================================================================================
#                     Define the spatial and temporal model
#====================================================================================
# Define gaussian_model and exponential_model function
def gaussian_model(t,s,p,x):
    result = t+s*(1-math.exp(-(x**2)/(p**2)))
    return result
def exponential_model(t,s,p,x):
    result = t+ s*(1-math.exp(-x/p))
    return result

# Define Covariance Function
if str(fit_model_type) == 'exponential_model':
    def spatial_exp_covariance(x):
        return (spatial_parameter_array[0] + spatial_parameter_array[1]) - exponential_model(spatial_parameter_array[0], spatial_parameter_array[1], spatial_parameter_array[2], x)
    def temporal_exp_covariance(x):
        return (temporal_parameter_array[0] + temporal_parameter_array[1]) - exponential_model(temporal_parameter_array[0], temporal_parameter_array[1], temporal_parameter_array[2], x)
elif str(fit_model_type) == 'gaussian_model':
    def spatial_exp_covariance(x):
        return (spatial_parameter_array[0] + spatial_parameter_array[1]) - gaussian_model(spatial_parameter_array[0], spatial_parameter_array[1], spatial_parameter_array[2], x)
    def temporal_exp_covariance(x):
        return (temporal_parameter_array[0] + temporal_parameter_array[1]) - gaussian_model(temporal_parameter_array[0], temporal_parameter_array[1], temporal_parameter_array[2], x)


#Spatial distance of fine resolution within fine image
def spatial_distance_fine(i,j):
    return math.sqrt((((j/27)-(i/27))**2)+(((j%27)-(i%27))**2))

#Spatial distance of fine resolution within coarse image
def spatial_distance_coarse(i,j):
    return math.sqrt((((j/3)-(i/3))**2)+(((j%3)-(i%3))**2))

#Spatial distance across fine and coarse resolution image
def spatial_distance_cross(i,j):
    return math.sqrt((((j/27)-((i/3)*9+4))**2)+((j%27)-((i%3)*9+4))**2)


#================================================================================================================
#                                      New matrix delete method
#================================================================================================================
#=================================================================================
#               CY=B   ==>   Y=(C**-1)B   C is the FinalMatrix 
#               B is the vector from predicted location every locations
#               Below is generating B vector
#=================================================================================
#Construct the left C
def left_matrix_C(totalFinePixels=range(729),totalCoaPixels=[range(9),range(9),range(9)]):
    tp=totalFinePixels#total fine pixels in the moving window
    cp=totalCoaPixels#total coarse pixels in the moving window
    matrix_rank=bool(tp)+bool(cp[0] or cp[1] or cp[2])#if all the pixels in moving window are valid value
    matrixSize=len(tp)+len(cp[0])+len(cp[1])+len(cp[2])+matrix_rank
    #print "The size of Water matrix is", len(tp),"+",len(cp[0]),'+',len(cp[1]),'+',len(cp[2]),'+ 2=',matrixSize
    # Create spatial(upper-left) sub-matrix(729*729)
    mUL=[[0 for i in tp] for j in tp]
    for i in tp:
        for j in tp:
            mUL[tp.index(j)][tp.index(i)]= spatial_exp_covariance(spatial_distance_fine(j,i))*temporal_exp_covariance(0)
    # Create Diagonal T1-T7 matrixs (9*9)
    mT11=[[0 for i in cp[0]] for j in cp[0]]
    mT22=[[0 for i in cp[1]] for j in cp[1]]
    mT33=[[0 for i in cp[2]] for j in cp[2]]
    for i in cp[0]:
        for j in cp[0]:
            mT11[cp[0].index(j)][cp[0].index(i)]= spatial_exp_covariance(spatial_distance_coarse(j,i)*9)*temporal_exp_covariance(0)
    for i in cp[1]:
        for j in cp[1]:
            mT22[cp[1].index(j)][cp[1].index(i)]= spatial_exp_covariance(spatial_distance_coarse(j,i)*9)*temporal_exp_covariance(0)
    for i in cp[2]:
        for j in cp[2]:
            mT33[cp[2].index(j)][cp[2].index(i)]= spatial_exp_covariance(spatial_distance_coarse(j,i)*9)*temporal_exp_covariance(0)
    # Create other T13-T31 matrixes (9*9)
    mT12=[[0 for i in cp[1]] for j in cp[0]]
    mT21=[[0 for i in cp[0]] for j in cp[1]] 
    mT13=[[0 for i in cp[2]] for j in cp[0]] 
    mT31=[[0 for i in cp[0]] for j in cp[2]] 
    mT23=[[0 for i in cp[2]] for j in cp[1]]
    mT32=[[0 for i in cp[1]] for j in cp[2]]
    for j in cp[1]:
        for i in cp[2]:
            mT23[cp[1].index(j)][cp[2].index(i)] = spatial_exp_covariance(spatial_distance_coarse(j,i)*9)*temporal_exp_covariance(1)
            mT32[cp[2].index(i)][cp[1].index(j)] = spatial_exp_covariance(spatial_distance_coarse(j,i)*9)*temporal_exp_covariance(1)
    for j in cp[0]:
        for i in cp[1]:
            mT12[cp[0].index(j)][cp[1].index(i)] = spatial_exp_covariance(spatial_distance_coarse(j,i)*9)*temporal_exp_covariance(1)
            mT21[cp[1].index(i)][cp[0].index(j)] = spatial_exp_covariance(spatial_distance_coarse(j,i)*9)*temporal_exp_covariance(1)
    for j in cp[0]:
        for i in cp[2]:      
            mT13[cp[0].index(j)][cp[2].index(i)] = spatial_exp_covariance(spatial_distance_coarse(j,i)*9)*temporal_exp_covariance(2)
            mT31[cp[2].index(i)][cp[0].index(j)] = spatial_exp_covariance(spatial_distance_coarse(j,i)*9)*temporal_exp_covariance(2)
    #Create Upper Right sub-matrix(9*729)
    mUR1=[[0 for i in cp[0]] for j in tp]
    mUR2=[[0 for i in cp[1]] for j in tp]
    mUR3=[[0 for i in cp[2]] for j in tp]
    mLL1=[[0 for i in tp] for j in cp[0]]
    mLL2=[[0 for i in tp] for j in cp[1]]
    mLL3=[[0 for i in tp] for j in cp[2]]
    for j in tp:
        for i in cp[0]:
            mUR1[tp.index(j)][cp[0].index(i)]=spatial_exp_covariance(spatial_distance_cross(i,j)/9.0)*temporal_exp_covariance(0)
            mLL1[cp[0].index(i)][tp.index(j)]=spatial_exp_covariance(spatial_distance_cross(i,j)/9.0)*temporal_exp_covariance(0)
    for j in tp:
        for i in cp[1]:
            mUR2[tp.index(j)][cp[1].index(i)]=spatial_exp_covariance(spatial_distance_cross(i,j)/9.0)*temporal_exp_covariance(1)
            mLL2[cp[1].index(i)][tp.index(j)]=spatial_exp_covariance(spatial_distance_cross(i,j)/9.0)*temporal_exp_covariance(1)
    for j in tp:
        for i in cp[2]:
            mUR3[tp.index(j)][cp[2].index(i)]=spatial_exp_covariance(spatial_distance_cross(i,j)/9.0)*temporal_exp_covariance(2)
            mLL3[cp[2].index(i)][tp.index(j)]=spatial_exp_covariance(spatial_distance_cross(i,j)/9.0)*temporal_exp_covariance(2)
    #Create Lower Left Quasi-diagonal matrix(747*4)
    mDiagLeft=[[0 for i in range(matrixSize-matrix_rank)] for j in range(matrix_rank)] 
    matrix_rankIndex=0
    if len(tp)!=0:
        for i in range(len(tp)):
            mDiagLeft[matrix_rankIndex][i]=1
        matrix_rankIndex=matrix_rankIndex+1
    if len(cp[0])!=0:
        for i in range(len(tp), len(tp)+len(cp[0])):
            mDiagLeft[matrix_rankIndex][i]=1
    if len(cp[1])!=0:
        for i in range(len(tp)+len(cp[0]), len(tp)+len(cp[0])+len(cp[1])):
            mDiagLeft[matrix_rankIndex][i]=1
#    matrix_rankIndex=matrix_rankIndex+1
    if len(cp[2])!=0:
        for i in range(len(tp)+len(cp[0])+len(cp[1]), len(tp)+len(cp[0])+len(cp[1])+len(cp[2])):
            mDiagLeft[matrix_rankIndex][i]=1
    # Create Upper Right Quasi-diagonal matrix(4*747)
    mDiagRight=[[0 for i in range(matrix_rank)] for j in range(matrixSize-matrix_rank)] 
    for j in range(matrixSize-matrix_rank):
        for i in range(matrix_rank):
            mDiagRight[j][i]=mDiagLeft[i][j]
    # Create Lower Right All zero matrix(4*4)
    mZero=[[0 for i in range(matrix_rank)] for j in range(matrix_rank)] 
    #Create The Final Matrix(749*749)
    mFinalMatrix=[[0 for i in range(matrixSize)] for j in range(matrixSize)] 
    for j in range(len(tp)):
        mFinalMatrix[j]=mUL[j]+mUR1[j] + mUR2[j] + mUR3[j] +mDiagRight[j]
    for j in range(len(cp[0])):
        mFinalMatrix[len(tp)+j]=mLL1[j] + mT11[j] + mT12[j] + mT13[j]+ mDiagRight[len(tp)+j]
    for j in range(len(cp[1])):
        mFinalMatrix[len(tp)+len(cp[0])+j]=mLL2[j] + mT21[j] + mT22[j] + mT23[j] + mDiagRight[len(tp)+len(cp[0])+j]
    for j in range(len(cp[2])):
        mFinalMatrix[len(tp)+len(cp[0])+len(cp[1])+j]=mLL3[j] + mT31[j] + mT32[j] +mT33[j] + mDiagRight[len(tp)+len(cp[0])+len(cp[1])+j]
    for j in range(matrix_rank):
        mFinalMatrix[len(tp)+len(cp[0])+len(cp[1])+len(cp[2])+j]=mDiagLeft[j]+mZero[j]
    return matrix(mFinalMatrix)

mleft_matrix_C = left_matrix_C(totalFinePixels=moving_window_pixels)

def right_matrix_B(totalFinePixels=range(729),totalCoaPixels=[range(9),range(9),range(9)],PixelInCenter=364):
    tp=totalFinePixels
    cp=totalCoaPixels
    matrix_rank=bool(tp)+bool(cp[0] or cp[1] or cp[2])
    matrixSize=len(tp)+len(cp[0])+len(cp[1])+len(cp[2])+matrix_rank
    matrixB = [0 for i in range(matrixSize)]
    matrixB[matrixSize-1]=1
    for i in tp:
        matrixB[tp.index(i)]=spatial_exp_covariance(spatial_distance_fine(i,364))*temporal_exp_covariance(2)
    for i in cp[0]:
        matrixB[len(tp)+cp[0].index(i)]=spatial_exp_covariance(spatial_distance_cross(i,PixelInCenter))*temporal_exp_covariance(2)
    for i in cp[1]:
        matrixB[len(tp)+len(cp[0])+cp[1].index(i)]=spatial_exp_covariance(spatial_distance_cross(i,PixelInCenter))*temporal_exp_covariance(1)
    for i in cp[2]:
        matrixB[len(tp)+len(cp[0])+len(cp[1])+cp[2].index(i)]=spatial_exp_covariance(spatial_distance_cross(i,PixelInCenter))*temporal_exp_covariance(0)
    return matrix(matrixB)

mright_matrix_B = [right_matrix_B(totalFinePixels=moving_window_pixels,totalCoaPixels=[range(9),range(9),range(9)],PixelInCenter = i) for i in range(729)]

def solve_del_matrix(fineDelete =np.array([]), coarseDelete = np.array([[],[],[]]),PixelInCenter=364):
    deleteNum = np.concatenate((fineDelete,coarseDelete[0]+len(moving_window_pixels),coarseDelete[1] + (9-len(coarseDelete[0]))+len(moving_window_pixels),coarseDelete[2] + (9-len(coarseDelete[1])) + (9-len(coarseDelete[0])+len(moving_window_pixels))))
    mleftMatrixDel = np.delete(np.delete(mleft_matrix_C,deleteNum,0),deleteNum,1)
    mrightMatrixDel = np.delete(mright_matrix_B[PixelInCenter],deleteNum,1)
    #delete the second last zero item
    mmleftMatrixDel = np.delete(np.delete(mleftMatrixDel,-2,0),-2,1)
    mmrightMatrixDel = np.delete(mrightMatrixDel,-2,1)
    return squeeze(asarray(linalg.solve(matrix(mmleftMatrixDel), matrix(mmrightMatrixDel).T)))

print ("Pre-calculating weight matrix")
full_pixel_weights = [solve_del_matrix(np.array([]),np.array([[],[],[]]),PixelInCenter = i) for i in range(729)]

#==========================================================================================
#                    fusion without covariable
#==========================================================================================
def fusion_no_covariable(Raster=matrix_fine_raster,coaRasters=matrix_coarse_raster):
    bad_pixel_count = 0
    for j in range(15,Raster.shape[0]-15):
        if (j % (Raster.shape[0]/100) == 0):
            print (str(int(j*1.0/Raster.shape[0]*100)) + "%, Row Number: "+ str(j) +", "+ str(round(bad_pixel_count*100.0/(j*Raster.shape[1]),2)) + "% missing pixels now")
        for i in range(15,Raster.shape[1]-15):
            #validPixelsValue = np.array([]) #the input data assigned the fine resolution as null
            validPixelsValue = Raster[j-13+windowBegin:j-13+windowEnd,i-13+windowBegin:i-13+windowEnd].flatten()
            validCoaPixelsValue = np.array([coaRasters[k][j/9-1:j/9-1+3,i/9-1:i/9-1+3].flatten() for k in range(3)])
            validPixels = np.array([]).tolist()
            valid_Pixels = np.where(validPixelsValue<500)[0]
            validCoaPixels0 = np.where(validCoaPixelsValue[0] > data_valid_range)[0].tolist()
            validCoaPixels1 = np.where(validCoaPixelsValue[1] > data_valid_range)[0].tolist()
            validCoaPixels2 = np.where(validCoaPixelsValue[2] > data_valid_range)[0].tolist()
            if len(validCoaPixels0)+len(validCoaPixels1)+len(validCoaPixels2)<3:
                continue
            valid_CoaPixels0 = np.where(validCoaPixelsValue[0] < data_valid_range)[0]
            valid_CoaPixels1 = np.where(validCoaPixelsValue[1] < data_valid_range)[0]
            valid_CoaPixels2 = np.where(validCoaPixelsValue[2] < data_valid_range)[0]
            returnValue=solve_del_matrix(valid_Pixels,[valid_CoaPixels0,valid_CoaPixels1,valid_CoaPixels2],((j%9+9)*27+(i%9+9)))#change here if want to apply changingsupport to uncertainty
            pixelValue32=np.float32(np.concatenate((validCoaPixelsValue[0][validCoaPixelsValue[0]>data_valid_range], validCoaPixelsValue[1][validCoaPixelsValue[1]>data_valid_range], validCoaPixelsValue[2][validCoaPixelsValue[2]>data_valid_range])))
            coefficient32=np.float32(returnValue[:len(validPixels)+len(validCoaPixels0)+len(validCoaPixels1)+len(validCoaPixels2)])
            bad_pixel_count+=1
            matrix_new_fine[j][i] = VectorMultiple(pixelValue32,coefficient32).sum()

#==========================================================================================
#                             fusion with co-variable
#==========================================================================================
def fusion_with_covariable(Raster=matrix_fine_raster,coaRasters = matrix_coarse_raster):
    bad_pixel_count = 0
    for j in range(15,Raster.shape[0]-15):
        if (j % (Raster.shape[0]/100) == 0):
            print (str(int(j*1.0/Raster.shape[0]*100)) + "%, Row Number: "+ str(j) +", "+ str(round(bad_pixel_count*100.0/(j*Raster.shape[1]),2)) + "% missing pixels now")
        for i in range(15,Raster.shape[1]-15):
            validPixelsValue = Raster[j-13+windowBegin:j-13+windowEnd,i-13+windowBegin:i-13+windowEnd].flatten()
            validCoaPixelsValue = np.array([coaRasters[k][j/9-1:j/9-1+3,i/9-1:i/9-1+3].flatten() for k in range(3)])
            if validPixelsValue[validPixelsValue<data_valid_range].any() or validCoaPixelsValue[validCoaPixelsValue<data_valid_range].any():
                validPixels = np.where(validPixelsValue > data_valid_range)[0].tolist()
                valid_Pixels = np.where(validPixelsValue < data_valid_range)[0]
                validCoaPixels0 = np.where(validCoaPixelsValue[0] > data_valid_range)[0].tolist()
                validCoaPixels1 = np.where(validCoaPixelsValue[1] > data_valid_range)[0].tolist()
                validCoaPixels2 = np.where(validCoaPixelsValue[2] > data_valid_range)[0].tolist()
                if len(validCoaPixels0)+len(validCoaPixels1)+len(validCoaPixels2)<3:
                    continue  
                valid_CoaPixels0 = np.where(validCoaPixelsValue[0] < data_valid_range)[0]
                valid_CoaPixels1 = np.where(validCoaPixelsValue[1] < data_valid_range)[0]
                valid_CoaPixels2 = np.where(validCoaPixelsValue[2] < data_valid_range)[0]
                returnValue= solve_del_matrix(valid_Pixels,[valid_CoaPixels0,valid_CoaPixels1,valid_CoaPixels2],((j%9+9)*27+(i%9+9)))#change here if want to apply changingsupport to uncertainty
                pixelValue32=np.float32(np.concatenate((validPixelsValue[validPixelsValue>data_valid_range], validCoaPixelsValue[0][validCoaPixelsValue[0]>data_valid_range], validCoaPixelsValue[1][validCoaPixelsValue[1]>data_valid_range], validCoaPixelsValue[2][validCoaPixelsValue[2]>data_valid_range])))
                coefficient32=np.float32(returnValue[:len(validPixels)+len(validCoaPixels0)+len(validCoaPixels1)+len(validCoaPixels2)])
                bad_pixel_count+=1
            else:
                returnValue = full_pixel_weights[((j%9+9)*27+(i%9+9))]
                pixelValue32 = np.float32(np.concatenate((validPixelsValue, validCoaPixelsValue[0], validCoaPixelsValue[1], validCoaPixelsValue[2])))
                coefficient32 = np.float32(returnValue[:len(validPixelsValue)+len(validCoaPixelsValue[0])+len(validCoaPixelsValue[1])+len(validCoaPixelsValue[2])])
            matrix_new_fine[j][i]=VectorMultiple(pixelValue32,coefficient32).sum()


#===================================================================================
#                 CUDA model can be injected here
#=================================================================================== 
def VectorMultiple(pixelValue,coefficient):
    return pixelValue * coefficient

#===================================================================================
#                    save raster
#===================================================================================
def save_raster(Raster = matrix_new_fine, Name = output_fusion_raster):
#save raster
    np.savetxt(Name,Raster,fmt = '%1.4f')

#===================================================================================
#                Images detrend (trend is the first day mean)
#===================================================================================

def fine_raster_detrend(Raster=matrix_fine_raster, trend = spatial_temporal_trend):
    for j in range(Raster.shape[0]):
        for i in range(Raster.shape[1]):
            if Raster[j][i] > data_valid_range:
                Raster[j][i] = Raster[j][i] - trend

def coarse_raster_detrend(Raster=matrix_coarse_raster,trend=spatial_temporal_trend):
    for j in range(matrix_coarse_raster[0].shape[0]):
        for i in range(matrix_coarse_raster[0].shape[1]):
            for k in range(len(matrix_coarse_raster)):
                if Raster[k][j][i] > data_valid_range:
                    Raster[k][j][i] = Raster[k][j][i]-trend

#======================================================================================
#               Run the fusion program over the fine raster
#======================================================================================

def final_run():
    print (time.strftime("%m-%d %X",time.localtime())+" Detrending the coarse images")
    coarse_raster_detrend(matrix_coarse_raster)  
    #
    print (time.strftime("%m-%d %X",time.localtime())+" Detrending the fine image")
    fine_raster_detrend(matrix_fine_raster)
    #
    print (time.strftime("%m-%d %X",time.localtime())+" Fusing the images")
    fusion_no_covariable()
    #
    print (time.strftime("%m-%d %X",time.localtime())+" Adding the trend to the residual image")
    fine_raster_detrend(matrix_new_fine, -(spatial_temporal_trend))
    #
    print (time.strftime("%m-%d %X",time.localtime())+" Saving the result")
    save_raster()
    #
    print (time.strftime("%m-%d %X",time.localtime())+output_fusion_raster+ " has finished")


final_run()