#Created by Theodore Endreny and his team at SUNY ESF for work with the i-Tree Cool Air model, te@esf.edu
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

#Model: p_itca_03_prepModel_Outputs_Mosaic.py

#Goal of script:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Mosaic together separate time-based map outputs in quad### folders generated by i-Tree HydroPlus CoolAir simulations

#How to use script:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 1a) If p_itca_02_PrepModel_Inputs.py generated quad### folders of inputs for area of interest (AOI), then ...
# 1b) ... have run HydroPlus CoolAir to generate Tair_K_YYYYMMDDHH.asc output for each quad### in AOI
# 2a) Change below the several list, variable, path, and file names to your project needs in the section below
# 2b) Note: city_name, grid_resolution_list should contain values used in HydroPlus simulations
# 2c) Note: variables RAM_max_cols and RAM_max_rows should contain values used in HydroPlus simulations
# 2d) Note: list or variable dates used for time-based outputs should contain values used in HydroPlus simulations
# 3) Use Python version 2.7 to launch this script
# 4) Find HydroPlus CoolAir mosaicked outputs in the variable path_prepModel_trunk_fs, defined below, outside of any quad### folder
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

#Define Python libraries and extensions# import numpy to use numpy arrays
import arcpy
from arcpy.sa import *
from arcpy.analysis import *
import os
import shutil
#Overwrite existing datasets option is enabled
from arcpy import env
env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Please update the following list, variable, path, and file names to your project needs
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Change: city abbreviation; first 3 letters of city name + "_" + 2 letters of state code (e.g. 'syr_ny' for Syracuse, New York)
city_list = ['fic_or'] #'was_dc' #'ric_va' #'bal_md' #'04001' #'41051' #['was_dc', 'bal_md', 'ric_va']
Central_meridian_ID = -96

#Flag_OverwriteOutput used to manage existing output; 0 = will not overwrite existing output; 1 = overwrite existing output
#Note: If Flag_OverwriteOutput = 0, then prior files should be manually deleted from path_prepData_trunk_fs + city_name + ires
Flag_OverwriteOutput = 1

#Change: data grid resolution for original NHD DEM and NLCD data
grid_resolution_base_m = 30
#Change: data cell size (the length of one edge of a square cell in meters)
grid_resolution_list = [grid_resolution_base_m] #Options include: [grid_resolution_base_m, 90, 180, 300]

#Change optional: Maximum number of columns and rows used in i-Tree HydroPlus CoolAir model simulation; affects RAM
RAM_max_cols = 150 #120 310
RAM_max_rows = 150 #120 310

#flag_delete_quad_tif_files = 1 will delete temporary files used in quad folders holding sub-AOI areas, = 0 will not delete
flag_delete_quad_tif_files = 1
#Change: flag_time_based_output_range = 1 to have time-based output for continuous range, or flag_time_based_output_range = 0 to have discrete periods
#cte 20211028 Construct list of continuous range
flag_time_based_output_range = 0
#Change: flag_set_Parameters_by_city_name, = 0 = simulation parameters set as default; 1 = simulation parameters based on city name conditional, useful for multiple areas
flag_set_Parameters_by_city_name = 1

#Change: If flag_time_based_output_range = 1, Start date and hour of continuous range for simulation time-based output map date and hour (YYYYMMDDHH), string
output_map_date_hour_start = '2021062913' #was_dc = 2018082800 #bal_md = 2018082900 #ric_va = 2017071300
#Change: If flag_time_based_output_range = 1, Start date and hour of continuous range for simulation time-based output map date and hour (YYYYMMDDHH), string
output_map_date_hour_stop = '2021062917' #was_dc = 2018082823 #bal_md = 2018082923 #ric_va = 2017071323
#Change: If flag_time_based_output_range = 0, List of discrete periods for simulation time-based output map date and hour (YYYYMMDDHH), strings
output_map_date_hour_list = ['2021062917'] #was_dc '2018082806','2018082815','2018082819' #bal_md '2018082906','2018082915','2018082919' #ric_va '2017071306','2017071315','2017071319'

#Change optional: Name map_mosaic folder to hold ArcGIS .tif outputs
folder_map_output_mosaic = 'map_mosaic'
#Change optional: Variable options are TAir_K or other specific outputs for time-based output (AbsHumidity_kg_p_m3, ..., WaterNR_W_p_m2) 
quad_map_type = 'TAir_K' #'SlopeGround_rad.asc' #'AspectGround_N_0_rad.asc' #'TAir_weighted_K' #'TAir_K'

#Change: Trunk path to NLCD and DEM source data, using distinct forward slash (fs), which must contain file lc_%s_%s.tif %(city_name,ires)
path_prepData_trunk_fs = 'C:/iTree/HydroPlus/TestingFilesAndScript/TestCases/CoolAir/MultipleStations/gis_data/' 
#Change: Target path to i-Tree HydroPlus CoolAir simulations, using distinct forward slash (fs) 
path_prepModel_trunk_fs = 'C:/iTree/HydroPlus/TestingFilesAndScript/TestCases/CoolAir/MultipleStations/' 
#Change optional: Folder name appended to 3 digit folder number within the simulation output folder
folder_quad_ID_str01 = 'quad' 

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Loop through all city_name in city_list
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for city_name in city_list:

    #if flag_set_Parameters_by_city_name equals 1 then
    if flag_set_Parameters_by_city_name == 1:
        if city_name == 'fic_or':
            #Change: If flag_time_based_output_range = 1, Start date and hour of continuous range for simulation time-based output map date and hour (YYYYMMDDHH), string
            output_map_date_hour_start = '2021062913'
            #Change: If flag_time_based_output_range = 1, Start date and hour of continuous range for simulation time-based output map date and hour (YYYYMMDDHH), string
            output_map_date_hour_stop = '2021062917'
            #Change: If flag_time_based_output_range = 0, List of discrete periods for simulation time-based output map date and hour (YYYYMMDDHH), strings
            output_map_date_hour_list = ['2021062917']

    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #Loop through all ires in grid_resolution_list
    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    for ires in grid_resolution_list:

        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #Step: Establish sub-AOI quads if needed fit row and column limits for RAM
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #create test to see if area of study has been sub-divided into quads of area of interest
        path_aoi = path_prepData_trunk_fs + '%s/' %(city_name)
        path_aoi_resolution = path_prepData_trunk_fs + '%s/res_%s/' %(city_name, ires)
     
        if not os.path.exists(path_aoi):
            print('Warning: The following directory is needed for the shapefile but does not exist.\n{}\n'.format(path_aoi))
        if not os.path.exists(path_aoi_resolution):
            print('Warning: The following directory is needed for the raster inputs but does not exist.\n{}\n'.format(path_aoi_resolution))
        
        constant_value = "1"
        data_type = "INTEGER"
        match_raster = path_aoi_resolution + 'lc_%s_%s.tif' %(city_name, ires)
        cell_size = match_raster
        extent = match_raster
        in_raster_check = CreateConstantRaster(constant_value, data_type, cell_size, extent)
        #get extent of raster
        in_raster_bottom = in_raster_check.extent.YMin
        in_raster_top = in_raster_check.extent.YMax
        in_raster_left = in_raster_check.extent.XMin
        in_raster_right = in_raster_check.extent.XMax
        
        #quads_intervals_cols and quads_intervals_rows are number of quadrants for the map grid to comply with RAM_max_cols and RAM_max_rows 
        quads_intervals_cols = float(in_raster_right - in_raster_left) / float(ires * RAM_max_cols)
        quads_intervals_rows = float(in_raster_top - in_raster_bottom) / float(ires * RAM_max_rows)

        #create empty list of quads_folders
        quad_ID = [] 
        quad_ID_value = 0
        #add a value of 1 to quads_intervals_cols and quads_intervals_rows to account for any remainder, remainder found using modulus %
        if (quads_intervals_cols % 1) > 0 or quads_intervals_cols == 0: quads_intervals_cols = int(quads_intervals_cols) + 1
        if (quads_intervals_rows % 1) > 0 or quads_intervals_rows == 0: quads_intervals_rows = int(quads_intervals_rows) + 1
        #add a value of 1 to quads_intervals_cols and quads_intervals_rows to account for any remainder
        quad_ID_last = quads_intervals_cols * quads_intervals_rows

        #create quad_ID list by appending quad numbers list until maximum folder is reached
        while(quad_ID_value < quad_ID_last): 
            quad_ID_value += 1
            quad_ID.append(quad_ID_value) 

        print('\n')
        if (quad_ID_value == 1):
            print('1. AOI at resolution {} is covered as a single AOI, without sub-dividing in quad folders.'.format(ires))
        else:
            print('1. AOI at resolution {} is sub-divided and covered by {} quad folders due to RAM map extent limits.'.format(ires, quad_ID_value))

        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #Step: Define paths and folders for i-Tree CoolAir model outputs that will mosaic 
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #path containing i-Tree CoolAir model .tif mosaic
        path_prepModel_mosaic = path_prepModel_trunk_fs + '%s/Outputs%s/' %(city_name, ires)

        #path containing new .tif maps; if no quad### folders this is just identifying one AOI folder
        path_map_mosaic_output = os.path.join(path_prepModel_mosaic, folder_map_output_mosaic)

        #query_folder_exists is result of Python os query if the path exists, either true or false
        query_folder_exists = os.path.exists(path_map_mosaic_output)
        #If the file exists, do not create or destroy
        if query_folder_exists == True:
            print("1a: {} folder exists, great!".format(path_map_mosaic_output))
        #If the file does not exist, create
        elif query_folder_exists == False:
            print("1a: Creating {} folder.".format(path_map_mosaic_output))
            os.makedirs(path_map_mosaic_output)               
        
        #Loop through all iquad_count in quad_ID list
        for iquad_count in quad_ID:
        
            #define strings: condition when >1 map is processed, when no quads were generated above
            if (len(quad_ID) > 1):
                if iquad_count < 10000: iquad_count_str = '%s' %(iquad_count)
                if iquad_count < 1000: iquad_count_str = '0%s' %(iquad_count)
                if iquad_count < 100: iquad_count_str = '00%s' %(iquad_count)
                if iquad_count < 10: iquad_count_str = '000%s' %(iquad_count)
                quad_ID_str01 = folder_quad_ID_str01 + iquad_count_str
                folder_str_fs = '/'
            #define strings: condition when 1 map is processed, when no quads were generated above
            else:
                iquad_count_str = '%s' %(iquad_count)
                quad_ID_str01 = ''
                folder_str_fs = ''            

            #If file_iTCA_output is TAir_K or TAir_weighted_K then add hour to string
            if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                #path containing i-Tree CoolAir model ASCII outputs, if no quad### folders this is just identifying one AOI folder
                path_prepModel_outputs = path_prepModel_trunk_fs + '%s/Outputs%s/%s%s' %(city_name, ires, quad_ID_str01, folder_str_fs)
            #Else if file_iTCA_output is AspectGround_N_0_rad.dat then do not add hour 
            elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                #path containing i-Tree CoolAir model ASCII outputs, if no quad### folders this is just identifying one AOI folder
                path_prepModel_outputs = path_prepModel_trunk_fs + '%s/Inputs%s/%s%s' %(city_name, ires, quad_ID_str01, folder_str_fs)
                    
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            #Step: Create Albers_Conical_Equal_Area projected .tif Raster from i-Tree CoolAir model outputs in ASCII format
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            #Loop through i-Tree CoolAir model output list 
            for date_hour_str in output_map_date_hour_list:
                
                #Define file_iTCA_output as ASCII .asc output map
                #If file_iTCA_output is TAir_K or TAir_weighted_K hen add hour to string
                if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                    file_iTCA_output = quad_map_type + '_' + date_hour_str + '.asc'
                #Else if file_iTCA_output is AspectGround_N_0_rad.asc then add .asc to end
                elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                    file_iTCA_output = quad_map_type + '.asc'
                    #os.rename takes argument of old_name and new_name, and converts
                    old_file = path_prepModel_outputs + quad_map_type
                    new_file = path_prepModel_outputs + file_iTCA_output
                    print("old_file is {}".format(old_file))
                    shutil.copyfile(old_file, new_file)

                #convert i-Tree CoolAir model output ASCII to ArcGIS raster
                #input ASCII file to convert to raster, from i-Tree CoolAir model output
                in_ASCII_file = path_prepModel_outputs + file_iTCA_output
                
                #If file_iTCA_output is TAir_K then add hour to string
                if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                    #path to time-based output file as .tif for each quad###, named by quad_map_type, date_hour_str, city_name, ires, and iquad_count_str
                    out_raster = path_map_mosaic_output + '/%s_%s_%s_%s_%s.tif' %(quad_map_type, date_hour_str, city_name, ires, iquad_count_str)
                #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour etc. 
                elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                    #path to time-based output file as .tif for each quad###, named by quad_map_type, city_name, ires, and iquad_count_str
                    out_raster = path_map_mosaic_output + '/%s_%s_%s_%s.tif' %(quad_map_type, city_name, ires, iquad_count_str)

                #Flag_OverwriteOutput used to manage existing output; 0 = will not overwrite existing output; 1 = overwrite existing output
                #If Flag_OverwriteOutput equals 1 then create files regardless if they exist
                if (Flag_OverwriteOutput == 1):
                    #If file_iTCA_output is TAir_K then add hour to string
                    if (quad_map_type == 'TAir_K'):
                        print("2a: Creating {}_{}_{}_{}_{}.tif file.".format(quad_map_type, date_hour_str, city_name, ires, iquad_count_str))
                    #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                    elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                        print("2a: Creating {}_{}_{}_{}.tif file.".format(quad_map_type, city_name, ires, iquad_count_str))
                        
                    #air temperatures are float decimal values
                    rasterType = 'FLOAT'
                    #arcpy function ASCIIToRaster_conversion with inputs as defined by ESRI
                    arcpy.ASCIIToRaster_conversion(in_ASCII_file, out_raster, rasterType)
                
                    #If file_iTCA_output is TAir_K then add hour to string
                    if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                        print("2b: Projecting to Albers_Conical_Equal_Area {}_{}_{}_{}_{}.tif file.".format(quad_map_type, date_hour_str, city_name, ires, iquad_count_str))
                    #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                    elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                        print("2b: Projecting to Albers_Conical_Equal_Area {}_{}_{}_{}.tif file.".format(quad_map_type, city_name, ires, iquad_count_str))

                    #define a coordinate system for the new raster
                    #the input raster is the output raster from the previous step
                    in_raster = out_raster
                    #project to the Albers_Conical_Equal_Area with variables Central_meridian_ID
                    coor_system = "PROJCS['Albers_Conical_Equal_Area',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['false_easting',0.0],PARAMETER['false_northing',0.0],PARAMETER['Central_Meridian',%s.0],PARAMETER['standard_parallel_1',29.5],PARAMETER['standard_parallel_2',45.5],PARAMETER['latitude_of_origin',23.0],UNIT['Meter',1.0]]" % (Central_meridian_ID)
                    #arcpy function DefineProjection_management with inputs as defined by ESRI
                    arcpy.DefineProjection_management(in_raster, coor_system)

                #Else If Flag_OverwriteOutput not equal to 1 then do not create files if they exist
                else:
                    #If the file exists, do not create or destroy
                    if arcpy.Exists(out_raster):
                        #If file_iTCA_output is TAir_K then add hour to string
                        if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                            print("2a: Existing {}_{}_{}_{}_{}.tif".format(quad_map_type, date_hour_str, city_name, ires, iquad_count_str))
                        #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                        elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                            print("2a: Existing {}_{}_{}_{}.tif".format(quad_map_type, city_name, ires, iquad_count_str))
                    #If the file does not exist, create
                    else:
                        #If file_iTCA_output is TAir_K then add hour to string
                        if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                            print("2a: Creating {}_{}_{}_{}_{}.tif file.".format(quad_map_type, date_hour_str, city_name, ires, iquad_count_str))
                        #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                        elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                            print("2a: Creating {}_{}_{}_{}.tif file.".format(quad_map_type, city_name, ires, iquad_count_str))
                            
                        #air temperatures are float decimal values
                        rasterType = 'FLOAT'
                        #arcpy function ASCIIToRaster_conversion with inputs as defined by ESRI
                        arcpy.ASCIIToRaster_conversion(in_ASCII_file, out_raster, rasterType)
                    
                        #If file_iTCA_output is TAir_K then add hour to string
                        if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                            print("2b: Projecting to Albers_Conical_Equal_Area {}_{}_{}_{}_{}.tif file.".format(quad_map_type, date_hour_str, city_name, ires, iquad_count_str))
                        #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                        elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                            print("2b: Projecting to Albers_Conical_Equal_Area {}_{}_{}_{}.tif file.".format(quad_map_type, city_name, ires, iquad_count_str))

                        #define a coordinate system for the new raster
                        #the input raster is the output raster from the previous step
                        in_raster = out_raster
                        #project to the  Albers_Conical_Equal_Area with variables Central_meridian_ID
                        coor_system = "PROJCS['Albers_Conical_Equal_Area',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['false_easting',0.0],PARAMETER['false_northing',0.0],PARAMETER['Central_Meridian',%s.0],PARAMETER['standard_parallel_1',29.5],PARAMETER['standard_parallel_2',45.5],PARAMETER['latitude_of_origin',23.0],UNIT['Meter',1.0]]" % (Central_meridian_ID)
                        #arcpy function DefineProjection_management with inputs as defined by ESRI
                        arcpy.DefineProjection_management(in_raster, coor_system)

            #close loop date_hour_str
        #close loop iquad_count in quad_ID

        #if quad### folders exist then mosaic time-based .tif files generated above
        if (len(quad_ID) > 1):

            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            #Step: Create single mosaic of AOI for each date_hour_str using Albers_Conical_Equal_Area projected .tif files
            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            #loop through date_hour_str in output_map_date_hour_list 
            for date_hour_str in output_map_date_hour_list:
                
                #reset quad rasters that will be assembled in mosaic
                m_iTCA_raster_list = ''
                #initialize counter
                i_cnt_quad = 0
                
                #Loop through all iquad_count in quad_ID list to create list of file names to be placed into mosaic
                for iquad_count in quad_ID:
                    
                    #construct quad string name
                    if iquad_count < 10000: iquad_count_str = '%s' %(iquad_count)
                    if iquad_count < 1000: iquad_count_str = '0%s' %(iquad_count)
                    if iquad_count < 100: iquad_count_str = '00%s' %(iquad_count)
                    if iquad_count < 10: iquad_count_str = '000%s' %(iquad_count)
                    quad_ID_str01 = folder_quad_ID_str01 + iquad_count_str
                    #Define forward slash used to construct paths
                    quad_ID_str_fs = '/'

                    #If file_iTCA_output is TAir_K then add hour to string
                    if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                        #path to time-based output file as .tif for each quad###, named by quad_map_type, date_hour_str, city_name, ires, and iquad_count_str
                        in_raster_quad = path_map_mosaic_output + '/%s_%s_%s_%s_%s.tif' %(quad_map_type, date_hour_str, city_name, ires, iquad_count_str)
                    #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                    elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                        #path to time-based output file as .tif for each quad###, named by quad_map_type, city_name, ires, and iquad_count_str
                        in_raster_quad = path_map_mosaic_output + '/%s_%s_%s_%s.tif' %(quad_map_type, city_name, ires, iquad_count_str)


                    #if iquad_count is not last iquad_count in quad_ID list
                    if (i_cnt_quad < len(quad_ID)):
                        #append list of file names going into ArcGIS mosaic with semicolon
                        m_iTCA_raster_list = m_iTCA_raster_list + in_raster_quad + ';'
                    #if iquad_count is last iquad_count in quad_ID list
                    elif (i_cnt_quad == len(quad_ID)):
                        #append list of file names going into ArcGIS mosaic without semicolon
                        m_iTCA_raster_list = m_iTCA_raster_list + in_raster_quad

                    #advance counter
                    i_cnt_quad = i_cnt_quad + 1
        
                #final mosaic output path to i-Tree CoolAir model output in Raster folder as .tif file
                #path to time-based output file as .tif for AOI, as temporary file that will be reconfigured by select by mask
                    #If file_iTCA_output is TAir_K then add hour to string
                    if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                        file_out_raster_mosaic_tmp = '/%s_%s_%s_%s_tmp.tif' %(quad_map_type, date_hour_str, city_name, ires)
                    #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                    elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                        file_out_raster_mosaic_tmp = '/%s_%s_%s_tmp.tif' %(quad_map_type, city_name, ires)

                #Flag_OverwriteOutput used to manage existing output; 0 = will not overwrite existing output; 1 = overwrite existing output
                #If Flag_OverwriteOutput equals 1 then create files regardless if they exist
                if (Flag_OverwriteOutput == 1):
                    #If file_iTCA_output is TAir_K then add hour to string
                    if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                        print("3a: Creating {}_{}_{}_{}_tmp.tif file.".format(quad_map_type, date_hour_str, city_name, ires))
                    #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                    elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                        print("3a: Creating {}_{}_{}_tmp.tif file.".format(quad_map_type, city_name, ires))
                        
                    #keep the same geographic coordinate system and cell size the same
                    coord_system = ""
                    cellsize= ""
                    pixel_type = "32_BIT_FLOAT"
                    number_of_bands = "1"
                    #default settings, input rasters should not overlap
                    mosaic_method = "LAST"
                    mosaic_colormap_mode = "FIRST"
                    #arcpy function MosaicToNewRaster_management with inputs as defined by ESRI
                    arcpy.MosaicToNewRaster_management(m_iTCA_raster_list, path_map_mosaic_output, file_out_raster_mosaic_tmp, coord_system, pixel_type, cellsize, number_of_bands, mosaic_method, mosaic_colormap_mode)

                #Else If Flag_OverwriteOutput not equal to 1 then do not create files if they exist
                else:
                    #If the file exists, do not create or destroy
                    if arcpy.Exists(path_map_mosaic_output + file_out_raster_mosaic_tmp):
                        #If file_iTCA_output is TAir_K then add hour to string
                        if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                            print("3a: Existing {}_{}_{}_{}_tmp.tif file".format(quad_map_type, date_hour_str, city_name, ires))
                        #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                        elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                            print("3a: Existing {}_{}_{}_tmp.tif".format(quad_map_type, city_name, ires))
                    #If the file does not exist, create
                    else:
                        #If file_iTCA_output is TAir_K then add hour to string
                        if (quad_map_type == 'TAir_K'):
                            print("3a: Creating {}_{}_{}_{}_tmp.tif file.".format(quad_map_type, date_hour_str, city_name, ires))
                        #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                        elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                            print("3a: Creating {}_{}_{}_tmp.tif file.".format(quad_map_type, city_name, ires))
                            
                        #keep the same geographic coordinate system and cell size the same
                        coord_system = ""
                        cellsize= ""
                        pixel_type = "32_BIT_FLOAT"
                        number_of_bands = "1"
                        #default settings, input rasters should not overlap
                        mosaic_method = "LAST"
                        mosaic_colormap_mode = "FIRST"
                        #arcpy function MosaicToNewRaster_management with inputs as defined by ESRI
                        arcpy.MosaicToNewRaster_management(m_iTCA_raster_list, path_map_mosaic_output, file_out_raster_mosaic_tmp, coord_system, pixel_type, cellsize, number_of_bands, mosaic_method, mosaic_colormap_mode)
            #close loop of date_hour_str
        #close loop of iquad_count in quad_ID
        #if no quad### folders exists then alternative to mosaic for time-based .tif files 
        else:
            #loop through date_hour_str in output_map_date_hour_list 
            for date_hour_str in output_map_date_hour_list:
                #If file_iTCA_output is TAir_K then add hour to string
                if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                    file_in_raster_quad = path_map_mosaic_output + '/%s_%s_%s_%s_%s.tif' %(quad_map_type, date_hour_str, city_name, ires, iquad_count_str)
                #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                    file_in_raster_quad = path_map_mosaic_output + '/%s_%s_%s_%s.tif' %(quad_map_type, city_name, ires, iquad_count_str)
                    
                #final mosaic output path to i-Tree CoolAir model output in Raster folder as .tif file
                #path to time-based output file as .tif for AOI, as temporary file that will be reconfigured by select by mask
                #If file_iTCA_output is TAir_K then add hour to string
                if (quad_map_type == 'TAir_K'):
                    file_out_raster_name_tmp = '/%s_%s_%s_%s_tmp.tif' %(quad_map_type, date_hour_str, city_name, ires)
                #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):                
                    file_out_raster_name_tmp = '/%s_%s_%s_tmp.tif' %(quad_map_type, city_name, ires)

                file_out_raster_tmp = path_map_mosaic_output + file_out_raster_name_tmp 
                    
                #Flag_OverwriteOutput used to manage existing output; 0 = will not overwrite existing output; 1 = overwrite existing output
                #If Flag_OverwriteOutput equals 1 then create files regardless if they exist
                if (Flag_OverwriteOutput == 1):
                    #If file_iTCA_output is TAir_K then add hour to string
                    if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                        print("3a: Creating {}_{}_{}_{}_tmp.tif file.".format(quad_map_type, date_hour_str, city_name, ires))
                    #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                    elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                        print("3a: Creating {}_{}_{}_tmp.tif file.".format(quad_map_type, city_name, ires))
                    #arcpy function management.Rename with inputs as defined by ESRI
                    arcpy.management.Rename(file_in_raster_quad, file_out_raster_tmp)

                #Else If Flag_OverwriteOutput not equal to 1 then do not create files if they exist
                else:
                    #If the file exists, do not create or destroy
                    if arcpy.Exists(path_map_mosaic_output + file_out_raster_name_tmp):
                        #If file_iTCA_output is TAir_K then add hour to string
                        if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                            print("3a: Existing {}_{}_{}_{}_tmp.tif".format(quad_map_type, date_hour_str, city_name, ires))
                        #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                        elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                            print("3a: Existing {}_{}_{}_tmp.tif".format(quad_map_type, city_name, ires))
                    #If the file does not exist, create
                    else:
                        #If file_iTCA_output is TAir_K then add hour to string
                        if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                            print("3a: Creating {}_{}_{}_{}_tmp.tif file.".format(quad_map_type, date_hour_str, city_name, ires))
                        #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                        elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                            print("3a: Creating {}_{}_{}_tmp.tif file.".format(quad_map_type, city_name, ires))
                        #arcpy function management.Rename with inputs as defined by ESRI
                        arcpy.management.Rename(file_in_raster_quad, file_out_raster_tmp)
            
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #Step: Extract by Mask mosaic raster to remove outlying reference station pixels from quad maps
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #loop through date_hour_str in output_map_date_hour_list 
        for date_hour_str in output_map_date_hour_list:
            #Define path_aoi to folder containing AOI shapefile 
            path_aoi = path_prepData_trunk_fs + '%s/' %(city_name)
            #Define file name for AOI shapefile in Albers_Conical_Equal_Area projection and at resolution 
            file_aoi_shapefile_albers = "%s_alb.shp" %(city_name)
            #Define in mask data 
            in_mask_data = path_aoi + file_aoi_shapefile_albers
            #Define input file of raster for select by mask
            #If file_iTCA_output is TAir_K then add hour to string
            if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                file_in_raster_mosaic_tmp = '/%s_%s_%s_%s_tmp.tif' %(quad_map_type, date_hour_str, city_name, ires)
            #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
            elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                file_in_raster_mosaic_tmp = '/%s_%s_%s_tmp.tif' %(quad_map_type, city_name, ires)
            #Define input path and file for select by mask
            in_raster = path_map_mosaic_output + file_in_raster_mosaic_tmp
            #Define file of raster that will be output for select by mask
            file_out_raster_mosaic = '/%s_%s_%s_%s.tif' %(quad_map_type, date_hour_str, city_name, ires)
            #Define output path and file for select by mask
            out_raster = path_map_mosaic_output + file_out_raster_mosaic

            #Flag_OverwriteOutput used to manage existing output; 0 = will not overwrite existing output; 1 = overwrite existing output
            #If Flag_OverwriteOutput equals 1 then create files regardless if they exist
            if (Flag_OverwriteOutput == 1):
                #If file_iTCA_output is TAir_K then add hour to string
                if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                    print("4: Creating {}_{}_{}_{}.tif file.".format(quad_map_type, date_hour_str, city_name, ires))
                #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                    print("4: Creating {}_{}_{}.tif file.".format(quad_map_type, city_name, ires))
                #arcpy function gp.ExtractByMask_sa with inputs as defined by ESRI
                arcpy.gp.ExtractByMask_sa(in_raster, in_mask_data, out_raster)

            #Else If Flag_OverwriteOutput not equal to 1 then do not create files if they exist
            else:
                #If the file exists, do not create or destroy
                if arcpy.Exists(out_raster):
                    #If file_iTCA_output is TAir_K then add hour to string
                    if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                        print("4: Existing {}_{}_{}_{}.tif".format(quad_map_type, date_hour_str, city_name, ires))
                    #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                    elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                        print("4: Existing {}_{}_{}.tif".format(quad_map_type, city_name, ires))
                #If the file does not exist, create
                else:
                    #If file_iTCA_output is TAir_K then add hour to string
                    if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                        print("4: Creating {}_{}_{}_{}.tif file.".format(quad_map_type, date_hour_str, city_name, ires))
                    #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                    elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                        print("4: Creating {}_{}_{}.tif file.".format(quad_map_type, city_name, ires))
                    #arcpy function gp.ExtractByMask_sa with inputs as defined by ESRI
                    arcpy.gp.ExtractByMask_sa(in_raster, in_mask_data, out_raster)
            
            #Delete temporary mosaic .tif files that contain reference station cells beyond AOI polygon 
            #Define input path and file for select by maskfile_in_raster_mosaic_tmp
            file_check = path_map_mosaic_output + file_in_raster_mosaic_tmp
            #cte 2023 debugging
            print("tmp file_check is {}".format(file_check))
            #If the file exists, do not create or destroy
            if arcpy.Exists(file_check):
                print("5: Deleting {}.".format(file_in_raster_mosaic_tmp))
                #arcpy function Delete_management with inputs as defined by ESRI
                arcpy.Delete_management(path_map_mosaic_output + file_in_raster_mosaic_tmp)
        
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #Step: Delete quad .tif files used to create mosaic 
        #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        #If temporary quad files are to be deleted flag == 1
        if (flag_delete_quad_tif_files == 1):
            #Loop through all iquad_count in quad_ID list
            for iquad_count in quad_ID:
                if (len(quad_ID) > 1):
                    if iquad_count < 10000: iquad_count_str = '%s' %(iquad_count)
                    if iquad_count < 1000: iquad_count_str = '0%s' %(iquad_count)
                    if iquad_count < 100: iquad_count_str = '00%s' %(iquad_count)
                    if iquad_count < 10: iquad_count_str = '000%s' %(iquad_count)
                #define strings: condition when 1 map is processed, when no quads were generated above
                else:
                    iquad_count_str = '%s' %(iquad_count)

                for date_hour_str in output_map_date_hour_list:
                    #If file_iTCA_output is TAir_K then add hour to string
                    if (quad_map_type == 'TAir_K') or (quad_map_type == 'TAir_weighted_K'):
                        #Define file_aoi_resolution_tmp for temporary clipped to AOI DEM raster in Albers projection, which is deleted below
                        file_quad_data_tif = '/%s_%s_%s_%s_%s.tif' %(quad_map_type, date_hour_str, city_name, ires, iquad_count_str)
                    #Else if file_iTCA_output is AspectGround_N_0_rad.asc then do not add hour 
                    elif (quad_map_type == 'AspectGround_N_0_rad.asc') or (quad_map_type == 'SlopeGround_rad.asc'):
                        #Define file_aoi_resolution_tmp for temporary clipped to AOI DEM raster in Albers projection, which is deleted below
                        file_quad_data_tif = '/%s_%s_%s_%s.tif' %(quad_map_type, city_name, ires, iquad_count_str)
                    #cte 2023 debugging
                    print("quad file_check is {}".format(file_check))

                    file_check = path_map_mosaic_output + file_quad_data_tif
                    #If the file exists, do not create or destroy
                    if arcpy.Exists(file_check):
                        print("5: Deleting {}.".format(file_quad_data_tif))
                        #arcpy function Delete_management with inputs as defined by ESRI
                        arcpy.Delete_management(path_map_mosaic_output + file_quad_data_tif)
                #close loop of date_hour_str
            #close loop of iquad_count in quad_ID
    #close loop grid resolutions
#close loop city_list

#arcpy function CheckInExtension checked back in to release feature
arcpy.CheckInExtension("Spatial")