Filling null on MXD13Q1 timeseries data with long-term mean

03_mxd13q1_fillnullwithstats.py

# -*- coding: utf-8 -*-
"""
NAME
    03_mxd13q1_fillnullwithstats.py
    Filling null on MXD13Q1 data with long-term mean
DESCRIPTION
    Input data for this script will use MXD13Q1 16-days data generate from GEE or downloaded from NASA
    This script can do filling null on MXD13Q1 data with long-term mean
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 03_mxd13q1_fillnullwithstats.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 os
import arcpy
from datetime import datetime, timedelta
import re
import calendar
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"  # Update this to your ISO3 Country Code

# Define input and output folders
input_folder = f"D:\\temp\\modis\\{iso3}\\gee\\02_positive\\evi_all"
stats_folder = f"D:\\temp\\modis\\{iso3}\\gee\\03_statistics\\evi_all_8day"
fnws_folder = f"D:\\temp\\modis\\{iso3}\\gee\\04_fillnullwithstats\\evi_all_8day"

# 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 process_raster(match, raster, input_folder, stats_folder, fnws_folder):
    global user_choice
    
    year = int(match.group(1))
    if len(match.group(2)) == 3:  # yyyydoy format
        doy = int(match.group(2))
        date = datetime(year, 1, 1) + timedelta(doy - 1)
    else:  # yyyymmdd format
        month = int(match.group(2))
        day = int(match.group(3))
        date = datetime(year, month, day)

    doy = date.strftime('%j')
    out_raster_name = f"{iso3}_phy_mxd13q1_evi_{date.strftime('%Y%m%d')}.tif"
    out_raster_path = os.path.join(fnws_folder, out_raster_name)

    if os.path.exists(out_raster_path):
        if user_choice is None:
            set_user_decision()
        if user_choice == 'S':
            print(f"Skipping existing file: {out_raster_path}")
            return
        elif user_choice == 'A':
            print("Aborting process.")
            exit()
    else:
        print(f"Processing {raster}...")

    stats_raster_name = f"{iso3}_phy_mxd13q1_20yr_avg_{doy.zfill(3)}.tif"
    stats_raster_path = os.path.join(stats_folder, stats_raster_name)

    if not arcpy.Exists(stats_raster_path):
        print(f"Error: {stats_raster_name} not found in {stats_folder}")
        return

    # 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)
    
    # Fill null values with the corresponding values from the stats raster
    arcpy.CheckOutExtension("spatial")
    in_raster = arcpy.Raster(os.path.join(input_folder, raster))
    stats_raster = arcpy.Raster(stats_raster_path)
    out_raster = arcpy.sa.Con(arcpy.sa.IsNull(in_raster), stats_raster, in_raster)

    # Save the output raster with LZW compression
    out_raster.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.")


def fill_null_with_stats(input_folder, stats_folder, fnws_folder):
    global user_choice
    
    # Check if output folder exists, create if not
    if not os.path.exists(fnws_folder):
        os.makedirs(fnws_folder)

    # Loop through the files in the input folder
    for raster in os.listdir(input_folder):
        if raster.endswith(".tif") or raster.endswith(".tiff"):
            match = re.search(r"(\d{4})(\d{2})(\d{2})\.tif$", raster) or re.search(r"(\d{4})(\d{3})\.tif$", raster)
            if match:
                process_raster(match, raster, input_folder, stats_folder, fnws_folder)
            else:
                print(f"Filename does not match expected formats: {raster}")

# Main function
def main():
    fill_null_with_stats(input_folder, stats_folder, fnws_folder)

if __name__ == '__main__':
    main()

print("Script completed.")