diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/process_tile_array.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/process_tile_array.py new file mode 100644 index 000000000..64aea765a --- /dev/null +++ b/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) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/readme.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/readme.txt new file mode 100644 index 000000000..b1ff2c6db --- /dev/null +++ b/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 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/stitch_tiles.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/stitch_tiles.py new file mode 100644 index 000000000..4e07675a0 --- /dev/null +++ b/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) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/submit_global_array.sh b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/submit_global_array.sh new file mode 100644 index 000000000..e354b0018 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/veg/lai_grn/MCD15A3H_smooth/submit_global_array.sh @@ -0,0 +1,109 @@ +#!/bin/bash +#SBATCH --job-name=MODIS_Global_10x10 +#SBATCH --output=logs/global_%A_%a.out # %A is master job ID, %a is array ID +#SBATCH --error=logs/global_%A_%a.err +#SBATCH --time=04:00:00 # Max time per node +#SBATCH --qos=allnccs +#SBATCH --nodes=1 # 1 node per array task +#SBATCH --ntasks=16 # 16 cores per node +#SBATCH --account=s1583 + +##Please run the script by: +##source g5_modules.sh +##mkdir -p logs temporary output +##sbatch --array=0-347%20 submit_global_array.sh +##when all the jobs are done, then: +##sbatch --array=348-647%20 submit_global_array.sh + +export OMP_NUM_THREADS=1 +export OPENBLAS_NUM_THREADS=1 +export MKL_NUM_THREADS=1 +export HDF5_USE_FILE_LOCKING=FALSE + +MODIS_DATA_DIR="/discover/nobackup/projects/lis/LS_PARAMETERS/MODIS/MCD15A3H.061" + +if [ -z "$SLURM_ARRAY_TASK_ID" ]; then + echo "FATAL ERROR: Missing array parameter!" + echo "You MUST submit this script with an array range, for example:" + echo "sbatch --array=0-347%20 submit_global_array.sh" + echo "sbatch --array=348-647%20 submit_global_array.sh" + exit 1 +fi + +mkdir -p logs +mkdir -p temporary +mkdir -p output + +# ============================================================================== +# [GLOBAL GRID MAPPING] Convert Array ID (0-647) to Lat & Lon +# Grid size: 10x10. Lat range: -90 to 90 (18 rows). Lon range: -180 to 180 (36 cols). +# ============================================================================== +LAT_IDX=$((SLURM_ARRAY_TASK_ID / 36)) +LON_IDX=$((SLURM_ARRAY_TASK_ID % 36)) + +# Calculate numerical values using bc (Starts from Bottom-Left: Lat -90, Lon -180) +CALC_LAT=$(echo "-90 + $LAT_IDX * 10" | bc) +CALC_LON=$(echo "-180 + $LON_IDX * 10" | bc) + +# Format to ensure .0 suffix (e.g., -90 -> -90.0) +BASE_LAT=$(printf "%.1f" $CALC_LAT) +BASE_LON=$(printf "%.1f" $CALC_LON) + +echo "==========================================================" +echo "Array Task ID: $SLURM_ARRAY_TASK_ID" +echo "Node assigned: $SLURM_NODELIST" +echo "Processing Global Grid at Base Lat: $BASE_LAT, Lon: $BASE_LON" +echo "==========================================================" + +# ============================================================================== +# BATCH 1: Execute sub-tasks 0 to 7 +# ============================================================================== +echo "--- Launching Batch 1 (Tasks 0 to 7) ---" +for task_id in {0..7} +do + python3 process_tile_array.py --task_id $task_id --base_lat $BASE_LAT --base_lon $BASE_LON --data_dir "$MODIS_DATA_DIR" & +done +wait +echo "Batch 1 completed successfully!" + +# ============================================================================== +# BATCH 2: Execute sub-tasks 8 to 15 +# ============================================================================== +echo "--- Launching Batch 2 (Tasks 8 to 15) ---" +for task_id in {8..15} +do + python3 process_tile_array.py --task_id $task_id --base_lat $BASE_LAT --base_lon $BASE_LON --data_dir "$MODIS_DATA_DIR" & +done +wait +echo "Batch 2 completed successfully!" + +# ============================================================================== +# STITCHING & COMPRESSION +# ============================================================================== +echo "--- Initiating Stitching and Compression Process ---" +python3 stitch_tiles.py --base_lat $BASE_LAT --base_lon $BASE_LON + +if [ $? -eq 0 ]; then + echo "Stitching script executed successfully!" + echo "Cleaning up 16 uncompressed temporary sub-tiles..." + + for row in {0..3}; do + for col in {0..3}; do + LAT_MIN=$(echo "$BASE_LAT + $row * 2.5" | bc) + LON_MIN=$(echo "$BASE_LON + $col * 2.5" | bc) + LAT_MIN_FMT=$(printf "%.1f" $LAT_MIN) + LON_MIN_FMT=$(printf "%.1f" $LON_MIN) + rm -f "temporary/MCD15A3H_2003_2025_lat_${LAT_MIN_FMT}_lon_${LON_MIN_FMT}.nc" + done + done + echo "Cleanup complete!" +else + echo "ERROR: Stitching failed! Temporary files preserved for debugging." + exit 1 +fi + + +echo "==========================================================" +echo "Global Task ID $SLURM_ARRAY_TASK_ID (Lat: $BASE_LAT, Lon: $BASE_LON) finished!" +echo "==========================================================" +