Add global MODIS LAI (MCD15A3H) smoothing pipeline in preprocessing

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/process_tile_array.py

@@ -0,0 +1,199 @@
+import xarray as xr
+import numpy as np
+import pandas as pd
+import glob
+import os
+import argparse
+from scipy.signal import savgol_filter
+import sys
+
+# ==============================================================================
+# 1. Core Algorithm & Preprocess Functions
+# ==============================================================================
+
+def timesat_upper_envelope_sg_nd(da, window_size=15, poly_order=2, iterations=3):
+    """
+    Applies an upper-envelope Savitzky-Golay filter to time-series data.
+    This mimics the TIMESAT algorithm (Jönsson and Eklundh 2004), 
+    which assumes negative biases in NDVI/LAI (due to clouds/aerosols) 
+    and iteratively fits the upper envelope to recover the true vegetation peak.
+
+    Jönsson, P., & Eklundh, L. (2004). TIMESAT—a program for analyzing 
+    time-series of satellite sensor data. Computers & geosciences, 30(8), 833-845.
+    """
+    y = da.values 
+    y_fit = np.copy(y)
+
+    # Iteratively apply SG filter and force the fitted curve towards the upper envelope
+    for _ in range(iterations):
+        y_smooth = savgol_filter(y_fit, window_size, poly_order, axis=0)
+        y_fit = np.maximum(y, y_smooth)
+
+    # Final smoothing step
+    y_final = savgol_filter(y_fit, window_size, poly_order, axis=0)
+
+    # Clip physically impossible values (LAI bounded between 0.0 and 8.0 for this context)
+    y_clipped = np.clip(y_final, 0.0, 8.0).astype(np.float32)
+
+    return xr.DataArray(y_clipped, dims=da.dims, coords=da.coords)
+
+def preprocess_modis(ds):
+    """
+    Extracts the date from the MODIS HDF/NetCDF filename and adds a 'time' dimension.
+    Expected filename format: MCD15A3H.061_YYYYDDD_*.nc4 (where DDD is day of year).
+    """
+    fname = ds.encoding["source"]
+    date_str = os.path.basename(fname).split('_')[1].split('.')[0]
+    time_index = pd.to_datetime([date_str], format='%Y%j')
+    ds = ds.expand_dims(time=time_index)
+    return ds
+
+def process_modis_chunk(ds):
+    """
+    Parses MODIS Quality Control (QC) flags via bitwise operations and applies masking.
+    Retains only 'perfect' pixels or conditionally fills snow-covered pixels.
+    """
+    # Scale LAI (MODIS scale factor is 0.1) and exclude fill values (255)
+    lai = ds['Lai_500m'].where(ds['Lai_500m'] != 255) * 0.1
+
+    qc = ds['FparLai_QC'].where(ds['FparLai_QC'] != 255)
+    exqc = ds['FparExtra_QC'].where(ds['FparExtra_QC'] != 255)
+
+    qc_int = qc.fillna(255).astype(int)
+    exqc_int = exqc.fillna(255).astype(int)
+
+    # Bitwise decoding of MODIS QC flags
+    is_good_main = (qc_int & 1) == 0                   # Bit 0: Good quality
+    is_snow = ((exqc_int >> 2) & 1) == 1               # Bit 2: Snow/Ice
+    is_internal_cloud = ((exqc_int >> 5) & 1) == 1     # Bit 5: Internal Cloud
+    is_shadow = ((exqc_int >> 6) & 1) == 1             # Bit 6: Cloud Shadow
+
+    # A pixel is 'perfect' if it has good quality and no clouds, snow, or shadows
+    is_perfect = is_good_main & (~is_internal_cloud) & (~is_snow) & (~is_shadow)
+
+    # Calculate a baseline for filling snow pixels (minimum valid LAI over time)
+    perfect_lai = lai.where(is_perfect & (lai > 0.05))
+    p_fillbase = perfect_lai.min(dim='time')
+    p_fillbase = p_fillbase.fillna(0.0) # Ensure water bodies get 0.0 base, not 0.1
+
+    # Apply baseline to snow pixels, keep original LAI otherwise
+    lai_filtered = xr.where(is_snow, p_fillbase, lai)
+
+    # Mask out bad pixels (not perfect and not snow)
+    bad_mask = (~is_perfect) & (~is_snow)
+    lai_filtered = lai_filtered.where(~bad_mask)
+
+    # Cap anomalously high LAI values
+    lai_filtered = lai_filtered.where(lai_filtered <= 25.0)
+
+    return lai_filtered
+
+# ==============================================================================
+# 2. Main Execution Block with Dynamic Arguments
+# ==============================================================================
+if __name__ == "__main__":
+    from dask.diagnostics import ProgressBar
+
+    # Parse the task ID and base coordinates passed by the Bash script
+    parser = argparse.ArgumentParser(description="Process a 2.5x2.5 degree MODIS tile.")
+    parser.add_argument('--task_id', type=int, required=True, help="Slurm Array Task ID (0-15)")
+    parser.add_argument('--base_lat', type=float, required=True, help="Base latitude (bottom-left corner) of the 10x10 tile")
+    parser.add_argument('--base_lon', type=float, required=True, help="Base longitude (bottom-left corner) of the 10x10 tile")
+    parser.add_argument('--data_dir', type=str, required=True, help="Path to the MODIS MCD15A3H.061 raw data directory")
+    args = parser.parse_args()
+
+    # --------------------------------------------------------------------------
+    # Master 10x10 Grid Definition & Sub-tile Mapping
+    # --------------------------------------------------------------------------
+    base_lat = args.base_lat
+    base_lon = args.base_lon
+    step = 2.5
+
+    # Calculate row (lat) and col (lon) index based on task_id (0-15)
+    # Maps a 1D Slurm Array ID to a 4x4 2D spatial grid: row = id // 4, col = id % 4
+    row = args.task_id // 4
+    col = args.task_id % 4
+
+    lat_min = base_lat + row * step
+    lat_max = lat_min + step
+    lon_min = base_lon + col * step
+    lon_max = lon_min + step
+
+    out_file = f"temporary/MCD15A3H_2003_2025_lat_{lat_min}_lon_{lon_min}.nc"
+
+    print(f"--- Task ID {args.task_id} Started ---", flush=True)
+    print(f"Master Tile Base: Lat={base_lat}, Lon={base_lon}", flush=True)
+    print(f"Targeting Sub-tile: Lat [{lat_min} to {lat_max}], Lon [{lon_min} to {lon_max}]", flush=True)
+    print(f"Output will be saved to: {out_file}", flush=True)
+
+    # --------------------------------------------------------------------------
+    # Data Loading
+    # --------------------------------------------------------------------------
+    data_dir = args.data_dir
+
+    files = []
+    for year in range(2002, 2027):
+        file_pattern = os.path.join(data_dir, str(year), f"MCD15A3H.061_{year}*.nc4")
+        files.extend(sorted(glob.glob(file_pattern)))
+
+    ds = xr.open_mfdataset(
+        files, 
+        preprocess=preprocess_modis,
+        combine='by_coords', 
+        chunks="auto",
+        parallel=True
+    )
+
+    # Slice target region (handle both ascending and descending latitude coords safely)
+    lat_arr = ds['lat'].values
+    if lat_arr[0] > lat_arr[-1]:
+        ds_region = ds.sel(lat=slice(lat_max, lat_min), lon=slice(lon_min, lon_max))
+    else:
+        ds_region = ds.sel(lat=slice(lat_min, lat_max), lon=slice(lon_min, lon_max))
+
+    # ==========================================================================
+    # Dask Workflow: QC -> Interpolation -> Smoothing
+    # ==========================================================================
+    print("Applying QC Filtering...", flush=True)
+    lai_qc = process_modis_chunk(ds_region)
+
+    print("Interpolating missing values...", flush=True)
+    # Rechunk across space to ensure full time-series is available per pixel for interpolation
+    lai_qc = lai_qc.chunk({'time': -1, 'lat': 'auto', 'lon': 'auto'})
+
+    # Temporal interpolation and edge filling
+    lai_filled = lai_qc.interpolate_na(dim='time', method='linear')
+    lai_filled = lai_filled.bfill(dim='time').ffill(dim='time').fillna(0.0) 
+
+    print("Applying S-G Smoothing...", flush=True)
+    # Map the TIMESAT algorithm block-by-block via Dask
+    lai_smoothed = lai_filled.map_blocks(
+        timesat_upper_envelope_sg_nd,
+        kwargs={'window_size': 15, 'poly_order': 2, 'iterations': 3},
+        template=lai_filled
+    )
+
+    core_time_slice = slice('2003-01-01', '2025-12-31')    
+    # Truncate to the core study period
+    lai_smoothed_core = lai_smoothed.sel(time=core_time_slice)
+
+    lai_smoothed_core.name = 'Lai_500m'
+    lai_smoothed_core.attrs['units'] = "m^2/m^2"
+    lai_smoothed_core.attrs['long_name'] = f"QC-Filtered Upper-Envelope Smoothed LAI ({lat_min} to {lat_max}, {lon_min} to {lon_max})"
+
+    # ==========================================================================
+    # I/O Execution & Teardown Hang Prevention
+    # ==========================================================================
+    print("Triggering lazy computation into RAM...", flush=True)
+    # .load() forces Dask to execute the graph into memory BEFORE I/O.
+    # This prevents Lustre filesystem lock conflicts during concurrent writes.
+    lai_smoothed_core = lai_smoothed_core.load() 
+
+    print(f"Writing physical data to {out_file}...", flush=True)
+    lai_smoothed_core.to_dataset().to_netcdf(out_file)
+
+    print(f"Task {args.task_id} successfully saved! Forcing immediate exit to prevent Dask hang...", flush=True)
+
+    # Force OS-level exit to prevent Dask thread-pool from hanging during garbage collection
+    import os
+    os._exit(0)

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/readme.txt

@@ -0,0 +1,68 @@
+================================================================================
+MODIS LAI Global Tile Processing Package
+================================================================================
+Author: Yujin Zeng
+Email: yujin.zeng@nasa.gov
+Date: 08/25/2026
+================================================================================
+
+## 1. Description
+This package provides a high-performance, parallelized workflow to process 
+global MODIS LAI data (MCD15A3H) from 2003-01-01 to 2025-12-31. 
+The pipeline applies Quality Control (QC) filtering, temporal interpolation, 
+and Savitzky-Golay (TIMESAT-style) smoothing (Jönsson and Eklundh 2004). 
+
+It processes data into localized 2.5x2.5 degree chunks and subsequently stitches 
+them into standard 10x10 degree MODIS tiles with extreme ubyte + zlib compression 
+to optimize storage on the Discover supercomputer.
+
+Reference:
+Jönsson, P., & Eklundh, L. (2004). TIMESAT—a program for analyzing 
+time-series of satellite sensor data. Computers & geosciences, 30(8), 833-845.
+
+## 2. File Manifest
+This package contains three primary scripts:
+
+* process_tile_array.py
+  The core Python processing script. It handles the I/O, QC filtering, and 
+  smoothing for a specific 2.5x2.5 degree sub-tile. Intermediate NetCDF outputs 
+  generated by this script are temporarily saved into the `temporary/` directory.
+
+* stitch_tiles.py
+  The post-processing Python script. It loads the 16 corresponding 2.5-degree 
+  sub-tiles, merges them into a unified 10x10 degree master tile, applies level-1 
+  zlib compression, and writes the final NetCDF product to the `output/` directory.
+
+* submit_global_array.sh
+  The master Slurm array bash script that orchestrates the entire workflow. It 
+  maps the Slurm array ID to geographic coordinates, distributes the processing 
+  tasks, executes the stitching script, and performs automated cleanup.
+
+## 3. Directory Structure Prerequisites
+Before executing the pipeline, the following directory structure must be 
+initialized in your working directory to handle logs and I/O routing properly:
+
+* `logs/`      : Stores the Slurm output and error logs (.out/.err) for each node.
+* `temporary/` : Stores the uncompressed 2.5x2.5 degree intermediate sub-tiles.
+* `output/`    : Stores the final, highly compressed 10x10 degree merged tiles.
+
+Note: Upon the successful execution of `stitch_tiles.py` for a 10x10 grid, 
+the master bash script will automatically remove the 16 intermediate files from 
+the `temporary/` directory to conserve filesystem quotas.
+
+## 4. Usage Instructions
+Due to the Slurm QOS max submission limits (QOSMaxSubmitJobPerUserLimit) on 
+Discover, the 648 global grid tasks must be submitted in two separate batches.
+
+Please run the workflow by executing the following commands sequentially:
+
+# Step 1: Source g5_module and create the required directories
+source g5_modules.sh
+mkdir -p logs temporary output
+
+# Step 2: Submit the first batch (Tasks 0 to 347)
+sbatch --array=0-347%20 submit_global_array.sh
+
+# Step 3: Wait for the first batch to complete (monitor via squeue)
+# When all jobs from the first batch are done, submit the second batch:
+sbatch --array=348-647%20 submit_global_array.sh

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/stitch_tiles.py

@@ -0,0 +1,103 @@
+import xarray as xr
+import os
+import sys
+import argparse
+from dask.diagnostics import ProgressBar
+
+if __name__ == "__main__":
+    # ==========================================================================
+    # 1. Parse Arguments (base_lat, base_lon)
+    # ==========================================================================
+    parser = argparse.ArgumentParser(description="Stitch 16 sub-tiles into a 10x10 master tile.")
+    parser.add_argument('--base_lat', type=float, required=True, help="Base latitude (bottom-left corner)")
+    parser.add_argument('--base_lon', type=float, required=True, help="Base longitude (bottom-left corner)")
+    args = parser.parse_args()
+
+    base_lat = args.base_lat
+    base_lon = args.base_lon
+    step = 2.5
+
+    print(f"Initiating stitching for Master Tile starting at Lat: {base_lat}, Lon: {base_lon}", flush=True)
+
+    # ==========================================================================
+    # 2. Strict File Validation & Size Checking
+    # ==========================================================================
+    expected_files = []
+
+    # Generate the expected 16 filenames based on the 4x4 grid math
+    for row in range(4):
+        for col in range(4):
+            lat_min = base_lat + row * step
+            lon_min = base_lon + col * step
+            expected_files.append(f"temporary/MCD15A3H_2003_2025_lat_{lat_min}_lon_{lon_min}.nc")
+
+    # Set the minimum size threshold to 2.8 GB (2.8 * 1024^3 bytes)
+    MIN_SIZE_BYTES = 2.8 * 1024 * 1024 * 1024 
+
+    print(f"Verifying existence and size (>= 2.8 GB) of {len(expected_files)} required sub-tiles...", flush=True)
+
+    for f in expected_files:
+        # Check 1: Verify file existence
+        if not os.path.exists(f):
+            print(f"\n FATAL ERROR: Required file is missing: {f}")
+            print("Aborting stitch process to prevent incomplete master tile.")
+            sys.exit(1)
+
+        # Check 2: Verify file size meets the threshold
+        f_size_bytes = os.path.getsize(f)
+        f_size_gb = f_size_bytes / (1024 ** 3)
+        if f_size_gb < 2.8:
+            print(f"\n FATAL ERROR: File too small! {f}")
+            print(f"Expected >= 2.8 GB, but found {f_size_gb:.2f} GB.")
+            print("This usually indicates a killed or incomplete Slurm task. Aborting.")
+            sys.exit(1)
+
+    print("All 16 sub-tiles exist and pass the 2.8GB size check!", flush=True)
+
+    # ==========================================================================
+    # 3. Memory Loading & Stitching
+    # ==========================================================================
+    print("Loading data entirely into RAM (Bypassing Dask Write Locks)...", flush=True)
+    ds_merged = xr.open_mfdataset(
+        expected_files, 
+        combine='by_coords'
+    )
+
+    ds_merged.load()
+
+    # ==========================================================================
+    # Calculate Grid Indices (H01-H36, V01-V18) and Format Filename
+    # ==========================================================================
+    h_idx = int((base_lon + 180.0) / 10.0) + 1
+    v_idx = int((base_lat + 90.0) / 10.0) + 1
+
+    # Format with leading zeros (e.g., 1 -> 01)
+    out_file = f"output/MCD15A3H_lai_2003-2025.H{h_idx:02d}V{v_idx:02d}.nc"
+
+    ds_merged.attrs['Description'] = f"Merged 10x10 degree MODIS LAI smoothed tile (Tile: H{h_idx:02d}V{v_idx:02d}, Base: Lat {base_lat}, Lon {base_lon})."
+
+    # ==========================================================================
+    # 4. Extreme Compression Encoding (ubyte + zlib)
+    # ==========================================================================
+    encoding_dict = {
+        'Lai_500m': {
+            'dtype': 'uint8',            # Maps exactly to netCDF 'ubyte'
+            'scale_factor': 0.1,         
+            'add_offset': 0.0,           
+            '_FillValue': 255,           
+            'zlib': True,                
+            'complevel': 1               # Level 1 for ultra-fast run-length compression
+        },
+        'time': {
+            'dtype': 'int32'
+        },
+        'lat': {'dtype': 'float32'},
+        'lon': {'dtype': 'float32'}
+    }
+
+    print(f"Saving highly compressed (ubyte + zlib) dataset to {out_file}...", flush=True)
+
+    #with ProgressBar():
+    ds_merged.to_netcdf(out_file, encoding=encoding_dict)
+
+    print(f"Stitching Complete! Master tile saved as: {out_file}", flush=True)