Extract positive value for EVI/NDVI data, exclude 0 values which associated with water/cloud

01_mxd13q1_extract_positive.py

# -*- coding: utf-8 -*-
"""
NAME
    01_mxd13q1_extract_positive.py
    Extract positive value for EVI/NDVI data, exclude 0 values which associated with water/cloud
DESCRIPTION
    Input data for this script will use MXD13Q1 16-days data generate from GEE or downloaded from NASA
    This script can do positive data extraction and keep the data free from 0 values.
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 01_mxd13q1_extract_positive.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

# 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"X:\\Temp\\modis\\{iso3}\\gee\\01_original"
bnd = f"X:\\Temp\\modis\\{iso3}\\cropland\\idn_bnd_subset_esa_cropland_250m_grid_diss.shp"
output_folder = f"X:\\Temp\\modis\\{iso3}\\gee\\02_positive\\all"

# 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_evi_files(input_folder, output_folder):
    # Check if output folder exists, create if not
    if not os.path.exists(output_folder):
        os.makedirs(output_folder)

    # Loop through each file in the input folder
    for file in os.listdir(input_folder):
        if file.endswith(".tif"):
            input_file_path = os.path.join(input_folder, file)
            output_file_path = os.path.join(output_folder, file)

            # Check if the output file already exists
            if os.path.exists(output_file_path):
                # Prompt for user decision when an existing file is encountered
                set_user_decision()
                if user_choice == 'S':
                    print(f"Skipping existing file: {output_file_path}")
                    continue
                elif user_choice == 'A':
                    print("Aborting process.")
                    exit()

            print(f"Processing file: {file}")

            # Read the input raster
            input_raster = arcpy.Raster(input_file_path)

            # Extract by attribute: keep pixels > 0
            query = "VALUE > 0"
            # Using in-memory workspace for temporary data
            temp_positive_path = "in_memory/positive"  
            positive = arcpy.sa.ExtractByAttributes(input_raster, query)
            # Save the intermediate raster to disk or in-memory workspace
            positive.save(temp_positive_path)  

            # Extract by mask using the boundary and directly save the output
            clip = arcpy.sa.ExtractByMask(temp_positive_path, bnd).save(output_file_path)
            
            # Delete the temporary "positive" raster to free up resources
            arcpy.Delete_management(temp_positive_path)
            
            print(f"Processed file saved: {output_file_path}")

            arcpy.CheckInExtension("spatial")

# Main function
def main():
    # Call the process_evi_files() function
    process_evi_files(input_folder, output_folder)

if __name__ == '__main__':
    main()

print("Script completed.")