Generate derivative product from Vegetation Indices per quarter

modis_viproducts_quarter.py

# -*- coding: utf-8 -*-
"""
NAME
    modis_viproducts_quarter.py
    Generate derivative product from Vegetation Indices
DESCRIPTION
    Input data for this script will use MXD13Q1 8-days data generate from GEE or downloaded 
    from NASA. This script can do calculation for ratio, difference, standardize anomaly 
    and vegetation condition index.
    The calculation required timeseries VI and the long-term statistics (min, mean, max, 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 modis_viproducts_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 os
import arcpy

# To avoid overwriting outputs, set this option to False.
arcpy.env.overwriteOutput = True

# ISO3 Country Code
iso3 = "syr"  # Update this to the country code you're working with

# Define input and output folders
input_folder = "X:\\Temp\\modis\\{}\\gee\\02_positive\\cropland_quarterly".format(iso3)
stats_folder = "X:\\Temp\\modis\\{}\\gee\\03_statistics\\cropland_quarterly_2012_2022".format(iso3)
ratioanom_folder = "X:\\Temp\\modis\\{}\\gee\\05_ratioanom\\cropland_quarterly_2012_2022".format(iso3)
diffanom_folder = "X:\\Temp\\modis\\{}\\gee\\06_diffanom\\cropland_quarterly_2012_2022".format(iso3)
stdanom_folder = "X:\\Temp\\modis\\{}\\gee\\07_stdanom\\cropland_quarterly_2012_2022".format(iso3)
vci_folder = "X:\\Temp\\modis\\{}\\gee\\08_vci\\cropland_quarterly_2012_2022".format(iso3)

# Quarters
quarters = ["JFM", "AMJ", "JAS", "OND"]

print("Script started...")

# Check if input folder exists
if not os.path.exists(input_folder):
    print(f"Input folder does not exist: {input_folder}")
    exit(1)

print("Input folder found, processing...")

def derivative_vi(input_folder, stats_folder, ratioanom_folder, diffanom_folder, stdanom_folder, vci_folder):
    # Loop through the files in the input folder
    for raster in os.listdir(input_folder):
        print(f"Processing raster: {raster}")
        if raster.endswith(".tif"):  # Only process .tif files
            print(f"Raster is a TIFF: {raster}")
            parts = raster.split("_")
            print(f"Filename parts: {parts}")  # Check the parts of the filename
            if len(parts) == 6:  # Adjusting to the actual number of parts in the filenames
                quarter = parts[3]
                year = parts[4]  # The year is now correctly the 5th part (index 4)
                print(f"Quarter: {quarter}, Year: {year}")
                if quarter in quarters:
                    print(f"Quarter is valid: {quarter}")
                    # Construct the corresponding filenames in the stats folder
                    avg_raster_name = f"{iso3}_phy_mxd13q1_{quarter}_2012_2022_avg.tif"
                    min_raster_name = f"{iso3}_phy_mxd13q1_{quarter}_2012_2022_min.tif"
                    max_raster_name = f"{iso3}_phy_mxd13q1_{quarter}_2012_2022_max.tif"
                    std_raster_name = f"{iso3}_phy_mxd13q1_{quarter}_2012_2022_std.tif"

                    print(f"Expected stats raster names: {avg_raster_name}, {min_raster_name}, {max_raster_name}, {std_raster_name}")

                    # Full paths to the statistics rasters
                    avg_raster = os.path.join(stats_folder, avg_raster_name)
                    min_raster = os.path.join(stats_folder, min_raster_name)
                    max_raster = os.path.join(stats_folder, max_raster_name)
                    std_raster = os.path.join(stats_folder, std_raster_name)

                    # Verify that the statistics rasters exist
                    stats_rasters_exist = True
                    for stats_raster in [avg_raster, min_raster, max_raster, std_raster]:
                        if not arcpy.Exists(stats_raster):
                            print(f"Error: {stats_raster} not found in {stats_folder}")
                            stats_rasters_exist = False
                            break
                        else:
                            print(f"Stats raster found: {stats_raster}")

                    if not stats_rasters_exist:
                        continue

                    # Define output paths for anomalies and VCI
                    ratioanom_raster = os.path.join(ratioanom_folder, f"{iso3}_phy_mxd13q1_ratioanom_{quarter}_{year}.tif")
                    diffanom_raster = os.path.join(diffanom_folder, f"{iso3}_phy_mxd13q1_diffanom_{quarter}_{year}.tif")
                    stdanom_raster = os.path.join(stdanom_folder, f"{iso3}_phy_mxd13q1_stdanom_{quarter}_{year}.tif")
                    vci_raster = os.path.join(vci_folder, f"{iso3}_phy_mxd13q1_vci_{quarter}_{year}.tif")

                    try:
                        arcpy.CheckOutExtension("spatial")

                        # Prepare the input raster
                        in_raster = arcpy.Raster(os.path.join(input_folder, raster))

                        # Calculate the indices
                        ratioanom = arcpy.sa.Divide(arcpy.sa.Float(in_raster), arcpy.sa.Float(avg_raster)) * 100
                        diffanom = arcpy.sa.Float(in_raster) - arcpy.sa.Float(avg_raster)
                        stdanom = arcpy.sa.Divide((arcpy.sa.Float(in_raster) - arcpy.sa.Float(avg_raster)), arcpy.sa.Float(std_raster))
                        vci_unclipped = arcpy.sa.Divide((arcpy.sa.Float(in_raster) - arcpy.sa.Float(min_raster)), (arcpy.sa.Float(max_raster) - arcpy.sa.Float(min_raster))) * 100
                        vci = arcpy.sa.Con(vci_unclipped > 100, 100, arcpy.sa.Con(vci_unclipped < 0, 0, vci_unclipped))

                        # Save the output rasters
                        ratioanom.save(ratioanom_raster)
                        print(f"{ratioanom_raster} completed, calculated from {raster} and {avg_raster_name}")

                        diffanom.save(diffanom_raster)
                        print(f"{diffanom_raster} completed, calculated from {raster} and {avg_raster_name}")

                        stdanom.save(stdanom_raster)
                        print(f"{stdanom_raster} completed, calculated from {raster}, {avg_raster_name} and {std_raster_name}")

                        vci.save(vci_raster)
                        print(f"{vci_raster} completed, calculated from {raster}, {min_raster_name} and {max_raster_name}")
                    
                    except Exception as e:
                        print(f"An error occurred: {e}")
                    finally:
                        arcpy.CheckInExtension("spatial")

            else:
                print(f"Filename does not have 5 parts: {raster}")
        else:
            print(f"Skipping non-TIFF file: {raster}")

# Main function
def main():
    # Call the derivative_vi() function for the input folder
    derivative_vi(input_folder, stats_folder, ratioanom_folder, diffanom_folder, stdanom_folder, vci_folder)

if __name__ == '__main__':
    main()

print("Script completed.")