MXD13Q1 quarter mean (JFM, AMJ, JAS, OND) from 8-day data

modis_quarter.py

# -*- coding: utf-8 -*-
"""
NAME
    modis_quarter.py
    MXD13Q1 quarter mean (JFM, AMJ, JAS, OND) from 8-day data
DESCRIPTION
    Input data for this script will use MXD13Q1 8-days data generate from GEE or downloaded from NASA
    This script can do MEAN calculation to derive quarter data (JFM, AMJ, JAS, OND)
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_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
from datetime import datetime

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

# ISO3 Country Code
iso3 = "syr" # Syria: syr, Myanmar: mmr, Lebanon: lbn

# Define the range of years
start_year = 2002
end_year = 2023

# Change the data and output folder
input_folder = "X:\\Temp\\modis\\{}\\gee\\02_positive\\cropland".format(iso3)
quarterly_output_folder = "X:\\Temp\\modis\\{}\\gee\\02_positive\\cropland_quarterly".format(iso3)

# Define the quarters
quarters = {
    "JFM": [(1, 1), (3, 23)],
    "AMJ": [(3, 28), (6, 27)],
    "JAS": [(7, 1), (9, 23)],
    "OND": [(9, 27), (12, 28)]
}

def date_in_range(date, start, end):
    return start <= date <= end

# Initialize the data structure for holding file paths
quarterly_groups = {year: {quarter: [] for quarter in quarters} for year in range(start_year, end_year + 1)}

for file in os.listdir(input_folder):
    if file.endswith(".tif") or file.endswith(".tiff"):
        parts = file.split("_")
        date_str = parts[-1].split(".")[0]
        date = datetime.strptime(date_str, "%Y%m%d")

        if not (start_year <= date.year <= end_year):
            continue

        for quarter, (start_month_day, end_month_day) in quarters.items():
            start_date = datetime(date.year, *start_month_day)
            end_date = datetime(date.year, *end_month_day)
            if date_in_range(date, start_date, end_date):
                quarterly_groups[date.year][quarter].append(os.path.join(input_folder, file))
                break

for year, quarters_files in quarterly_groups.items():
    for quarter, files in quarters_files.items():
        if files:
            out_raster_name = f"{iso3}_phy_mxd13q1_{quarter}_{year}_mean.tif"
            out_raster_path = os.path.join(quarterly_output_folder, out_raster_name)
            
            if not arcpy.Exists(out_raster_path):
                arcpy.CheckOutExtension("spatial")
                outCellStatistics = arcpy.sa.CellStatistics(files, "MEAN", "DATA")
                outCellStatistics.save(out_raster_path)
                arcpy.CheckInExtension("spatial")
                print(f"{out_raster_name} completed")
            else:
                print(f"{out_raster_name} exists")