TIMESAT output processing from different scenario, 2 year and 3 year

modis_pheno_remap_2y.py

import arcpy
import os

# Input and output directories
input_folder = "X:\\Temp\\modis\\lbn\\gee\\11_timesat_raw\\geotiff\\step01_raw\\2y"
output_folder = "X:\\Temp\\modis\\lbn\\gee\\11_timesat_raw\\geotiff\\step02_remap\\2y"
# Temporary folder to store intermediate rasters
temp_folder = "X:\\Temp\\modis\\lbn\\gee\\11_timesat_raw\\geotiff\\temp"  

# To avoid overwriting outputs, change overwriteOutput option to False.
arcpy.env.overwriteOutput = True

# Check out the Spatial Analyst extension
arcpy.CheckOutExtension("spatial")

# Create output folder if it doesn't exist
if not os.path.exists(output_folder):
    os.makedirs(output_folder)

# Define reclassification range for 'start' files
reclass_start = arcpy.sa.RemapValue([
    [47, 1], [48, 2], [49, 3], [50, 4], [51, 5], [52, 6],
    [53, 7], [54, 8], [55, 9], [56, 10], [57, 11], [58, 12],
    [59, 13], [60, 14], [61, 15], [62, 16], [63, 17], [64, 18],
    [65, 19], [66, 20], [67, 21], [68, 22], [69, 23]
])

# Process each 'start' file
for filename in os.listdir(input_folder):
    if filename.endswith('.tif'):
        print(f"\nProcessing: {filename}")
        input_file_path = os.path.join(input_folder, filename)

        # Build raster attribute table if not present or if it needs to be refreshed
        if not arcpy.Exists(input_file_path + ".vat.dbf"):
            print(f"Building attribute table for {filename}...")
        else:
            # Delete the existing raster attribute table
            print(f"Deleting existing attribute table for {filename}...")
            arcpy.Delete_management(input_file_path + ".vat.dbf")
            
        # Build a new raster attribute table
        arcpy.BuildRasterAttributeTable_management(input_file_path)

        # Extract values based on filename type
        if "start" in filename:
            print(f"Extracting attributes for 'start' file {filename}.")
            extracted = arcpy.sa.ExtractByAttributes(input_file_path, "VALUE >= 47 AND VALUE <= 69")
            extracted_start_path = os.path.join(temp_folder, "extracted_start_" + filename)  # Save to temp folder
            extracted.save(extracted_start_path)
        elif "end" in filename:
            print(f"Extracting attributes for 'end' file {filename}.")
            extracted = arcpy.sa.ExtractByAttributes(input_file_path, "VALUE >= 24 AND VALUE <= 46")
            extracted_end_path = os.path.join(temp_folder, "extracted_end_" + filename)  # Save to temp folder
            extracted.save(extracted_end_path)
        else:
            print(f"Skipping file {filename}.")
            continue  # Skip files that don't match 'start' or 'end'

        # Reclassify 'start' file if it's a 'start' file
        if "start" in filename:
            reclassified_start = arcpy.sa.Reclassify(extracted, "Value", reclass_start, "NODATA")
            reclassified_start.save(os.path.join(temp_folder, "reclassified_" + filename))  # Save to temp folder

            # Find corresponding 'end' file
            end_filename = filename.replace("start", "end")
            extracted_end_path = os.path.join(temp_folder, "extracted_end_" + end_filename)

            # Check if corresponding 'end' file exists
            if os.path.exists(extracted_end_path):
                # Mosaic 'start' and 'end' rasters
                output_filename = filename.replace("start_", "").replace("end_", "")
                output_file_path = os.path.join(output_folder, output_filename)

                # The ~ operator inverts the boolean raster output from IsNull, effectively giving us the result for "IsNotNull".
                combined_raster = arcpy.sa.Con(~arcpy.sa.IsNull(reclassified_start), reclassified_start, extracted_end_path)
                combined_raster.save(output_file_path)
                print(f"Combined '{reclassified_start}' and '{os.path.basename(extracted_end_path)}' to '{output_filename}'.")
            else:
                print(f"Corresponding 'end' file '{end_filename}' does not exist. Skipping combine process for '{filename}'.")

# Check in the Spatial Analyst extension
arcpy.CheckInExtension("spatial")

print("Processing complete!")

modis_pheno_remap_3y.py

import arcpy
import os

# Input and output directories
input_folder = "X:\\Temp\\modis\\lbn\\gee\\11_timesat_raw\\geotiff\\step01_raw\\3y"
output_folder = "X:\\Temp\\modis\\lbn\\gee\\11_timesat_raw\\geotiff\\step02_remap\\3y"

# To avoid overwriting outputs, change overwriteOutput option to False.
arcpy.env.overwriteOutput = True

# Check out the Spatial Analyst extension
arcpy.CheckOutExtension("spatial")

# Create output folder if it doesn't exist
if not os.path.exists(output_folder):
    os.makedirs(output_folder)

# Define reclassification range
reclass_values = arcpy.sa.RemapValue([
    [47, 1], [48, 2], [49, 3], [50, 4], [51, 5], [52, 6], [53, 7],
    [54, 8], [55, 9], [56, 10], [57, 11], [58, 12], [59, 13], [60, 14],
    [61, 15], [62, 16], [63, 17], [64, 18], [65, 19], [66, 20], [67, 21],
    [68, 22], [69, 23], [70, 24], [71, 25], [72, 26], [73, 27], [74, 28],
    [75, 29], [76, 30], [77, 31], [78, 32], [79, 33], [80, 34], [81, 35],
    [82, 36], [83, 37], [84, 38], [85, 39], [86, 40], [87, 41], [88, 42],
    [89, 43], [90, 44], [91, 45], [92, 46]
])

# Process each file
for filename in os.listdir(input_folder):
    if filename.endswith('.tif'):
        print(f"\nProcessing: {filename}")
        input_file_path = os.path.join(input_folder, filename)

        # Build raster attribute table if not present or if it needs to be refreshed
        if arcpy.Exists(input_file_path + ".vat.dbf"):
            # Delete the existing raster attribute table
            print(f"Deleting existing attribute table for {filename}...")
            arcpy.Delete_management(input_file_path + ".vat.dbf")

        # Build a new raster attribute table
        print(f"Building attribute table for {filename}...")
        arcpy.BuildRasterAttributeTable_management(input_file_path)

        # Extract values
        print(f"Extracting attributes for {filename}.")
        extracted = arcpy.sa.ExtractByAttributes(input_file_path, "VALUE >= 47 AND VALUE <= 92")
        
        # Reclassify
        reclassified = arcpy.sa.Reclassify(extracted, "Value", reclass_values, "NODATA")
        
        # Save the reclassified raster
        reclassified.save(os.path.join(output_folder, filename))
        print(f"Saved reclassified raster as '{filename}' in the output folder.")

# Check in the Spatial Analyst extension
arcpy.CheckInExtension("spatial")

print("Processing complete!")

TIMESAT_output_remap.md

This script required GeoTIFF file from TIMESAT seas2img binary output conversion process for 2-year scenario. It will remap the raster value and combined between two simulation exercise.