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

02_mxd13q1_stats_8day.py

# -*- coding: utf-8 -*-
"""
NAME
    02_mxd13q1_stats_8day.py
    MXD13Q1 8-days statistics data, long-term average, max, min and stdev
DESCRIPTION
    Input data for this script will use MXD13Q1 8-days data generate from GEE or downloaded from NASA
    This script can do 8-days 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 02_mxd13q1_stats_8day.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 = "X:\\Temp\\modis\\{}\\gee\\02_positive\\all".format(iso3)
output_folder = "X:\\Temp\\modis\\{}\\gee\\03_statistics\\cropland_8day".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
    # Only ask the user if a choice hasn't been made yet
    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 is_leap_year(year):
    return calendar.isleap(year)

def adjust_doy_for_leap_year(year, doy):
    if is_leap_year(year) and doy >= 60:
        # For leap years, after Feb 28 (DOY 59), adjust DOY by -1 to match non-leap year DOY
        return doy - 1
    return doy

def date_from_doy(year, doy):
    return datetime(year, 1, 1) + timedelta(doy - 1)

def process_files(input_folder, output_folder):
    global user_choice

    groups = defaultdict(list)

    for file in os.listdir(input_folder):
        if file.endswith((".tif", ".tiff")):
            parts = file.split("_")
            date_str = parts[-1].split(".")[0]  # Extract date string
            
            # Detect date format and adjust accordingly
            if len(date_str) == 7:  # yyyydoy format
                year, doy = int(date_str[:4]), int(date_str[4:])
            elif len(date_str) == 8:  # yyyymmdd format
                year, month, day = int(date_str[:4]), int(date_str[4:6]), int(date_str[6:])
                date = datetime(year, month, day)
                doy = date.timetuple().tm_yday
            
            # Adjust DOY for leap years to match non-leap year DOYs
            doy = adjust_doy_for_leap_year(year, doy)
            
            # Skip years outside the specified range
            if not start_year <= year <= end_year:
                continue
            
            groups[doy].append(os.path.join(input_folder, file))
    

    # 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.
    stat_names = {"MAXIMUM": "max", "MINIMUM": "min", "MEAN": "avg", "MEDIAN": "med", "STD": "std"}

    # Process each DOY group
    for doy, files in groups.items():
        for stat_type, suffix in stat_names.items():
            out_raster_name = f"{iso3}_phy_mxd13q1_20yr_{suffix}_{str(doy).zfill(3)}.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(files, stat_type, "DATA")
            outCellStatistics.save(temp_file_path)
            
            # Use CopyRaster to save the output raster with specific pixel type and LZW 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 in-memory raster
            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_files() function for the input folder
	process_files(input_folder, output_folder)

if __name__ == '__main__':
	main()

print("Script completed.")