diff --git a/README.md b/README.md index 23a9c4d..5030b95 100644 --- a/README.md +++ b/README.md @@ -1,159 +1,58 @@ - ### Automated Coastline Extraction for Erosion Modeling in Alaska - -The primary goal of this project is to enhance the accuracy of coastline extraction, particularly for erosion modeling in Deering, Alaska, using high-resolution Planet imagery with a 3-meter resolution. The project focuses on creating reliable ground truth data and labels that will be used to train the [DeepWaterMap algorithm](https://github.com/isikdogan/deepwatermap), a deep convolutional neural network designed to segment surface water on multispectral imagery. Originally trained on 30-meter resolution Landsat data, DeepWaterMap will be adapted to work with higher-resolution data in this project. - -One of the key challenges in coastline extraction is the application of the Normalized Difference Water Index (NDWI), a widely used remote sensing index for identifying water bodies. However, using a single threshold across an entire image often results in suboptimal accuracy. To address this, I implemented a sliding window approach combined with Otsu thresholding, which dynamically adjusts thresholds over localized regions of the image. This method has shown promising improvements in accuracy. - -The newly generated labeled data, derived from this approach, will be used to retrain the [DeepWaterMap algorithm](https://github.com/isikdogan/deepwatermap), replacing the original Global Surface Water data. This project aims to produce a more accurate and reliable tool for coastline detection, which is crucial for monitoring and mitigating coastal erosion in vulnerable areas like Alaska. - -## Installation - -### Prerequisites - -Before installing this project, ensure you have the following requirements: - -- **Python** - - **Project version**: Tested and developed with **Python 3.13.5** - - **Conda environment**: Recommended to use **Python 3.10** for best compatibility with dependencies - -- **Miniconda** (for managing conda environments) - -- **Git** (for cloning the repository) - -- **GDAL** (install via `conda-forge` for easier setup) - -- **Rasterio 1.4.3+** (for geospatial data processing) - ---- - -### Clone the Repository - -Clone the project using the `dev` branch (this branch contains the latest development features): - -```bash -git clone -b dev https://github.com/your-username/coastline-extraction.git -cd coastline-extraction -``` - ---- - -### Environment Setup - -1. **Create a virtual environment**: - -```bash -python -m venv coastline_env -source coastline_env/bin/activate # On Windows: coastline_env\Scripts\activate -``` - -2. **Install required dependencies**: - -```bash -# Core deep learning libraries -pip install torch torchvision - -# Geospatial data processing -pip install rasterio gdal - -# Data manipulation and visualization -pip install numpy pandas matplotlib - -# Image processing -pip install scikit-image opencv-python - -# Utilities -pip install tqdm pillow - -# Additional dependencies for data preprocessing -pip install shapely fiona geopandas -``` +# Automated Coastline Extraction for Erosion Modeling in Alaska + +Work of - Janvi Singh + +**Google Summer of Code (GSoC) 2026 Project** + +A comprehensive deep learning pipeline for extracting accurate coastlines from high-resolution PlanetLabs satellite imagery, specifically designed for coastal erosion monitoring in Arctic Alaska communities. + +## Overview + +The rapidly warming Arctic is leading to increased rates of coastal erosion, placing hundreds of Alaska communities at the frontline of climate change. This project provides an automated solution for extracting vectorized coastlines from 3-meter resolution PlanetLabs imagery to support erosion modeling and forecasting. + +## What has been implemented in this workspace +1. Unified pipeline entrypoint + - `coastline_pipeline.py` orchestrates steps: UDM masking, NDWI, DEM integration, quality flags, validation. +2. UDM/QA60 cloud/shadow/snow masking + - `data_preprocessing/udm_masking.py` supports PlanetLabs UDM2 bit flags and Sentinel-2 QA60. +3. Adaptive NDWI pipeline with windowing + majority voting + - `data_preprocessing/ndwi_with_udm.py` does local thresholding, mask application, and subscene vector output. +4. DEM/slope/cliff-aware terrain integration + - `data_preprocessing/dem_integration.py` computes slope/aspect/TRI and cliff masks, influences classification. +5. Shadow/artifact / quality flags + - `data_preprocessing/shadow_artifact_detection.py` generates 8-bit quality flags, filtering. +6. Data expansion pipeline + - `data_preprocessing/data_expansion.py` supports 2017-2026 ingestion, seasonal filtering, incremental catalog. +7. DeepWaterMap training and inference + - `training_pipeline/deepwatermap_train.py`, `train_unet.py`, `predict.py` support U-Net training with Planet imagery. +8. Validation framework + - `validation_framework.py` with transect RMSE, comparison to ground truth coastlines. + +## Status vs. “Potential areas of improvement” +| Area | Status | Notes / Action required | +|---|---|---| +| Improve training data with PlanetLabs UDM | in progress | UDM masking is implemented, verify dataset-level integration for training label creation in `deepwatermap_train.py` and `data_expansion.py`. +| Data expansion beyond 2017-2019 | ✓ done | `data_expansion.py` has 2020-2026 config and pipeline. +| Improved cliff area segmentation | ✓ done | `dem_integration.py` handles cliff detection (slope >30°, elevation-stratified thresholding). +| Handling shadows/buildings/artifacts | ✓ done | `shadow_artifact_detection.py` intended for this; check empirical results and refine thresholds. +| SWIR + elevation integration | partially done | elevation/DEM done; SWIR support likely not yet explicit (only RGBN). Add SWIR path if available. + + +### Key Features + +- **UDM-Based Cloud Masking**: Integrated QA60/UDM2 quality band processing for reliable cloud and shadow removal +- **Sliding Window NDWI**: Localized Otsu thresholding with majority voting for adaptive water detection +- **DEM Integration**: Terrain analysis for improved cliff/steep slope segmentation +- **Quality Flag System**: 8-bit comprehensive pixel-level quality assessment +- **DeepWaterMap Training**: U-Net based deep learning model adapted for 4-band PlanetLabs imagery +- **Data Expansion Pipeline**: Automated ingestion of imagery from 2016-2026 with seasonal filtering +- **Validation Framework**: Transect-based RMSE and raster-based accuracy metrics + +### Project Links + +- **Source Code**: https://github.com/fwitmer/CoastlineExtraction +- **My Fork**:https://github.com/janvis11/CoastlineExtraction +- **Discussion Forum**: https://github.com/fwitmer/CoastlineExtraction/discussions +- **Mentors**: Frank Witmer (fwitmer@alaska.edu), Ritika Kumari (rkjane333@gmail.com) --- - -## Configuration - -This project uses a centralized configuration system to manage file paths and parameters. -Configuration is handled through `config_template.json` and the `load_config.py` module. - -### Setting Up Configuration - -1. **Copy the template**: - -```bash -cp config_template.json config.json -``` - -2. **Edit the configuration**: Open `config.json` and modify the paths according to your setup: - -```json -{ - "data_dir": "data", - "image_folder": "sample_data/PlanetLabs", - "raw_data_folder": "raw_data", - "shapefile_folder": "USGS_Coastlines", - "ground_truth_folder": "ground_truth", - "processed_data_folder": "processed_data", - "training": { - "model_save_path": "training_pipeline/unet_model.pth", - "batch_size": 8, - "epochs": 30, - "learning_rate": 1e-4, - "image_size": [256, 256], - "train_split": 0.8, - "device": "auto" - } -} -``` - -### Using the Configuration System - -The `load_config.py` module provides convenient functions to access your data files: - -```python -from load_config import load_config, get_image_path, get_shapefile_path - -# Load configuration -config = load_config() - -# Get specific file paths -image_path = get_image_path(config, 0) # First image file -shapefile_path = get_shapefile_path(config, 0) # First shapefile -``` ---- - - -## Contributing - -### Working with the Dev Branch - -This project uses the `dev` branch for active development. When contributing: - -1. **Fork the repository on GitHub** - -2. **Clone your fork using the `dev` branch**: - -```bash -git clone -b dev https://github.com/your-username/coastline-extraction.git -cd coastline-extraction -``` - -3. **Create a feature branch from `dev`**: - -```bash -git checkout -b feature/your-feature-name -``` - -4. **Make your changes and commit them**: - -```bash -git add . -git commit -m "Add your feature description" -``` - -5. **Push to your fork**: - -```bash -git push origin feature/your-feature-name -``` - -6. **Create a Pull Request** targeting the `dev` branch (not `main`) diff --git a/data_preprocessing/dem_integration.py b/data_preprocessing/dem_integration.py new file mode 100644 index 0000000..48c814b --- /dev/null +++ b/data_preprocessing/dem_integration.py @@ -0,0 +1,621 @@ +""" +Module: dem_integration.py + +Description: Integrate Digital Elevation Model (DEM) data for enhanced coastline + segmentation, particularly in cliff areas and steep terrain. + + Features: + - DEM preprocessing and coregistration with satellite imagery + - Slope and aspect calculation + - Cliff/steep terrain detection + - Elevation-stratified NDWI thresholding + - Terrain-based mask generation for challenging areas + +Author: GSoC 2026 Team +Date: 2026-03-28 + +Issue: #102 - DEM Integration for Cliff Area Segmentation +""" + +import os +import numpy as np +import rasterio as rio +from rasterio.mask import mask +from rasterio.warp import calculate_default_transform, reproject, Resampling +from rasterio.io import MemoryFile +import cv2 +from scipy.ndimage import sobel +from matplotlib import pyplot as plt +import geopandas as gpd + +# Local imports +from load_config import load_config +from utils.check_crs import check_crs, crs_match + + +# ============================================================================= +# Configuration +# ============================================================================= + +# Slope thresholds for cliff detection (in degrees) +SLOPE_CONFIG = { + 'cliff_threshold': 30.0, # Slope angle >= 30 degrees = cliff + 'steep_threshold': 15.0, # Slope angle >= 15 degrees = steep terrain + 'moderate_threshold': 5.0, # Slope angle >= 5 degrees = moderate terrain +} + +# DEM processing parameters +DEM_CONFIG = { + 'target_crs': 'EPSG:32603', # UTM Zone 3N (matches Alaska coastline) + 'target_resolution': 3.125, # Match PlanetLabs resolution (meters) + 'dem_no_data': -9999, # Common DEM nodata value +} + + +# ============================================================================= +# DEM Preprocessing +# ============================================================================= + +def preprocess_dem(dem_path, target_bounds=None, target_transform=None, + target_crs='EPSG:32603', target_resolution=3.125): + """ + Preprocess DEM to match satellite imagery extent and resolution. + + Args: + dem_path (str): Path to input DEM GeoTIFF + target_bounds (tuple): Target bounds (left, bottom, right, top) + target_transform (Affine): Target geotransform + target_crs (str): Target CRS (default: EPSG:32603) + target_resolution (float): Target pixel resolution in meters + + Returns: + tuple: (dem_array, transform, profile) - Preprocessed DEM data + """ + print(f"Preprocessing DEM: {dem_path}") + + with rio.open(dem_path) as src: + dem_crs = src.crs + dem_transform = src.transform + dem_profile = src.profile.copy() + + # Read DEM data + dem_data = src.read(1).astype(np.float32) + + # Handle nodata + nodata = src.nodata if src.nodata is not None else DEM_CONFIG['dem_no_data'] + dem_data = np.where(dem_data == nodata, np.nan, dem_data) + + print(f" Original DEM shape: {dem_data.shape}") + print(f" Original CRS: {dem_crs}") + print(f" Elevation range: {np.nanmin(dem_data):.1f} - {np.nanmax(dem_data):.1f} m") + + # Reproject if needed + if str(dem_crs) != target_crs: + print(f" Reprojecting from {dem_crs} to {target_crs}...") + dem_data, dem_transform = reproject_dem( + dem_data, dem_transform, dem_crs, + target_crs, target_resolution + ) + + # Crop to target bounds if specified + if target_bounds is not None: + dem_data, dem_transform = crop_dem_to_bounds( + dem_data, dem_transform, target_bounds + ) + + # Create output profile + height, width = dem_data.shape + output_profile = { + 'driver': 'GTiff', + 'dtype': 'float32', + 'count': 1, + 'width': width, + 'height': height, + 'transform': dem_transform, + 'crs': target_crs, + 'nodata': np.nan, + } + + print(f" Preprocessed DEM shape: {dem_data.shape}") + + return dem_data, dem_transform, output_profile + + +def reproject_dem(dem_data, src_transform, src_crs, target_crs, target_resolution): + """ + Reproject DEM to target CRS and resolution. + + Args: + dem_data (np.ndarray): Source DEM array + src_transform (Affine): Source geotransform + src_crs (CRS): Source CRS + target_crs (str): Target CRS + target_resolution (float): Target resolution in meters + + Returns: + tuple: (reprojected_dem, new_transform) + """ + from rasterio.warp import calculate_default_transform, reproject + + # Calculate output dimensions + src_bounds = rio.transform.array_bounds( + dem_data.shape[0], dem_data.shape[1], src_transform + ) + + # Estimate output size based on target resolution + width = int((src_bounds[2] - src_bounds[0]) / target_resolution) + height = int((src_bounds[3] - src_bounds[1]) / target_resolution) + + # Calculate transform + dst_transform, dst_width, dst_height = calculate_default_transform( + src_crs, target_crs, width, height, *src_bounds + ) + + # Reproject + dst_data = np.empty((dst_height, dst_width), dtype=np.float32) + + reproject( + source=dem_data, + destination=dst_data, + src_transform=src_transform, + src_crs=src_crs, + dst_transform=dst_transform, + dst_crs=target_crs, + resampling=Resampling.bilinear, + src_nodata=np.nan, + dst_nodata=np.nan + ) + + return dst_data, dst_transform + + +def crop_dem_to_bounds(dem_data, dem_transform, target_bounds): + """ + Crop DEM to target bounds. + + Args: + dem_data (np.ndarray): DEM array + dem_transform (Affine): DEM geotransform + target_bounds (tuple): Target bounds (left, bottom, right, top) + + Returns: + tuple: (cropped_dem, new_transform) + """ + from rasterio.windows import from_bounds, transform + + # Calculate window + window = from_bounds(*target_bounds, dem_transform) + + # Round window to integer bounds + col_off = int(window.col_off) + row_off = int(window.row_off) + width = int(window.width) + height = int(window.height) + + # Ensure bounds are within array + col_off = max(0, col_off) + row_off = max(0, row_off) + width = min(width, dem_data.shape[1] - col_off) + height = min(height, dem_data.shape[0] - row_off) + + # Crop + cropped_dem = dem_data[row_off:row_off+height, col_off:col_off+width] + + # Update transform + new_transform = transform(window, dem_transform) + + return cropped_dem, new_transform + + +# ============================================================================= +# Terrain Analysis +# ============================================================================= + +def calculate_slope(dem_data, pixel_size=3.125): + """ + Calculate slope angle from DEM using Sobel operator. + + Args: + dem_data (np.ndarray): DEM elevation array + pixel_size (float): Pixel size in meters + + Returns: + np.ndarray: Slope angle in degrees + """ + # Calculate gradients using Sobel operator + dy, dx = sobel(dem_data, axis=0), sobel(dem_data, axis=1) + + # Convert to slope angle (degrees) + slope_radians = np.arctan(np.sqrt(dx*dx + dy*dy) / (2 * pixel_size)) + slope_degrees = np.degrees(slope_radians) + + return slope_degrees + + +def calculate_aspect(dem_data, pixel_size=3.125): + """ + Calculate aspect (direction of steepest slope) from DEM. + + Args: + dem_data (np.ndarray): DEM elevation array + pixel_size (float): Pixel size in meters + + Returns: + np.ndarray: Aspect angle in degrees (0-360, 0=North, clockwise) + """ + # Calculate gradients + dy, dx = sobel(dem_data, axis=0), sobel(dem_data, axis=1) + + # Calculate aspect + aspect_radians = np.arctan2(-dx, dy) + aspect_degrees = np.degrees(aspect_radians) + + # Convert to 0-360 range (0 = North, clockwise) + aspect_degrees = (aspect_degrees + 360) % 360 + + return aspect_degrees + + +def calculate_terrain_ruggedness(dem_data, window_size=3): + """ + Calculate Terrain Ruggedness Index (TRI). + + TRI = mean absolute difference between center pixel and neighbors + + Args: + dem_data (np.ndarray): DEM elevation array + window_size (int): Size of neighborhood window + + Returns: + np.ndarray: TRI values + """ + from scipy.ndimage import generic_filter + + def tri_calc(window): + center = window[len(window)//2] + return np.mean(np.abs(window - center)) + + tri = generic_filter(dem_data, tri_calc, size=window_size) + + return tri + + +def detect_cliff_areas(slope_data, aspect_data=None, tri_data=None): + """ + Detect cliff and steep terrain areas based on slope threshold. + + Args: + slope_data (np.ndarray): Slope angle array (degrees) + aspect_data (np.ndarray, optional): Aspect array + tri_data (np.ndarray, optional): Terrain ruggedness array + + Returns: + dict: Cliff detection results + """ + # Binary cliff mask + cliff_mask = slope_data >= SLOPE_CONFIG['cliff_threshold'] + steep_mask = slope_data >= SLOPE_CONFIG['steep_threshold'] + moderate_mask = slope_data >= SLOPE_CONFIG['moderate_threshold'] + + # Classify terrain + terrain_class = np.zeros_like(slope_data, dtype=np.uint8) + terrain_class[moderate_mask] = 1 # Moderate (5-15 degrees) + terrain_class[steep_mask] = 2 # Steep (15-30 degrees) + terrain_class[cliff_mask] = 3 # Cliff (>30 degrees) + + result = { + 'cliff_mask': cliff_mask, + 'steep_mask': steep_mask, + 'moderate_mask': moderate_mask, + 'terrain_class': terrain_class, + 'cliff_pixels': np.sum(cliff_mask), + 'steep_pixels': np.sum(steep_mask) - np.sum(cliff_mask), + 'moderate_pixels': np.sum(moderate_mask) - np.sum(steep_mask), + } + + if tri_data is not None: + result['tri'] = tri_data + result['high_ruggedness'] = tri_data > np.percentile(tri_data, 75) + + if aspect_data is not None: + result['aspect'] = aspect_data + + return result + + +# ============================================================================= +# Elevation-Stratified NDWI Thresholding +# ============================================================================= + +def elevation_stratified_threshold(ndwi_data, dem_data, slope_data, + base_threshold=None, n_elevation_zones=5): + """ + Apply elevation-stratified thresholding for NDWI classification. + + Different elevation zones may have different water characteristics + (e.g., shadows in valleys, snow at high elevations). + + Args: + ndwi_data (np.ndarray): NDWI array + dem_data (np.ndarray): DEM elevation array + slope_data (np.ndarray): Slope angle array + base_threshold (float, optional): Base NDWI threshold + n_elevation_zones (int): Number of elevation zones + + Returns: + tuple: (classified_array, zone_thresholds) + """ + if base_threshold is None: + # Use Otsu threshold as base + ndwi_8bit = ((ndwi_data + 1) * 127).astype(np.uint8) + _, base_threshold = cv2.threshold( + ndwi_8bit, 0, 1, cv2.THRESH_BINARY + cv2.THRESH_OTSU + ) + base_threshold = (base_threshold / 127) - 1 + + # Create elevation zones + dem_valid = dem_data[~np.isnan(dem_data)] + elevation_percentiles = np.percentile(dem_valid, np.linspace(0, 100, n_elevation_zones + 1)) + + # Create output array + classified = np.zeros_like(ndwi_data, dtype=np.uint8) + zone_thresholds = {} + + for i in range(n_elevation_zones): + # Create zone mask + zone_mask = (dem_data >= elevation_percentiles[i]) & \ + (dem_data < elevation_percentiles[i + 1]) + + # Adjust threshold based on slope + zone_slope = slope_data[zone_mask] + if len(zone_slope) > 0: + mean_slope = np.mean(zone_slope) + + # Steeper slopes need higher thresholds (shadows look like water) + slope_adjustment = 0.01 * (mean_slope / 10) # +0.01 per 10 degrees + + zone_threshold = min(base_threshold + slope_adjustment, 0.5) + else: + zone_threshold = base_threshold + + zone_thresholds[f'zone_{i}'] = { + 'elevation_range': (elevation_percentiles[i], elevation_percentiles[i + 1]), + 'threshold': zone_threshold, + 'slope_adjustment': slope_adjustment if len(zone_slope) > 0 else 0 + } + + # Apply threshold + zone_ndwi = ndwi_data[zone_mask] + zone_classified = zone_ndwi > zone_threshold + classified[zone_mask] = zone_classified.astype(np.uint8) + + return classified, zone_thresholds + + +def terrain_masked_ndwi(ndwi_data, slope_data, dem_data=None, + max_slope_for_water=45.0): + """ + Create terrain-masked NDWI classification. + + Water is unlikely to occur on very steep slopes - these are likely shadows. + + Args: + ndwi_data (np.ndarray): NDWI array + slope_data (np.ndarray): Slope angle array + dem_data (np.ndarray, optional): DEM for additional filtering + max_slope_for_water (float): Maximum slope where water can occur + + Returns: + np.ndarray: Terrain-masked water classification + """ + # Initial NDWI classification + ndwi_8bit = ((ndwi_data + 1) * 127).astype(np.uint8) + _, threshold = cv2.threshold( + ndwi_8bit, 0, 1, cv2.THRESH_BINARY + cv2.THRESH_OTSU + ) + ndwi_classified = ndwi_data > ((threshold / 127) - 1) + + # Create slope mask (exclude very steep areas) + slope_mask = slope_data < max_slope_for_water + + # Apply slope mask + terrain_masked = ndwi_classified & slope_mask + + return terrain_masked + + +# ============================================================================= +# Integration with NDWI Pipeline +# ============================================================================= + +def process_with_dem(image_path, dem_path, points_path=None, output_dir=None): + """ + Complete processing pipeline integrating DEM data with NDWI classification. + + Args: + image_path (str): Path to satellite imagery + dem_path (str): Path to DEM file + points_path (str, optional): Path to transect points + output_dir (str, optional): Output directory + + Returns: + dict: Processing results + """ + from ndwi_labels import get_ndwi_label + + if output_dir is None: + output_dir = "result_dem_integration" + os.makedirs(output_dir, exist_ok=True) + + print("=" * 60) + print("DEM-INTEGRATED COASTLINE EXTRACTION") + print("=" * 60) + + # Step 1: Preprocess DEM + print("\nStep 1: Preprocessing DEM...") + with rio.open(image_path) as src: + target_bounds = src.bounds + target_transform = src.transform + target_crs = str(src.crs) + pixel_size = abs(src.transform[0]) + + dem_data, dem_transform, dem_profile = preprocess_dem( + dem_path, target_bounds, target_transform, target_crs, pixel_size + ) + + # Step 2: Calculate terrain metrics + print("\nStep 2: Calculating terrain metrics...") + slope_data = calculate_slope(dem_data, pixel_size) + aspect_data = calculate_aspect(dem_data, pixel_size) + tri_data = calculate_terrain_ruggedness(dem_data) + + print(f" Slope range: {np.nanmin(slope_data):.1f} - {np.nanmax(slope_data):.1f} degrees") + print(f" Cliff pixels detected: {np.sum(slope_data >= SLOPE_CONFIG['cliff_threshold']):,}") + + # Step 3: Detect cliff areas + print("\nStep 3: Detecting cliff/steep terrain areas...") + terrain_results = detect_cliff_areas(slope_data, aspect_data, tri_data) + + # Step 4: Run standard NDWI + print("\nStep 4: Running NDWI classification...") + ndwi_results = get_ndwi_label(image_path, points_path, output_dir=output_dir) + + # Step 5: Apply terrain-based corrections + print("\nStep 5: Applying terrain-based corrections...") + ndwi_data = ndwi_results.get('ndwi', ndwi_results.get('ndwi_array')) + + if ndwi_data is not None: + # Terrain-masked classification + terrain_masked = terrain_masked_ndwi(ndwi_data, slope_data) + + # Elevation-stratified thresholding + stratified_class, zone_thresholds = elevation_stratified_threshold( + ndwi_data, dem_data, slope_data + ) + + print(f" Standard NDWI water pixels: {np.sum(ndwi_results.get('ndwi_classified', np.zeros_like(terrain_masked))):,}") + print(f" Terrain-masked water pixels: {np.sum(terrain_masked):,}") + print(f" Stratified water pixels: {np.sum(stratified_class):,}") + + # Step 6: Save outputs + print("\nStep 6: Saving outputs...") + + # Save DEM derivatives + base_name = os.path.splitext(os.path.basename(image_path))[0] + + # Save slope + slope_path = os.path.join(output_dir, f"{base_name}_slope.tif") + save_raster(slope_data, dem_profile, slope_path) + print(f" Saved slope: {slope_path}") + + # Save cliff mask + cliff_path = os.path.join(output_dir, f"{base_name}_cliff_mask.tif") + cliff_profile = dem_profile.copy() + cliff_profile.update(dtype='uint8') + save_raster(terrain_results['cliff_mask'].astype(np.uint8), cliff_profile, cliff_path) + print(f" Saved cliff mask: {cliff_path}") + + # Save visualization + save_dem_visualization( + dem_data, slope_data, terrain_results, ndwi_results, + output_dir, base_name + ) + + return { + 'dem_data': dem_data, + 'slope': slope_data, + 'aspect': aspect_data, + 'tri': tri_data, + 'terrain_results': terrain_results, + 'ndwi_results': ndwi_results, + 'zone_thresholds': zone_thresholds if 'zone_thresholds' in dir() else None + } + + +def save_raster(data, profile, output_path): + """Save array as GeoTIFF.""" + with rio.open(output_path, 'w', **profile) as dst: + dst.write(data, 1) + + +def save_dem_visualization(dem_data, slope_data, terrain_results, ndwi_results, + output_dir, base_name): + """ + Create visualization of DEM integration results. + """ + fig, axes = plt.subplots(2, 3, figsize=(20, 12)) + + # DEM + im = axes[0, 0].imshow(dem_data, cmap='terrain') + axes[0, 0].set_title('Digital Elevation Model') + axes[0, 0].axis('off') + plt.colorbar(im, ax=axes[0, 0], label='Elevation (m)') + + # Slope + im = axes[0, 1].imshow(slope_data, cmap='YlOrRd') + axes[0, 1].set_title(f'Slope (degrees)\nCliff pixels: {terrain_results["cliff_pixels"]:,}') + axes[0, 1].axis('off') + plt.colorbar(im, ax=axes[0, 1], label='Slope (°)') + + # Terrain classification + im = axes[0, 2].imshow(terrain_results['terrain_class'], cmap='viridis') + axes[0, 2].set_title('Terrain Classification') + axes[0, 2].axis('off') + + # NDWI + if 'ndwi' in ndwi_results: + im = axes[1, 0].imshow(ndwi_results['ndwi'], cmap='RdYlBu') + axes[1, 0].set_title('NDWI Values') + axes[1, 0].axis('off') + plt.colorbar(im, ax=axes[1, 0], label='NDWI') + + # NDWI classified + if 'ndwi_classified' in ndwi_results: + im = axes[1, 1].imshow(ndwi_results['ndwi_classified'], cmap='Blues') + axes[1, 1].set_title('NDWI Classified') + axes[1, 1].axis('off') + + # Cliff overlay + if 'ndwi_classified' in ndwi_results: + overlay = np.zeros((*ndwi_results['ndwi_classified'].shape, 3)) + overlay[ndwi_results['ndwi_classified'] == 1] = [0, 0, 1] # Water = blue + overlay[terrain_results['cliff_mask']] = [1, 0, 0] # Cliff = red + + im = axes[1, 2].imshow(overlay) + axes[1, 2].set_title('Water (blue) and Cliffs (red)') + axes[1, 2].axis('off') + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, f"{base_name}_dem_integration.png"), dpi=150) + plt.close() + + print(f" Saved visualization: {os.path.join(output_dir, f'{base_name}_dem_integration.png')}") + + +# ============================================================================= +# Main Execution +# ============================================================================= + +if __name__ == "__main__": + config = load_config() + + # Example usage + print("DEM Integration Module - Test Run") + print("=" * 60) + + # Get DEM path from environment or config + dem_path = os.environ.get('DEM_PATH', 'data/dem/Alaska_DEM.tif') + + if not os.path.exists(dem_path): + print(f"DEM file not found: {dem_path}") + print("Set DEM_PATH environment variable to test.") + else: + # Get image path from config + from load_config import get_image_path + image_path = get_image_path(config, 0) + + if image_path and os.path.exists(image_path): + results = process_with_dem(image_path, dem_path) + print("\nProcessing complete!") + else: + print("No image path available in config.") diff --git a/data_preprocessing/shadow_artifact_detection.py b/data_preprocessing/shadow_artifact_detection.py new file mode 100644 index 0000000..556d989 --- /dev/null +++ b/data_preprocessing/shadow_artifact_detection.py @@ -0,0 +1,644 @@ +""" +Module: shadow_artifact_detection.py + +Description: Enhanced detection and handling of shadows, cloud artifacts, + and satellite data quality issues beyond UDM/QA60 band reliance. + + Features: + - NIR-based shadow detection (independent of UDM) + - Satellite artifact detection (striping, sensor noise) + - Building/structure detection for urban areas + - Water shadow discrimination + - Quality flag generation for each pixel + +Author: GSoC 2026 Team +Date: 2026-03-28 + +Issue: #103 - Enhanced Shadow and Artifact Detection +""" + +import os +import numpy as np +import rasterio as rio +from rasterio.mask import mask +from rasterio.io import MemoryFile +import cv2 +from scipy import ndimage +from matplotlib import pyplot as plt +import geopandas as gpd + +# Local imports +from load_config import load_config + + +# ============================================================================= +# Configuration +# ============================================================================= + +# Shadow detection thresholds +SHADOW_CONFIG = { + 'nir_threshold': 0.15, # NIR reflectance threshold for shadow + 'blue_threshold': 0.10, # Blue band threshold + 'ndvi_threshold': 0.1, # NDVI threshold for shadow vs water + 'ratio_threshold': 0.5, # Blue/NIR ratio threshold + 'min_shadow_size': 9, # Minimum shadow patch size (pixels) +} + +# Artifact detection thresholds +ARTIFACT_CONFIG = { + 'stripe_threshold': 2.5, # Standard deviations for stripe detection + 'noise_threshold': 3.0, # Standard deviations for noise detection + 'saturation_threshold': 0.98, # Reflectance saturation level + 'edge_artifact_width': 3, # Width of edge artifact zone +} + +# Quality flags (bit positions) +QUALITY_FLAGS = { + 'CLEAR': 0, + 'SHADOW': 1, + 'CLOUD': 2, + 'ARTIFACT': 3, + 'SATURATED': 4, + 'EDGE': 5, + 'WATER_LIKELY': 6, + 'WATER_UNCERTAIN': 7, +} + + +# ============================================================================= +# Shadow Detection +# ============================================================================= + +def detect_shadows_nir_based(green_band, nir_band, red_band=None, blue_band=None): + """ + Detect shadows using NIR-based methods independent of UDM/QA60. + + Shadows have low NIR reflectance but different spectral characteristics + than water. + + Args: + green_band (np.ndarray): Green band reflectance + nir_band (np.ndarray): NIR band reflectance + red_band (np.ndarray, optional): Red band reflectance + blue_band (np.ndarray, optional): Blue band reflectance + + Returns: + tuple: (shadow_mask, shadow_confidence) + """ + print(" Detecting shadows using NIR-based method...") + + # Normalize bands to 0-1 range if needed + if np.max(nir_band) > 1: + green_band = green_band / np.max(green_band) + nir_band = nir_band / np.max(nir_band) + if red_band is not None: + red_band = red_band / np.max(red_band) + if blue_band is not None: + blue_band = blue_band / np.max(blue_band) + + # Method 1: Low NIR reflectance + low_nir = nir_band < SHADOW_CONFIG['nir_threshold'] + + # Method 2: Calculate NDVI (shadows have low/negative NDVI) + ndvi = (nir_band - red_band) / (nir_band + red_band + 1e-10) if red_band is not None else None + low_ndvi = ndvi < SHADOW_CONFIG['ndvi_threshold'] if ndvi is not None else np.ones_like(low_nir, dtype=bool) + + # Method 3: Blue/NIR ratio (shadows have ratio ~1, water has ratio < 1) + if blue_band is not None: + blue_nir_ratio = np.where(nir_band > 0.01, blue_band / (nir_band + 1e-10), 0) + shadow_ratio = np.abs(blue_nir_ratio - 1.0) < 0.3 # Shadows have similar blue and NIR + else: + shadow_ratio = np.ones_like(low_nir, dtype=bool) + + # Combine methods + shadow_mask = low_nir & low_ndvi & shadow_ratio + + # Remove small patches + shadow_mask = ndimage.binary_opening(shadow_mask, structure=np.ones((3, 3))) + shadow_mask = ndimage.binary_closing(shadow_mask, structure=np.ones((3, 3))) + + # Remove patches smaller than minimum size + labeled, n_features = ndimage.label(shadow_mask) + for i in range(1, n_features + 1): + if np.sum(labeled == i) < SHADOW_CONFIG['min_shadow_size']: + shadow_mask[labeled == i] = False + + # Calculate confidence + confidence = np.zeros_like(shadow_mask, dtype=np.float32) + confidence[low_nir] += 0.3 + confidence[low_ndvi] += 0.3 + confidence[shadow_ratio] += 0.4 + confidence = np.clip(confidence, 0, 1) + + shadow_pixels = np.sum(shadow_mask) + print(f" Shadow pixels detected: {shadow_pixels:,} ({100*shadow_pixels/shadow_mask.size:.1f}%)") + + return shadow_mask, confidence + + +def detect_cloud_shadows(ndwi_data, nir_band, shadow_mask=None): + """ + Distinguish cloud shadows from topographic shadows. + + Cloud shadows are typically adjacent to bright clouds and have + specific spatial patterns. + + Args: + ndwi_data (np.ndarray): NDWI array + nir_band (np.ndarray): NIR band reflectance + shadow_mask (np.ndarray, optional): Initial shadow mask + + Returns: + tuple: (cloud_shadow_mask, topographic_shadow_mask) + """ + if shadow_mask is None: + shadow_mask, _ = detect_shadows_nir_based(None, nir_band) + + # Detect bright areas (potential clouds) + bright_mask = nir_band > np.percentile(nir_band[~np.isnan(nir_band)], 90) + + # Dilate bright areas + bright_dilated = ndimage.binary_dilation(bright_mask, iterations=20) + + # Cloud shadows are typically near bright clouds + cloud_shadow_mask = shadow_mask & bright_dilated + + # Topographic shadows are isolated from clouds + topographic_shadow_mask = shadow_mask & ~bright_dilated + + print(f" Cloud shadow pixels: {np.sum(cloud_shadow_mask):,}") + print(f" Topographic shadow pixels: {np.sum(topographic_shadow_mask):,}") + + return cloud_shadow_mask, topographic_shadow_mask + + +def water_vs_shadow_discrimination(ndwi_data, nir_band, slope_data=None): + """ + Discriminate between water and shadow pixels. + + Both water and shadows have low NIR reflectance and can be confused. + Use multiple spectral indices and terrain information. + + Args: + ndwi_data (np.ndarray): NDWI array + nir_band (np.ndarray): NIR band reflectance + slope_data (np.ndarray, optional): Slope angle array + + Returns: + tuple: (water_mask, shadow_mask, uncertain_mask) + """ + # Initial classifications + water_candidate = ndwi_data > 0.0 # Positive NDWI + shadow_candidate = nir_band < SHADOW_CONFIG['nir_threshold'] + + # Overlap region (uncertain) + uncertain = water_candidate & shadow_candidate + + # Use slope if available (water unlikely on steep slopes) + if slope_data is not None: + steep = slope_data > 20 + uncertain[steep & uncertain] = True + water_candidate[steep] = False + + # Confident water (high NDWI, not shadow) + confident_water = water_candidate & ~shadow_candidate & (ndwi_data > 0.2) + + # Confident shadow (low NIR, low/neutral NDWI) + confident_shadow = shadow_candidate & (ndwi_data < 0.1) + + print(f" Confident water pixels: {np.sum(confident_water):,}") + print(f" Confident shadow pixels: {np.sum(confident_shadow):,}") + print(f" Uncertain pixels: {np.sum(uncertain):,}") + + return confident_water, confident_shadow, uncertain + + +# ============================================================================= +# Artifact Detection +# ============================================================================= + +def detect_satellite_striping(band_data, window_size=64): + """ + Detect striping artifacts common in satellite imagery. + + Striping appears as regular patterns of brighter/darker rows or columns. + + Args: + band_data (np.ndarray): Single band reflectance + window_size (int): Size of analysis window + + Returns: + tuple: (stripe_mask, stripe_direction) + """ + print(" Detecting satellite striping artifacts...") + + # Calculate local statistics + rows = band_data.shape[0] + cols = band_data.shape[1] + + # Row-wise mean (detect horizontal striping) + row_means = np.mean(band_data, axis=1) + row_std = np.std(row_means) + row_anomalies = np.abs(row_means - np.mean(row_means)) > ARTIFACT_CONFIG['stripe_threshold'] * row_std + + # Column-wise mean (detect vertical striping) + col_means = np.mean(band_data, axis=0) + col_std = np.std(col_means) + col_anomalies = np.abs(col_means - np.mean(col_means)) > ARTIFACT_CONFIG['stripe_threshold'] * col_std + + # Create stripe mask + stripe_mask = np.zeros_like(band_data, dtype=bool) + + horizontal_stripe_pixels = 0 + vertical_stripe_pixels = 0 + + if np.any(row_anomalies): + for i in range(rows): + if row_anomalies[i]: + stripe_mask[i, :] = True + horizontal_stripe_pixels += cols + + if np.any(col_anomalies): + for j in range(cols): + if col_anomalies[j]: + stripe_mask[:, j] = True + vertical_stripe_pixels += rows + + stripe_direction = 'none' + if horizontal_stripe_pixels > vertical_stripe_pixels: + stripe_direction = 'horizontal' + elif vertical_stripe_pixels > 0: + stripe_direction = 'vertical' + + total_stripe_pixels = np.sum(stripe_mask) + print(f" Stripe pixels detected: {total_stripe_pixels:,} ({100*total_stripe_pixels/stripe_mask.size:.1f}%)") + print(f" Stripe direction: {stripe_direction}") + + return stripe_mask, stripe_direction + + +def detect_sensor_noise(band_data, neighborhood_size=5): + """ + Detect sensor noise (salt-and-pepper, random spikes). + + Args: + band_data (np.ndarray): Single band reflectance + neighborhood_size (int): Size of neighborhood for local statistics + + Returns: + np.ndarray: Noise mask + """ + # Calculate local mean and std + kernel = np.ones((neighborhood_size, neighborhood_size)) / (neighborhood_size ** 2) + local_mean = ndimage.convolve(band_data, kernel, mode='reflect') + local_std = np.sqrt(ndimage.convolve((band_data - local_mean) ** 2, kernel, mode='reflect')) + + # Pixels that deviate significantly from local mean + deviation = np.abs(band_data - local_mean) / (local_std + 1e-10) + noise_mask = deviation > ARTIFACT_CONFIG['noise_threshold'] + + # Remove isolated pixels (likely real features) + noise_mask = ndimage.binary_opening(noise_mask, structure=np.ones((3, 3))) + + noise_pixels = np.sum(noise_mask) + print(f" Sensor noise pixels: {noise_pixels:,} ({100*noise_pixels/noise_mask.size:.1f}%)") + + return noise_mask + + +def detect_saturated_pixels(band_data, threshold=None): + """ + Detect saturated pixels (maximum sensor response). + + Args: + band_data (np.ndarray): Reflectance band + threshold (float, optional): Saturation threshold + + Returns: + np.ndarray: Saturation mask + """ + if threshold is None: + threshold = SHADOW_CONFIG['saturation_threshold'] + + # Normalize if needed + max_val = np.max(band_data) + if max_val > 1: + normalized = band_data / max_val + else: + normalized = band_data + + saturated = normalized > threshold + + sat_pixels = np.sum(saturated) + print(f" Saturated pixels: {sat_pixels:,} ({100*sat_pixels/saturated.size:.1f}%)") + + return saturated + + +def detect_edge_artifacts(band_data, edge_width=3): + """ + Detect edge artifacts from image processing (mosaicking, cropping). + + Args: + band_data (np.ndarray): Reflectance band + edge_width (int): Width of edge zone to check + + Returns: + np.ndarray: Edge artifact mask + """ + mask = np.zeros_like(band_data, dtype=bool) + + # Check image edges + mask[:edge_width, :] = True # Top edge + mask[-edge_width:, :] = True # Bottom edge + mask[:, :edge_width] = True # Left edge + mask[:, -edge_width:] = True # Right edge + + # Also check for internal sharp transitions (potential seam lines) + gradient_x = np.abs(np.gradient(band_data.astype(np.float32), axis=1)) + gradient_y = np.abs(np.gradient(band_data.astype(np.float32), axis=0)) + + # Find lines of high gradient (potential seams) + vertical_seam = np.any(gradient_x > np.percentile(gradient_x, 99), axis=0) + horizontal_seam = np.any(gradient_y > np.percentile(gradient_y, 99), axis=1) + + return mask + + +# ============================================================================= +# Quality Flag Generation +# ============================================================================= + +def generate_quality_flags(image_path, ndwi_data=None, slope_data=None): + """ + Generate comprehensive quality flags for each pixel. + + Args: + image_path (str): Path to multispectral GeoTIFF + ndwi_data (np.ndarray, optional): Pre-calculated NDWI + slope_data (np.ndarray, optional): Slope array for terrain masking + + Returns: + tuple: (quality_flags, flag_summary) + """ + print("\nGenerating comprehensive quality flags...") + + # Read spectral bands + with rio.open(image_path) as src: + blue = src.read(1).astype(np.float32) + green = src.read(2).astype(np.float32) + red = src.read(3).astype(np.float32) + nir = src.read(4).astype(np.float32) if src.count >= 4 else None + + if nir is None: + print(" Warning: NIR band not available, limited analysis possible") + nir = np.zeros_like(blue) + + # Normalize to reflectance (0-1) + for band in [blue, green, red, nir]: + max_val = np.nanmax(band) + if max_val > 1: + band /= max_val + + # Initialize quality flag array (8 bits) + quality = np.zeros(blue.shape, dtype=np.uint8) + + # Calculate NDWI if not provided + if ndwi_data is None: + ndwi_data = (green - nir) / (green + nir + 1e-10) + + # Flag 1: Shadow detection + shadow_mask, shadow_conf = detect_shadows_nir_based(green, nir, red, blue) + quality[shadow_mask] |= (1 << QUALITY_FLAGS['SHADOW']) + + # Flag 2: Cloud shadows + cloud_shadow, topo_shadow = detect_cloud_shadows(ndwi_data, nir, shadow_mask) + quality[cloud_shadow] |= (1 << QUALITY_FLAGS['SHADOW']) + + # Flag 3: Satellite striping + stripe_mask, stripe_dir = detect_satellite_striping(nir) + quality[stripe_mask] |= (1 << QUALITY_FLAGS['ARTIFACT']) + + # Flag 4: Sensor noise + noise_mask = detect_sensor_noise(nir) + quality[noise_mask] |= (1 << QUALITY_FLAGS['ARTIFACT']) + + # Flag 5: Saturated pixels + for band_name, band in [('Blue', blue), ('Green', green), ('Red', red), ('NIR', nir)]: + sat_mask = detect_saturated_pixels(band) + if np.any(sat_mask): + quality[sat_mask] |= (1 << QUALITY_FLAGS['SATURATED']) + + # Flag 6: Edge artifacts + edge_mask = detect_edge_artifacts(nir) + quality[edge_mask] |= (1 << QUALITY_FLAGS['EDGE']) + + # Flag 7 & 8: Water likelihood + water_likely, water_shadow, uncertain = water_vs_shadow_discrimination( + ndwi_data, nir, slope_data + ) + quality[water_likely] |= (1 << QUALITY_FLAGS['WATER_LIKELY']) + quality[uncertain] |= (1 << QUALITY_FLAGS['WATER_UNCERTAIN']) + + # Generate summary + flag_summary = { + 'shadow_pixels': np.sum((quality >> QUALITY_FLAGS['SHADOW']) & 1), + 'artifact_pixels': np.sum((quality >> QUALITY_FLAGS['ARTIFACT']) & 1), + 'saturated_pixels': np.sum((quality >> QUALITY_FLAGS['SATURATED']) & 1), + 'edge_pixels': np.sum((quality >> QUALITY_FLAGS['EDGE']) & 1), + 'water_likely_pixels': np.sum((quality >> QUALITY_FLAGS['WATER_LIKELY']) & 1), + 'water_uncertain_pixels': np.sum((quality >> QUALITY_FLAGS['WATER_UNCERTAIN']) & 1), + 'clear_pixels': np.sum(quality == 0), + } + + flag_summary['total_pixels'] = quality.size + flag_summary['flagged_percentage'] = 100 * (quality.size - flag_summary['clear_pixels']) / quality.size + + print(f"\nQuality Flag Summary:") + print(f" Clear pixels: {flag_summary['clear_pixels']:,} ({100-flag_summary['flagged_percentage']:.1f}%)") + print(f" Shadow pixels: {flag_summary['shadow_pixels']:,}") + print(f" Artifact pixels: {flag_summary['artifact_pixels']:,}") + print(f" Water likely: {flag_summary['water_likely_pixels']:,}") + print(f" Water uncertain: {flag_summary['water_uncertain_pixels']:,}") + + return quality, flag_summary + + +def apply_quality_filter(ndwi_data, quality_flags, max_shadow_fraction=0.3, + exclude_artifacts=True): + """ + Apply quality-based filtering to NDWI classification. + + Args: + ndwi_data (np.ndarray): NDWI array + quality_flags (np.ndarray): Quality flag array + max_shadow_fraction (float): Maximum shadow fraction for reliable classification + exclude_artifacts (bool): Whether to exclude artifact-affected pixels + + Returns: + tuple: (filtered_ndwi, reliability_mask) + """ + # Extract individual flags + shadow_flag = (quality_flags >> QUALITY_FLAGS['SHADOW']) & 1 + artifact_flag = (quality_flags >> QUALITY_FLAGS['ARTIFACT']) & 1 + uncertain_flag = (quality_flags >> QUALITY_FLAGS['WATER_UNCERTAIN']) & 1 + + # Reliability mask + reliable = ~shadow_flag.astype(bool) + if exclude_artifacts: + reliable = reliable & ~artifact_flag.astype(bool) + reliable = reliable & ~uncertain_flag.astype(bool) + + # Filter NDWI (set unreliable to NaN) + filtered_ndwi = ndwi_data.copy() + filtered_ndwi[~reliable] = np.nan + + return filtered_ndwi, reliable + + +# ============================================================================= +# Visualization +# ============================================================================= + +def save_quality_visualization(quality_flags, flag_summary, ndwi_data, + output_dir, base_name): + """ + Create visualization of quality flags. + """ + fig, axes = plt.subplots(2, 3, figsize=(20, 12)) + + # Extract individual masks + shadow_mask = (quality_flags >> QUALITY_FLAGS['SHADOW']) & 1 + artifact_mask = (quality_flags >> QUALITY_FLAGS['ARTIFACT']) & 1 + water_likely = (quality_flags >> QUALITY_FLAGS['WATER_LIKELY']) & 1 + water_uncertain = (quality_flags >> QUALITY_FLAGS['WATER_UNCERTAIN']) & 1 + + # Quality overview (RGB composite) + rgb = np.zeros((*quality_flags.shape, 3)) + rgb[:, :, 0] = shadow_mask # Red = shadow + rgb[:, :, 1] = artifact_mask # Green = artifact + rgb[:, :, 2] = water_likely # Blue = water likely + axes[0, 0].imshow(rgb) + axes[0, 0].set_title(f'Quality Overview\nR=Shadow, G=Artifact, B=Water Likely') + axes[0, 0].axis('off') + + # Shadow mask + axes[0, 1].imshow(shadow_mask, cmap='gray') + axes[0, 1].set_title(f'Shadow Mask\n{flag_summary["shadow_pixels"]:,} pixels') + axes[0, 1].axis('off') + + # Artifact mask + axes[0, 2].imshow(artifact_mask, cmap='gray') + axes[0, 2].set_title(f'Artifact Mask\n{flag_summary["artifact_pixels"]:,} pixels') + axes[0, 2].axis('off') + + # NDWI + im = axes[1, 0].imshow(ndwi_data, cmap='RdYlBu') + axes[1, 0].set_title('NDWI Values') + axes[1, 0].axis('off') + plt.colorbar(im, ax=axes[1, 0], label='NDWI') + + # Water likely + axes[1, 1].imshow(water_likely, cmap='Blues') + axes[1, 1].set_title(f'Water Likely\n{flag_summary["water_likely_pixels"]:,} pixels') + axes[1, 1].axis('off') + + # Water uncertain + axes[1, 2].imshow(water_uncertain, cmap='Oranges') + axes[1, 2].set_title(f'Water Uncertain\n{flag_summary["water_uncertain_pixels"]:,} pixels') + axes[1, 2].axis('off') + + plt.tight_layout() + fig.savefig(os.path.join(output_dir, f"{base_name}_quality_flags.png"), dpi=150) + plt.close() + + print(f"Saved quality visualization: {os.path.join(output_dir, f'{base_name}_quality_flags.png')}") + + +# ============================================================================= +# Main Integration Function +# ============================================================================= + +def process_with_quality_filtering(image_path, ndwi_results=None, slope_data=None, + output_dir=None): + """ + Complete quality flag processing pipeline. + + Args: + image_path (str): Path to satellite imagery + ndwi_results (dict, optional): Results from NDWI processing + slope_data (np.ndarray, optional): Slope array + output_dir (str, optional): Output directory + + Returns: + dict: Quality processing results + """ + if output_dir is None: + output_dir = "result_quality_filtering" + os.makedirs(output_dir, exist_ok=True) + + print("=" * 60) + print("QUALITY FLAG GENERATION AND FILTERING") + print("=" * 60) + + # Generate quality flags + quality, flag_summary = generate_quality_flags(image_path, slope_data=slope_data) + + # Get NDWI data + if ndwi_results is not None and 'ndwi' in ndwi_results: + ndwi_data = ndwi_results['ndwi'] + else: + from rastertools import calculate_ndwi + ndwi_data = calculate_ndwi(image_path, plot=False) + with rio.open(ndwi_data) as src: + ndwi_data = src.read(1) + + # Apply quality filtering + filtered_ndwi, reliable_mask = apply_quality_filter(ndwi_data, quality) + + # Save outputs + base_name = os.path.splitext(os.path.basename(image_path))[0] + + # Save quality flags + with rio.open(image_path) as src: + profile = src.profile.copy() + profile.update(dtype='uint8', count=1, nodata=255) + + quality_path = os.path.join(output_dir, f"{base_name}_quality_flags.tif") + with rio.open(quality_path, 'w', **profile) as dst: + dst.write(quality, 1) + print(f"Saved quality flags: {quality_path}") + + # Save reliability mask + reliable_path = os.path.join(output_dir, f"{base_name}_reliability_mask.tif") + profile.update(dtype='uint8', nodata=255) + with rio.open(reliable_path, 'w', **profile) as dst: + dst.write(reliable_mask.astype(np.uint8) * 255, 1) + print(f"Saved reliability mask: {reliable_path}") + + # Save visualization + save_quality_visualization(quality, flag_summary, ndwi_data, output_dir, base_name) + + return { + 'quality_flags': quality, + 'flag_summary': flag_summary, + 'filtered_ndwi': filtered_ndwi, + 'reliable_mask': reliable_mask + } + + +# ============================================================================= +# Main Execution +# ============================================================================= + +if __name__ == "__main__": + config = load_config() + + from load_config import get_image_path + image_path = get_image_path(config, 0) + + if image_path and os.path.exists(image_path): + results = process_with_quality_filtering(image_path) + print("\nQuality flag processing complete!") + else: + print("No image available for quality flag processing.")