MXD13Q1 annual statistics data, long-term average, max, min and stdev

09_mxd13q1_stats_annual.py

# -*- coding: utf-8 -*-
"""
NAME
    09_mxd13q1_stats_annual.py
    MXD13Q1 annual statistics data, long-term average, max, min and stdev
DESCRIPTION
    Input data for this script will use MXD13Q1 annual data generate from modis_8day2annual.py
    This script can do annual statistics calculation (AVERAGE, MAXIMUM, MINIMUM and STD)
REQUIREMENT
    ArcGIS must installed before using this script, as it required arcpy module.
EXAMPLES
    C:\\Program Files\\ArcGIS\\Pro\\bin\\Python\\envs\\arcgispro-py3\\python 09_mxd13q1_stats_annual.py
NOTES
    This script is designed to work with MODIS naming convention
    If using other data, some adjustment are required: parsing filename and directory
CONTACT
    Benny Istanto
    Climate Geographer
    GOST/DECAT/DECDG, The World Bank
LICENSE
    This script is in the public domain, free from copyrights or restrictions.
VERSION
    $Id$
TODO
    xx
"""
import calendar
import os
import arcpy
from collections import defaultdict
from datetime import datetime, timedelta
import uuid

# To avoid overwriting outputs, set overwriteOutput option to True for this script execution context
arcpy.env.overwriteOutput = True
# Raster environment settings for output consistency and optimization
arcpy.env.compression = "LZ77"
arcpy.env.pyramid = "PYRAMIDS -1 NEAREST LZ77 NO_SKIP"

# ISO3 Country Code
iso3 = "idn"  # Example codes: Syria: syr, Myanmar: mmr, Lebanon: lbn

# Define the range of years
start_year = 2003
end_year = 2022

# Define input and output folders
input_folder = "D:\\temp\\modis\\{}\\gee\\04_fillnullwithstats\\evi_all_annual".format(iso3)
output_folder = "D:\\temp\\modis\\{}\\gee\\03_statistics\\evi_all_annual".format(iso3)

# Global variable to store user's choice
user_choice = None

def set_user_decision():
    """Prompt user for decision on existing files and store it globally."""
    global user_choice
    if user_choice is None:
        decision = input("An output file already exists. Choose an action - Replace (R), Skip (S), Abort (A): ").upper()
        while decision not in ['R', 'S', 'A']:
            print("Invalid choice. Please choose again.")
            decision = input("Choose an action - Replace (R), Skip (S), Abort (A): ").upper()
        user_choice = decision

def process_annualstats(input_folder, output_folder):
    global user_choice

    # Define statistic names
    # Statistics type.
        # MEAN — The mean (average) of the inputs will be calculated.
        # MAJORITY — The majority (value that occurs most often) of the inputs will be determined.
        # MAXIMUM — The maximum (largest value) of the inputs will be determined.
        # MEDIAN — The median of the inputs will be calculated. Note: The input must in integers
        # MINIMUM — The minimum (smallest value) of the inputs will be determined.
        # MINORITY — The minority (value that occurs least often) of the inputs will be determined.
        # RANGE — The range (difference between largest and smallest value) of the inputs will be calculated.
        # STD — The standard deviation of the inputs will be calculated.
        # SUM — The sum (total of all values) of the inputs will be calculated.
        # VARIETY — The variety (number of unique values) of the inputs will be calculated.

    # Define statistic names
    stat_names = {"MAXIMUM": "max", "MINIMUM": "min", "MEAN": "avg", "MEDIAN": "med", "STD": "std"}

    annual_files = []

    # Gather all the annual files for each year
    for year in range(start_year, end_year + 1):
        annual_file_name = f"{iso3}_phy_mxd13q1_evi_mean_{year}.tif"
        annual_file_path = os.path.join(input_folder, annual_file_name)
        
        if arcpy.Exists(annual_file_path):
            annual_files.append(annual_file_path)

    # Calculate the long-term statistics for the annual files
    for stat_type, suffix in stat_names.items():
        out_raster_name = f"{iso3}_phy_mxd13q1_evi_annual_{start_year}_{end_year}_{suffix}.tif"
        out_raster_path = os.path.join(output_folder, out_raster_name)
        
        if arcpy.Exists(out_raster_path):
            if user_choice is None:
                set_user_decision()
            if user_choice == 'S':
                print(f"Skipping existing file: {out_raster_path}")
                continue
            elif user_choice == 'A':
                print("Aborting process.")
                exit()
            elif user_choice == 'R':
                print(f"Replacing existing file: {out_raster_path}")
        else:
            print(f"Processing {out_raster_name}...")

        # Generate a unique temporary filename for intermediate storage
        temp_filename = f"temp_{uuid.uuid4().hex}.tif"
        temp_file_path = os.path.join(arcpy.env.scratchFolder, temp_filename)
        
        # Proceed to save the output
        arcpy.CheckOutExtension("spatial")
        outCellStatistics = arcpy.sa.CellStatistics(annual_files, stat_type, "DATA")
        outCellStatistics.save(temp_file_path)
        
        # Use CopyRaster to save the output raster with specific pixel type and LZ77 compression
        arcpy.management.CopyRaster(
            in_raster=arcpy.sa.Int(temp_file_path), 
            out_rasterdataset=out_raster_path, 
            config_keyword="", 
            background_value="", 
            nodata_value="", 
            onebit_to_eightbit="NONE", 
            colormap_to_RGB="NONE", 
            pixel_type="16_BIT_SIGNED", 
            scale_pixel_value="NONE", 
            RGB_to_Colormap="NONE", 
            format="TIFF", 
            transform="NONE", 
            process_as_multidimensional="CURRENT_SLICE", 
            build_multidimensional_transpose="NO_TRANSPOSE"
        )
        # Attempt to delete the temporary file
        try:
            arcpy.Delete_management(temp_file_path)
            print(f"Temporary file {temp_file_path} deleted.")
        except Exception as e:
            print(f"Warning: Unable to delete temporary file {temp_file_path}")

        # Clear any locks or residual files
        arcpy.ClearWorkspaceCache_management()

        arcpy.CheckInExtension("spatial")
        print(f"{out_raster_name} completed with LZ77 compression.")

# Main function
def main():
    # Call the process_quarters() function for the input folder
    process_annualstats(input_folder, output_folder)

if __name__ == '__main__':
    main()

print("Script completed.")