Generate monthly derivative product from Vegetation Indices

modis_viproducts_monthly.py

# -*- coding: utf-8 -*-
"""
NAME
    modis_viproducts_monthly.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_monthly.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 = "tur"  # Update this to the country code you're working with

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

# Monthlys
monthlys = ["01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"]

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):
        if raster.endswith(".tif"):
            print(f"Processing raster: {raster}")
            parts = raster.split("_")
            if len(parts) == 6: # tur_phy_mxd13q1_evi_mean_202311.tif
                # Extract year and month from the last part of the filename
                date_part = parts[5].split(".")[0]
                year = date_part[:4]
                monthly = date_part[4:]

                if monthly in monthlys:
                    # Construct the corresponding filenames in the stats folder
                    avg_raster_name = f"{iso3}_phy_mxd13q1_{monthly}_2003_2022_avg.tif"
                    min_raster_name = f"{iso3}_phy_mxd13q1_{monthly}_2003_2022_min.tif"
                    max_raster_name = f"{iso3}_phy_mxd13q1_{monthly}_2003_2022_max.tif"
                    std_raster_name = f"{iso3}_phy_mxd13q1_{monthly}_2003_2022_std.tif"

                    # 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_{year}{monthly}.tif")
                    diffanom_raster = os.path.join(diffanom_folder, f"{iso3}_phy_mxd13q1_diffanom_{year}{monthly}.tif")
                    stdanom_raster = os.path.join(stdanom_folder, f"{iso3}_phy_mxd13q1_stdanom_{year}{monthly}.tif")
                    vci_raster = os.path.join(vci_folder, f"{iso3}_phy_mxd13q1_vci_{year}{monthly}.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"Invalid month in filename: {raster}")
            else:
                print(f"Unexpected filename format: {raster}")
        else:
            print(f"Skipping non-TIFF file: {raster}")

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.")