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

#the data causes this scripot to have memory issues on laptops
#use this to get around it if you have enough hard drive space to
#substitute as ram
#Source:https://stackoverflow.com/questions/57507832/unable-to-allocate-array-with-shape-and-data-type
#Windows 10
#1. Press the Windows key
#2. Type SystemPropertiesAdvanced
#3. Click run as administrator
#4. Under performance click settings
#5. Select advanced tab
#6. Select change
#7. Uncheck automatically managing paging file size for all drives
#8. Then select Custom size and fill in the appropriate value
#9. Then press OK and exit everything and reboot
#10. Then run script

#rasterio relevant links
#https://gist.github.com/cordmaur/b588da93e2ca7f5c23ed248a4fb03cdd#file-pfg_05_04-py
#https://github.com/mapbox/rasterio/issues/2048

# I ran gdal translate on the nlcd rasters to get them into geotiff
# like this
#gdal_translate -of GTiff  C:\temp\input\a.img   C:\temp\output\a.tif

#gdal_translate -of GTiff  Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2001_land_cover_l48_20210604.img   Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2001_land_cover_l48_20210604.tif
#gdal_translate -of GTiff  Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2004_land_cover_l48_20210604.img   Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2004_land_cover_l48_20210604.tif
#gdal_translate -of GTiff  Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2006_land_cover_l48_20210604.img   Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2006_land_cover_l48_20210604.tif
#gdal_translate -of GTiff  Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2008_land_cover_l48_20210604.img   Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2008_land_cover_l48_20210604.tif
#gdal_translate -of GTiff  Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2011_land_cover_l48_20210604.img   Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2011_land_cover_l48_20210604.tif
#gdal_translate -of GTiff  Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2013_land_cover_l48_20210604.img   Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2013_land_cover_l48_20210604.tif
#gdal_translate -of GTiff  Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2016_land_cover_l48_20210604.img   Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2016_land_cover_l48_20210604.tif
#gdal_translate -of GTiff  Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2019_land_cover_l48_20210604.img   Y:\UHI\data\phase2_temp_storage\RASTERS_NLCD\nlcd_2019_land_cover_l48_20210604.tif

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

import json
import math
import os
import shutil
import time

import rasterio
from rasterio.mask import mask
from rasterio.io import MemoryFile
from rasterio.merge import merge
import gdal_merge as gm

from osgeo import osr
import pycrs
import pyproj
import numpy as np
import fiona

from joblib import Parallel, delayed

import datetime as dt
import argparse


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

lcmap_pri_base_path = "/caldera/projects/usgs/eros/urban_heat_islands/Phase2/section_332/RASTERS_LCMAP_PRI"
lcmap_conf_base_path = "/caldera/projects/usgs/eros/urban_heat_islands/Phase2/section_332/RASTERS_LCMAP_CONF"
nlcd_base_path = "/caldera/projects/usgs/eros/urban_heat_islands/Phase2/section_332/RASTERS_NLCD"
base_output_directory = "/caldera/projects/usgs/eros/urban_heat_islands/Phase2/section_332/output"
temp_directory = "/caldera/projects/usgs/eros/urban_heat_islands/Phase2/section_332/temp"
full_path_to_shapefile = "/caldera/projects/usgs/eros/urban_heat_islands/Phase2/section_332/ard_tile_shapefile/ard_grd_conus_final_pjr.shp"

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

def write_urban_nonurban_raster(output_directory_path,output_filename,temp_directory,raster_path_nlcd,raster_path_lcmap_pri,raster_path_lcmap_conf,hv_filtered_features_dict):

    dataset_lcmap_pri = rasterio.open(raster_path_lcmap_pri)
    dataset_lcmap_conf = rasterio.open(raster_path_lcmap_conf)
    dataset_nlcd = rasterio.open(raster_path_nlcd)

    tile_paths_to_merge = []

    for (k,v) in hv_filtered_features_dict.items():

        tile_name = k
        hv_geometries = []
        feature = v
        #only do one at a time
        hv_geometries.append(feature["geometry"])
        print("attempting "+tile_name)#remove later

        #out_lcmap_pri, out_transform_lcmap_pri = rasterio.mask.mask(dataset_lcmap_pri, hv_geometries, crop=True)

        '''
        out_lcmap_pri, out_transform_lcmap_pri = mask(dataset_lcmap_pri, hv_geometries, crop=True,pad=False, pad_width=0.0,all_touched=True)
        out_meta_lcmap_pri = dataset_lcmap_pri.meta
        out_lcmap_conf, out_transform_lcmap_conf = mask(dataset_lcmap_conf, hv_geometries, crop=True,pad=False, pad_width=0.0,all_touched=True)
        out_meta_lcmap_conf = dataset_lcmap_conf.meta
        out_nlcd, out_transform_nlcd = mask(dataset_nlcd, hv_geometries, crop=True,pad=False, pad_width=0.0,all_touched=True)
        out_meta_nlcd = dataset_nlcd.meta
        '''

        out_lcmap_pri, out_transform_lcmap_pri = mask(dataset_lcmap_pri, hv_geometries, crop=True,pad=False, pad_width=0.0,all_touched=False)
        out_meta_lcmap_pri = dataset_lcmap_pri.meta
        out_lcmap_conf, out_transform_lcmap_conf = mask(dataset_lcmap_conf, hv_geometries, crop=True,pad=False, pad_width=0.0,all_touched=False)
        out_meta_lcmap_conf = dataset_lcmap_conf.meta
        # out_nlcd, out_transform_nlcd = mask(dataset_nlcd, hv_geometries, crop=True,pad=False, pad_width=0.0,all_touched=False)
        # out_meta_nlcd = dataset_nlcd.meta
        out_nlcd, out_transform_nlcd = mask(dataset_nlcd, hv_geometries, crop=True,pad=True, pad_width=1.0,all_touched=True)
        out_meta_nlcd = dataset_nlcd.meta
        
        '''
        out_lcmap_pri, out_transform_lcmap_pri = mask(dataset_lcmap_pri, hv_geometries, crop=True,pad=True, pad_width=0.1,all_touched=False)
        out_meta_lcmap_pri = dataset_lcmap_pri.meta
        out_lcmap_conf, out_transform_lcmap_conf = mask(dataset_lcmap_conf, hv_geometries, crop=True,pad=True, pad_width=0.1,all_touched=False)
        out_meta_lcmap_conf = dataset_lcmap_conf.meta
        out_nlcd, out_transform_nlcd = mask(dataset_nlcd, hv_geometries, crop=True,pad=True, pad_width=0.1,all_touched=False)
        out_meta_nlcd = dataset_nlcd.meta
        '''


        designated_output_nodata_val = -9999

        new_image = np.full(out_lcmap_pri.shape,designated_output_nodata_val)
        new_out_meta = out_meta_lcmap_pri
        new_out_meta['dtype'] = 'int32'

        for n1 in range(new_image.shape[1]):
            for n2 in range(new_image.shape[2]):

                lcmap_pri_val = out_lcmap_pri[0][n1][n2]
                lcmap_conf_val = out_lcmap_conf[0][n1][n2]

                # nlcd has nodata spatially distributed in a way that
                # results in it having smaller grid
                # for certain areas like the edges in miami
                # ["H27V18","H27V19"]
                # we dont have all sets of pixels to compare here
                # so if it fails this check the output pixel is set to nodata
                out_of_bounds_for_nlcd = False
                if (n1 > out_nlcd.shape[1]-1) or (n2 > out_nlcd.shape[2]-1):
                    out_of_bounds_for_nlcd = True
                    #print("out_of_bounds_for_nlcd")

                lcmap_pri_nodata_val = 0
                nlcd_nodata_val = 0

                o_val = None
                if not out_of_bounds_for_nlcd:

                    #+1 bc we pad nlcd, and padding by 1 adds 1 on both sides
                    #so it adds two total
                    #but we want to align at the same place as lcmap
                    #so we pad
                    #nlcd has smaller pixels, so it comes up short on the geometries edges sometimes
                    #so we padd it a little, then use offsets here so stuff still aligns
                    nlcd_val = out_nlcd[0][n1+1][n2+1]


                    #==============================================================================================

                    """
                    Chase's Table
                    nlcd_value | lcmap_value | output
                    developed open space 21 | developed 1 | 21
                    developed low intensity 22 | developed 1 | 22
                    developed medium intensity 23 | developed 1| 23
                    developed high intensity 24 | developed 1 | 24
                    not in [21,23,23,24] | developed 1 | 25
                    hay/pasture cultivated crops 81-82 | cropland 2 | 2
                    shrub/scrub, grassland/herbaceous 52,71 | grass/shrub 3 | 3
                    deciduous forest, evergreen forest, mixed forest 41,42,43 | treecover 4 | 4
                    water 11 | water 5 | 5
                    woody wetland, emergent herbaceous wetland 90,95 | wetland 6 | 6
                    ice/snow 12 | ice and snow 7 | 7
                    barren 31 | barren 8| 8
                    
                    addenum 11/8
                    [9:32 AM] Mueller, Chase (Contractor)
                    Hey Reza, I'm still finishing up a tool to do a full stack verification (so we can easily do this for C2 forward), 
                    but I did find the following
                    -Non Urban NLCD with Urban LCMAP with Conf < 99 was labelled as 25 when it should have been 2. 
                    It should only be 25 when Conf is 99 or 100 with Urban LCMAP and non urban NLCD
                    
                    addenum 11/15
                    When NLCD is non urban and LCMAP is urban, only the 99-100 is addressed. 
                    If less than 99, it should return the NLCD non urban class. 
                    So I believe that needs to be added.
                    
                    The only scenario where non-urban NLCD gets cross-mapped to LCMAP values is when conf is < 99 and lcmap is Urban. 
                    So you could just update those last set of conditions so that nlcd_pri_val=[what you currently have] and 
                    lcmap_pri_val = 1 and conf < 99
                    
                    and if both are non urban, then they simply use the LCMAP version.  
                    There's no condition for that. Conf is irrelevant here.
          
                    if NLCD is urban and LCMAP is non urban with conf < 99, set it to LCMAP val. 
          
                    """

                    if nlcd_val == 21 and lcmap_pri_val == 1:
                        o_val = 21
                    if nlcd_val == 22 and lcmap_pri_val == 1:
                        o_val = 22
                    if nlcd_val == 23 and lcmap_pri_val == 1:
                        o_val = 23
                    if nlcd_val == 24 and lcmap_pri_val == 1:
                        o_val = 24

                    if nlcd_val not in [21, 22, 23, 24] and lcmap_pri_val == 1 and lcmap_conf_val > 98:
                        o_val = 25

                    #new
                    if nlcd_val not in [21, 22, 23, 24] and lcmap_pri_val == 1 and lcmap_conf_val < 99:
                        o_val = nlcd_val


                    if nlcd_val in [81, 82] and lcmap_pri_val == 1 and lcmap_conf_val < 99:
                        o_val = 2
                    if nlcd_val in [52,71] and lcmap_pri_val == 1 and lcmap_conf_val < 99:
                        o_val = 3
                    if nlcd_val in [41, 42, 43] and lcmap_pri_val == 1 and lcmap_conf_val < 99:
                        o_val = 4
                    if nlcd_val == 11 and lcmap_pri_val == 1 and lcmap_conf_val < 99:
                        o_val = 5
                    if nlcd_val in [90, 95] and lcmap_pri_val == 1 and lcmap_conf_val < 99:
                        o_val = 6
                    if nlcd_val == 12 and lcmap_pri_val == 1 and lcmap_conf_val < 99:
                        o_val = 7
                    if nlcd_val == 31 and lcmap_pri_val == 1 and lcmap_conf_val < 99:
                        o_val = 8

                    if nlcd_val not in [21, 22, 23, 24] and lcmap_pri_val != 1:
                        o_val = lcmap_pri_val

                    if nlcd_val in [21, 22, 23, 24] and lcmap_pri_val != 1 and lcmap_conf_val < 99:
                        o_val = lcmap_pri_val

                    if o_val is None:
                        o_val = lcmap_pri_val

                    if o_val == 1:
                        print("1o_val is 1 " + str(lcmap_pri_val) + " " + str(lcmap_conf_val) + " " + str(nlcd_val))
                        ABC = None

                    #===============================================================================================

                else:
                    o_val = designated_output_nodata_val

                if o_val == 1:
                    print("2o_val is 1 " + str(lcmap_pri_val) + " " + str(lcmap_conf_val) + " " + str(nlcd_val))
                    ABC = None
                
                if o_val == lcmap_pri_nodata_val or o_val == nlcd_nodata_val:
                    o_val = designated_output_nodata_val

                new_image[0][n1][n2] = o_val
        
        if 0 in new_image:
            print("0_is_in_new_image_premerge")

        new_out_meta.update({"driver": "GTiff",
                         "height": out_lcmap_pri.shape[1],
                         "width": out_lcmap_pri.shape[2],
                         "transform": out_transform_lcmap_pri,
                         "nodata":designated_output_nodata_val,
                         "compress": "LZW"})

        tile_output_path = os.path.join(temp_directory,tile_name+"_"+output_filename)

        with rasterio.open(tile_output_path, "w", **new_out_meta) as dest:
            dest.write(new_image.astype(np.int32))
            print(tile_output_path)
            tile_paths_to_merge.append(tile_output_path)
            dest.close()

    merged_path = os.path.join(output_directory_path,output_filename)
    
    '''
    datasets_to_merge = []
    for tile_path in tile_paths_to_merge:
        tile_dataset = rasterio.open(tile_path)
        datasets_to_merge.append(tile_dataset)

    
    print(merged_path)
    merge(datasets_to_merge,dst_path=merged_path)

    for a_dataset in datasets_to_merge:
        a_dataset.close()
    '''

    merge_command = ['', "-o", merged_path,'-a_nodata',"-9999"]
    merge_command.extend(tile_paths_to_merge)
    gm.main(merge_command)


    dataset_lcmap_pri.close()
    dataset_lcmap_conf.close()
    dataset_nlcd.close()



#----------------------------------------------------------------------------------------------------------|
sites_tuples = [
    ("/caldera/projects/usgs/eros/urban_heat_islands/Austin", ["H15V16", "H16V16", "H15V17"], "AUS", "Austin"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/SiouxFalls", ["H16V06"], "FSD", "SiouxFalls"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Missoula", ["H07V03", "H08V03"], "MSO", "Missoula"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Albuquerque", ["H10V12", "H10V13"], "ABQ", "Albuquerque"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Atlanta", ["H23V13", "H23V14", "H24V13", "H24V14"], "ATL","Atlanta"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Birmingham", ["H22V13", "H22V14"], "BHX", "Birmingham"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Nashville", ["H22V11", "H22V12"], "BNA", "Nashville"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Boise", ["H06V05"], "BOI", "Boise"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Boston", ["H30V05", "H30V06"], "BOS", "Boston"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Baltimore", ["H27V08", "H27V09", "H28V08", "H28V09"], "BWI","Baltimore"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Charlotte", ["H25V12", "H26V12", "H26V11"], "CLT", "Charlotte"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/ColoradoSprings", ["H11V10", "H12V10", "H11V09", "H12V09"], "COS","ColoradoSprings"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Columbus", ["H24V08", "H24V09"], "CMH", "Columbus"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Cheyenne", ["H12V08"], "CYS", "Cheyenne"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Washington", ["H27V09", "H28V09", "H27V08"], "DCA", "Washington"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Denver", ["H11V09", "H12V09"], "DEN", "Denver"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Dallas", ["H16V14", "H16V15"], "DFW", "Dallas"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/DesMoines", ["H18V08"], "DSM", "DesMoines"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Detroit", ["H23V06", "H24V06", "H23V07", "H24V07"], "DTW","Detroit"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/ElPaso", ["H10V15"], "ELP", "ElPaso"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Fargo", ["H16V04"], "FAR", "Fargo"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/SiouxFalls", ["H16V06"], "FSD", "SiouxFalls"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Houston", ["H17V16", "H17V17"], "IAH", "Houston"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Indianapolis", ["H22V09"], "IND", "Indianapolis"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Jacksonville", ["H26V15", "H26V16"], "JAX", "Jacksonville"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/LasVegas", ["H05V11"], "LAS", "LasVegas"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/LosAngeles", ["H03V12", "H04V12"], "LAX", "LosAngeles"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/NewYork", ["H28V07", "H29V06", "H29V07", "H29V08"], "LGA","NewYork"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/LittleRock", ["H19V13"], "LIT", "LittleRock"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Cincinnati", ["H23V09", "H23V10"], "LUK", "Cincinnati"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/KansasCity", ["H17V09", "H17V10", "H18V10"], "MCI", "KansasCity"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Memphis", ["H20V12", "H20V13"], "MEM", "Memphis"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Miami",["H27V18","H27V19"],"MIA","Miami"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Milwaukee", ["H21V06", "H21V07"], "MKE", "Milwaukee"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Missoula", ["H07V03", "H08V03"], "MSO", "Missoula"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Minneapolis", ["H18V05", "H18V06"], "MSP", "Minneapolis"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/NewOrleans", ["H20V16", "H20V17", "H21V16", "H21V17"], "MSY","NewOrleans"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/OklahomaCity", ["H16V12", "H16V13"], "OKC", "OklahomaCity"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Omaha", ["H16V08", "H17V08"], "OMA", "Omaha"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Chicago", ["H21V07", "H21V08", "H22V08"], "ORD", "Chicago"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Portland", ["H03V03"], "PDX", "Portland"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Philadelphia", ["H28V07", "H28V08", "H29V08"], "PHL","Philadelphia"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Phoenix", ["H06V13", "H07V13"], "PHX", "Phoenix"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Pittsburgh", ["H25V08", "H26V08"], "PIT", "Pittsburgh"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Raleigh", ["H27V11"], "RDU", "Raleigh"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/SanDiego", ["H03V13", "H04V12", "H04V13"], "SAN", "SanDiego"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/SanAntonio", ["H15V17"], "SAT", "SanAntonio"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Louisville", ["H22V10", "H23V10"], "SDF", "Louisville"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Seattle", ["H03V02", "H03V01", "H04V01", "H04V02"], "SEA","Seattle"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/SanFrancisco", ["H01V08", "H01V09", "H02V08", "H02V09"], "SFO","SanFrancisco"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/SaltLakeCity", ["H08V07", "H08V08"], "SLC", "SaltLakeCity"),    
    ("/caldera/projects/usgs/eros/urban_heat_islands/Sacramento", ["H02V08", "H03V08"], "SMF", "Sacramento"),
    ("/caldera/projects/usgs/eros/urban_heat_islands/Cincinnati", ["H23V09", "H23V10"], "CVG", "Cincinnati")

]

"""

    
    


"""


raster_groupings = [
    (1985,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1985_V11_LCPRI.tif", "LCMAP_CU_1985_V11_LCPCONF.tif")

    ,
(1986,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1986_V11_LCPRI.tif", "LCMAP_CU_1986_V11_LCPCONF.tif"),
(1987,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1987_V11_LCPRI.tif", "LCMAP_CU_1987_V11_LCPCONF.tif"),
(1988,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1988_V11_LCPRI.tif", "LCMAP_CU_1988_V11_LCPCONF.tif"),
(1989,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1989_V11_LCPRI.tif", "LCMAP_CU_1989_V11_LCPCONF.tif"),
(1990,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1990_V11_LCPRI.tif", "LCMAP_CU_1990_V11_LCPCONF.tif"),
(1991,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1991_V11_LCPRI.tif", "LCMAP_CU_1991_V11_LCPCONF.tif"),
(1992,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1992_V11_LCPRI.tif", "LCMAP_CU_1992_V11_LCPCONF.tif"),
(1993,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1993_V11_LCPRI.tif", "LCMAP_CU_1993_V11_LCPCONF.tif"),
(1994,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1994_V11_LCPRI.tif", "LCMAP_CU_1994_V11_LCPCONF.tif"),
(1995,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1995_V11_LCPRI.tif", "LCMAP_CU_1995_V11_LCPCONF.tif"),
(1996,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1996_V11_LCPRI.tif", "LCMAP_CU_1996_V11_LCPCONF.tif"),
(1997,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1997_V11_LCPRI.tif", "LCMAP_CU_1997_V11_LCPCONF.tif"),
(1998,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1998_V11_LCPRI.tif", "LCMAP_CU_1998_V11_LCPCONF.tif"),
(1999,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_1999_V11_LCPRI.tif", "LCMAP_CU_1999_V11_LCPCONF.tif"),
(2000,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_2000_V11_LCPRI.tif", "LCMAP_CU_2000_V11_LCPCONF.tif"),
(2001,"nlcd_2001_land_cover_l48_20210604.tif", "LCMAP_CU_2001_V11_LCPRI.tif", "LCMAP_CU_2001_V11_LCPCONF.tif"),


(2002,"nlcd_2004_land_cover_l48_20210604.tif", "LCMAP_CU_2002_V11_LCPRI.tif", "LCMAP_CU_2002_V11_LCPCONF.tif"),

(2003,"nlcd_2004_land_cover_l48_20210604.tif", "LCMAP_CU_2003_V11_LCPRI.tif", "LCMAP_CU_2003_V11_LCPCONF.tif"),

(2004,"nlcd_2004_land_cover_l48_20210604.tif", "LCMAP_CU_2004_V11_LCPRI.tif", "LCMAP_CU_2004_V11_LCPCONF.tif"),

#using 2006
(2005,"nlcd_2006_land_cover_l48_20210604.tif", "LCMAP_CU_2005_V11_LCPRI.tif", "LCMAP_CU_2005_V11_LCPCONF.tif"),

(2006,"nlcd_2006_land_cover_l48_20210604.tif", "LCMAP_CU_2006_V11_LCPRI.tif", "LCMAP_CU_2006_V11_LCPCONF.tif"),

#using 2008
(2007,"nlcd_2008_land_cover_l48_20210604.tif", "LCMAP_CU_2007_V11_LCPRI.tif", "LCMAP_CU_2007_V11_LCPCONF.tif"),

(2008,"nlcd_2008_land_cover_l48_20210604.tif", "LCMAP_CU_2008_V11_LCPRI.tif", "LCMAP_CU_2008_V11_LCPCONF.tif"),
(2009,"nlcd_2011_land_cover_l48_20210604.tif", "LCMAP_CU_2009_V11_LCPRI.tif", "LCMAP_CU_2009_V11_LCPCONF.tif"),
(2010,"nlcd_2011_land_cover_l48_20210604.tif", "LCMAP_CU_2010_V11_LCPRI.tif", "LCMAP_CU_2010_V11_LCPCONF.tif"),
(2011,"nlcd_2011_land_cover_l48_20210604.tif", "LCMAP_CU_2011_V11_LCPRI.tif", "LCMAP_CU_2011_V11_LCPCONF.tif"),
(2012,"nlcd_2013_land_cover_l48_20210604.tif", "LCMAP_CU_2012_V11_LCPRI.tif", "LCMAP_CU_2012_V11_LCPCONF.tif"),
(2013,"nlcd_2013_land_cover_l48_20210604.tif", "LCMAP_CU_2013_V11_LCPRI.tif", "LCMAP_CU_2013_V11_LCPCONF.tif"),
(2014,"nlcd_2016_land_cover_l48_20210604.tif", "LCMAP_CU_2014_V11_LCPRI.tif", "LCMAP_CU_2014_V11_LCPCONF.tif"),
(2015,"nlcd_2016_land_cover_l48_20210604.tif", "LCMAP_CU_2015_V11_LCPRI.tif", "LCMAP_CU_2015_V11_LCPCONF.tif"),
(2016,"nlcd_2016_land_cover_l48_20210604.tif", "LCMAP_CU_2016_V11_LCPRI.tif", "LCMAP_CU_2016_V11_LCPCONF.tif"),
(2017,"nlcd_2019_land_cover_l48_20210604.tif", "LCMAP_CU_2017_V11_LCPRI.tif", "LCMAP_CU_2017_V11_LCPCONF.tif"),
(2018,"nlcd_2019_land_cover_l48_20210604.tif", "LCMAP_CU_2018_V11_LCPRI.tif", "LCMAP_CU_2018_V11_LCPCONF.tif"),
(2019,"nlcd_2019_land_cover_l48_20210604.tif", "LCMAP_CU_2019_V11_LCPRI.tif", "LCMAP_CU_2019_V11_LCPCONF.tif"),
(2020, "nlcd_2019_land_cover_l48_20210604.tif", "LCMAP_CU_2020_V12_LCPRI.tif", "LCMAP_CU_2020_V12_LCPCONF.tif")
]

"""
raster_groupings = [
    (2020, "nlcd_2019_land_cover_l48_20210604.tif", "LCMAP_CU_2020_V12_LCPRI.tif", "LCMAP_CU_2020_V12_LCPCONF.tif")
]
"""

#where there is not an nlcd matching year for the lcmap
#use the nearest future nlcd date


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

parser = argparse.ArgumentParser()
parser.add_argument("--site_tuples_index", "-s", help="site_tuples_index")
args = parser.parse_args()
#do_step_344(args.raster_step_332_1985_path,args.raster_step_332_2020_path,args.final_output_directory,args.full_path_to_dissolved_city_boundaries_shapefile)


sites_tuples = [sites_tuples[int(args.site_tuples_index)]]

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


new_site_tuples = []

for (path,hv_list,airport_code,city) in sites_tuples:
    hv_required_list = []
    for hv in hv_list:
        h = hv[1:3]
        v = hv[4:6]
        new_format_hv = "h"+str(int(h))+"v"+str(int(v))
        hv_required_list.append(new_format_hv)

    hv_filtered_features_dict = {}
    with fiona.open(full_path_to_shapefile, "r") as shapefile:
        for feature in shapefile:
            h_v = feature["properties"]["h_v"]
            if h_v in hv_required_list:
                hv_filtered_features_dict[h_v] = feature

    new_tuple = (hv_list,airport_code,city,hv_filtered_features_dict)
    new_site_tuples.append(new_tuple)


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

start_time = time.time()



#if os.path.exists(base_output_directory):
#    shutil.rmtree(base_output_directory)
#    time.sleep(10)

if not os.path.exists(base_output_directory):
    os.mkdir(base_output_directory)

if not os.path.exists(temp_directory):
    #shutil.rmtree(temp_directory)
    #time.sleep(10)
    os.mkdir(temp_directory)

process_tuples = []

def make_rasters_for_site(process_tuple):
    (output_directory,output_filename,temp_directory,raster_path_nlcd,raster_path_lcmap_pri,raster_path_lcmap_conf,hv_filtered_features_dict) = process_tuple
    write_urban_nonurban_raster(output_directory,output_filename,temp_directory,raster_path_nlcd,raster_path_lcmap_pri,raster_path_lcmap_conf,hv_filtered_features_dict)


for (year,raster_file_name_ncld,raster_file_name_lcmap_pri,raster_file_name_lcmap_conf) in raster_groupings:

    raster_path_lcmap_pri = os.path.join(lcmap_pri_base_path, raster_file_name_lcmap_pri)
    raster_path_lcmap_conf = os.path.join(lcmap_conf_base_path, raster_file_name_lcmap_conf)
    raster_path_nlcd = os.path.join(nlcd_base_path, raster_file_name_ncld)

    for (hv_list,airport_code,city,hv_filtered_features_dict) in new_site_tuples:

        #UHI_CU_FSD_1984_20211026_C01_LC.tif
        ard_region = "CU"
        processing_date_label = dt.datetime.today().strftime('%Y%m%d')
        ard_collection_version = "C01"

        output_filename = "UHI_CU_"+str(airport_code)+"_"+str(year)+"_"+str(processing_date_label)+"_C01_LC.tif"
        output_directory = os.path.join(base_output_directory, airport_code)

        temp_city_directory = os.path.join(temp_directory, airport_code)

        if os.path.exists(output_directory):
            #shutil.rmtree(output_directory)
            time.sleep(1)

        if os.path.exists(temp_city_directory):
            #shutil.rmtree(temp_city_directory)
            time.sleep(1)

        if not os.path.exists(output_directory):
            os.mkdir(output_directory)

        if not os.path.exists(temp_city_directory):
            time.sleep(1)
            os.mkdir(temp_city_directory)

        a_process_tuple = (output_directory,output_filename,temp_city_directory,raster_path_nlcd,raster_path_lcmap_pri,raster_path_lcmap_conf,hv_filtered_features_dict)
        process_tuples.append(a_process_tuple)

#for a_process_tuple in process_tuples:
#    make_rasters_for_site(a_process_tuple)

#if you do 40 it gets out of memory
Parallel(n_jobs=20)(delayed(make_rasters_for_site)(a_site_tuple) for a_site_tuple in process_tuples)






end_time = time.time()
time_taken = (end_time-start_time)/60

print(time_taken)