MXD13Q1 monthly mean from 8-day data

04_mxd13q1_8day2monthly.py

# -*- coding: utf-8 -*-
"""
NAME
    08_mxd13q1_stats_quarter.py
    MXD13Q1 quarter statistics data, long-term average, max, min and stdev
DESCRIPTION
    Input data for this script will use MXD13Q1 quarter data generate from modis_quarter.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 08_mxd13q1_stats_quarter.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, change overwriteOutput option to False.
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" # Syria: syr, Myanmar: mmr, Lebanon: lbn

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

# Change the data and output folder
input_folder = "D:\\temp\\modis\\{}\\gee\\04_fillnullwithstats\\evi_all_quarterly".format(iso3)
output_folder = "D:\\temp\\modis\\{}\\gee\\03_statistics\\evi_all_quarterly".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_quarterstats(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 the statistic names for file naming and their corresponding arcpy names for processing
    stat_names = {"MAXIMUM": "max", "MINIMUM": "min", "MEAN": "avg", "MEDIAN": "med", "STD": "std"}

    # Quarters definition
    quarters_def = ["010203", "020304", "030405", "040506", "050607", "060708", "070809", "080910", "091011", "101112", "111201", "120102"]

    for quarter in quarters_def:
        quarterly_files = []

        # Gather all the quarterly mean files for each quarter across all years
        for year in range(start_year, end_year + 1):
            quarterly_mean_name = f"{iso3}_phy_mxd13q1_evi_mean_{year}_{quarter}.tif"
            quarterly_mean_path = os.path.join(input_folder, quarterly_mean_name)
            
            if arcpy.Exists(quarterly_mean_path):
                quarterly_files.append(quarterly_mean_path)

        # Calculate the long-term statistics for the quarter
        for stat_type, suffix in stat_names.items():
            out_raster_name = f"{iso3}_phy_mxd13q1_evi_{quarter}_{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(quarterly_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_quarterstats(input_folder, output_folder)

if __name__ == '__main__':
    main()

print("Script completed.")