FireVulnerability.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
FireVulnerability.py: Python implementation of FireVulnerability.js using Earth Engine Python API.

This script maps fire vulnerability based on long-term trends, static factors, 
and historical fire data using a Random Forest classifier.

Steps:
1. Load study area boundary.
2. Load and prepare historical fire points.
3. Load and process auxiliary datasets (NICFI, MODIS LST, DEM, Roads).
4. Calculate trends for NICFI and MODIS LST.
5. Load pre-computed trend layers from TrendFire.
6. Combine all input features into a single multi-band image.
7. Resample the combined feature image to a consistent grid (e.g., 30m EPSG:3857).
8. Prepare training data by sampling features at fire and non-fire locations.
9. Train a Random Forest classifier.
10. Evaluate the classifier.
11. Apply the classifier to the entire study area to create the vulnerability map.
12. Export the final vulnerability map.
"""

import ee
import os
import datetime
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import geopandas as gpd
from tqdm import tqdm
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, confusion_matrix, classification_report
import seaborn as sns

# Define the Goa study area boundary
def get_study_boundary():
    """
    Get the study area boundary for Goa, India.
    Returns an ee.Geometry.
    """
    # Method 1: Load from a shapefile if available locally
    try:
        boundary_path = os.path.join('data', 'pa_boundary.shp')
        if os.path.exists(boundary_path):
            boundary_gdf = gpd.read_file(boundary_path)
            # Convert to GeoJSON
            boundary_geojson = boundary_gdf.geometry.__geo_interface__
            # Create an EE geometry
            return ee.Geometry.Polygon(boundary_geojson['features'][0]['geometry']['coordinates'])
    except Exception as e:
        print(f"Could not load boundary from shapefile: {e}")
    
    # Method 2: Use predefined feature from Earth Engine
    try:
        # Get the Goa district boundary from Earth Engine's administrative boundaries
        goa = ee.FeatureCollection("FAO/GAUL/2015/level1").filter(ee.Filter.eq('ADM1_NAME', 'Goa'))
        return goa.geometry()
    except Exception as e:
        print(f"Could not load boundary from Earth Engine: {e}")
        
    # Method 3: Define manually as fallback
    # These coordinates would need to be adjusted to match the study area
    goa_coords = [
        [73.6765, 15.7560],
        [74.3161, 15.7560],
        [74.3161, 14.8922],
        [73.6765, 14.8922],
        [73.6765, 15.7560]
    ]
    return ee.Geometry.Polygon(goa_coords)

# Function to load fire events data
def load_fire_events(pa):
    """
    Load fire events data from 2013-2023 and merge into a single feature collection.
    
    Args:
        pa: ee.Geometry representing the study area
    
    Returns:
        ee.FeatureCollection of fire events
    """
    print("Loading fire events data...")
    
    # Load fire events data (2013-2019)
    try:
        # If the data is available in Earth Engine, you can load it directly
        fire_13_19 = ee.FeatureCollection("users/your_username/fire13_19")
    except:
        # Otherwise, try to find and upload it, or use a placeholder
        print("Fire data 2013-2019 not available in Earth Engine. Using placeholder.")
        fire_13_19 = ee.FeatureCollection([])
    
    # Load fire events data (2020-2023)
    try:
        # If the data is available in Earth Engine, you can load it directly
        fire_20_23 = ee.FeatureCollection("users/your_username/fire20_23")
    except:
        # Otherwise, try to find and upload it, or use a placeholder
        print("Fire data 2020-2023 not available in Earth Engine. Using placeholder.")
        fire_20_23 = ee.FeatureCollection([])
    
    # Merge the two collections
    fire_all = fire_13_19.merge(fire_20_23)
    
    # Filter by the study area boundary
    fire_all = fire_all.filterBounds(pa)
    
    # Add a 'year' property based on the acquisition date
    def add_year(feature):
        acq_date = ee.String(feature.get('ACQ_DATE'))
        year = ee.Number.parse(acq_date.slice(0, 4))
        return feature.set('year', year)
    
    fire_all = fire_all.map(add_year)
    
    # Filter to include only 2013-2022 for training (2023 will be used for validation)
    fire_training = fire_all.filter(ee.Filter.lt('year', 2023))
    
    # Count the number of fire points
    training_count = fire_training.size().getInfo()
    print(f"Number of fire events for training (2013-2022): {training_count}")
    
    # Get 2023 fire events for validation
    fire_validation = fire_all.filter(ee.Filter.eq('year', 2023))
    validation_count = fire_validation.size().getInfo()
    print(f"Number of fire events for validation (2023): {validation_count}")
    
    return fire_all, fire_training, fire_validation

# Function to load road data
def load_roads(pa):
    """
    Load roads data and clip to the study area.
    
    Args:
        pa: ee.Geometry representing the study area
    
    Returns:
        ee.Image with roads data
    """
    print("Loading roads data...")
    
    try:
        # Try to load the road data from Earth Engine assets
        roads = ee.Image("users/your_username/road").clip(pa)
        return roads
    except:
        print("Road data not available in Earth Engine. Using OSM roads as fallback.")
        # Use OpenStreetMap data as a fallback
        roads = ee.FeatureCollection('OSM/2020/roads').filterBounds(pa)
        
        # Convert to an image
        roads_image = roads.reduceToImage(
            properties=['highway'],
            reducer=ee.Reducer.first()
        ).clip(pa)
        
        return roads_image

# Function to generate random non-fire points
def generate_non_fire_points(pa, fire_training, num_points=100, buffer_distance=500):
    """
    Generate random points for non-fire locations within the study area.
    
    Args:
        pa: ee.Geometry representing the study area
        fire_training: ee.FeatureCollection of fire points
        num_points: Number of random points to generate
        buffer_distance: Buffer distance (in meters) to avoid fire locations
    
    Returns:
        ee.FeatureCollection of non-fire points
    """
    print(f"Generating {num_points} random non-fire points...")
    
    # Create a buffer around fire points to avoid
    fire_buffer = fire_training.geometry().buffer(buffer_distance)
    
    # Create a mask to avoid sampling in fire buffer zones
    mask = ee.Image(1).clip(pa).updateMask(
        ee.Image(1).clip(pa).mask().subtract(
            ee.Image(1).clip(fire_buffer).mask()
        )
    )
    
    # Generate random points in the valid sampling area
    non_fire_points = ee.FeatureCollection.randomPoints({
        'region': pa,
        'points': num_points,
        'seed': 42,
        'mask': mask
    })
    
    # Add a risk property (0 for non-fire)
    non_fire_points = non_fire_points.map(lambda f: f.set('RiskNumeric', 0))
    
    # Check if we have enough points
    point_count = non_fire_points.size().getInfo()
    print(f"Generated {point_count} non-fire points")
    
    return non_fire_points

# Function to calculate Land Surface Temperature trends from MODIS data
def calculate_lst_trends(pa, start_date='2013-01-01', end_date='2022-12-31'):
    """
    Calculate Land Surface Temperature (LST) trends from MODIS data.
    
    Args:
        pa: ee.Geometry representing the study area
        start_date: Start date for trend calculation
        end_date: End date for trend calculation
    
    Returns:
        ee.Image with LST trend bands
    """
    print("Calculating Land Surface Temperature trends...")
    
    # Function to mask clouds in MODIS LST imagery
    def maskLstClouds(image):
        qc = image.select('QC_Day')
        # Bit 0-1: LST produced, good quality
        mask = qc.bitwiseAnd(3).eq(0)
        return image.updateMask(mask)
    
    # Load the MODIS LST data
    modis_lst = ee.ImageCollection('MODIS/061/MOD11A1') \
        .filterDate(start_date, end_date) \
        .filterBounds(pa)
    
    # Apply cloud masking
    modis_lst_masked = modis_lst.map(maskLstClouds)
    
    # Convert LST to Celsius
    modis_lst_celsius = modis_lst_masked.map(
        lambda image: image.select(['LST_Day_1km', 'LST_Night_1km'])
                          .multiply(0.02)
                          .subtract(273.15)
                          .copyProperties(image, ['system:time_start'])
    )
    
    # Function to add time band for trend calculation
    def addTime(image):
        # Get the timestamp and convert to years since start
        year = ee.Date(image.get('system:time_start')).difference(
            ee.Date(start_date), 'year')
        return image.addBands(ee.Image(year).rename('t'))
    
    # Add time band to the collection
    modis_lst_with_time = modis_lst_celsius.map(addTime)
    
    # Calculate trends for day and night LST
    day_trend = modis_lst_with_time.select(['t', 'LST_Day_1km']) \
        .reduce(ee.Reducer.linearFit()) \
        .select(['scale'], ['lst_day_trend'])
    
    night_trend = modis_lst_with_time.select(['t', 'LST_Night_1km']) \
        .reduce(ee.Reducer.linearFit()) \
        .select(['scale'], ['lst_night_trend'])
    
    # Calculate delta LST (day-night)
    def calculateDeltaLST(image):
        day = image.select('LST_Day_1km')
        night = image.select('LST_Night_1km')
        delta = day.subtract(night).rename('delta_lst')
        return image.addBands(delta)
    
    modis_lst_with_delta = modis_lst_celsius.map(calculateDeltaLST)
    
    # Calculate delta LST trend
    delta_trend = modis_lst_with_delta.map(addTime) \
        .select(['t', 'delta_lst']) \
        .reduce(ee.Reducer.linearFit()) \
        .select(['scale'], ['delta_lst_trend'])
    
    # Merge all LST trend layers
    lst_trends = day_trend.addBands(night_trend).addBands(delta_trend)
    
    return lst_trends.clip(pa)

# Function to load trend layers from TrendFire output
def load_trend_layers(pa):
    """
    Load trend layers previously calculated using TrendFire.py.
    
    Args:
        pa: ee.Geometry representing the study area
    
    Returns:
        ee.Image with all trend bands
    """
    print("Loading trend layers...")
    
    try:
        # Try to load the trend layers from Earth Engine assets
        trends = ee.Image("users/your_username/goa_fire_trends").clip(pa)
        print("Successfully loaded trend layers from Earth Engine asset.")
        return trends
    except:
        print("Trend layers not found in Earth Engine assets. Calculating LST trends as fallback.")
        # Calculate LST trends as a fallback
        lst_trends = calculate_lst_trends(pa)
        return lst_trends

# Function to categorize fire risk based on temperature difference
def categorize_fire_risk(fire_training, lst_data, pa):
    """
    Categorize fire points by risk level based on temperature difference.
    
    Args:
        fire_training: ee.FeatureCollection of fire points
        lst_data: ee.Image with LST data
        pa: ee.Geometry representing the study area
    
    Returns:
        ee.FeatureCollection of categorized fire points
    """
    print("Categorizing fire risk levels...")
    
    # Sample the delta LST values at fire locations
    fire_samples = lst_data.select('delta_lst_trend').sampleRegions({
        'collection': fire_training,
        'properties': ['ACQ_DATE', 'year'],
        'scale': 1000,  # MODIS LST resolution
        'geometries': True
    })
    
    # Define threshold values for risk categories
    # These thresholds would need to be calibrated based on local conditions
    def categorize(feature):
        delta_lst = feature.get('delta_lst_trend')
        
        # Define risk categories based on delta LST
        # 1: Low risk, 2: Moderate risk, 3: High risk
        risk = ee.Algorithms.If(ee.Number(delta_lst).lt(0.1), 1,
               ee.Algorithms.If(ee.Number(delta_lst).lt(0.2), 2, 3))
        
        return feature.set('RiskNumeric', risk)
    
    # Apply categorization to all fire points
    categorized_fire = fire_samples.map(categorize)
    
    # Count points in each risk category
    low_risk = categorized_fire.filter(ee.Filter.eq('RiskNumeric', 1)).size().getInfo()
    moderate_risk = categorized_fire.filter(ee.Filter.eq('RiskNumeric', 2)).size().getInfo()
    high_risk = categorized_fire.filter(ee.Filter.eq('RiskNumeric', 3)).size().getInfo()
    
    print(f"Fire risk categorization: {low_risk} Low, {moderate_risk} Moderate, {high_risk} High risk")
    
    return categorized_fire

# Function to prepare training data
def prepare_training_data(fire_points_input, boundary, predictor_image, label_property, \
                            num_random_points=90, random_points_seed=42, \
                            moderate_risk_limit=70, use_risk_categories=True):
    """
    Prepares training and validation data by combining categorized fire points 
    with generated non-fire points and sampling predictor values.

    Args:
        fire_points_input: ee.FeatureCollection of historical fire points. 
                           Should contain 'Delta T' if use_risk_categories=True.
        boundary: ee.Geometry defining the study area.
        predictor_image: ee.Image containing all predictor bands.
        label_property: String, the name of the property to use for classification label (e.g., 'RiskNumeric').
        num_random_points: Integer, number of random non-fire points to generate.
        random_points_seed: Integer, seed for random point generation.
        moderate_risk_limit: Integer, threshold for Delta T to limit moderate risk points (if used).
        use_risk_categories: Boolean, if True, categorize fires (0=NoFire, 1=Low, 2=Moderate, 3=High).
                             If False, label all fires as 1 (binary fire/no-fire).

    Returns:
        An ee.FeatureCollection containing sampled points with predictor values and the class label,
        or None on error.
    """
    print("\nPreparing training data...")
    try:
        # --- Categorize Fire Points ---
        print("Categorizing fire points...")
        no_fire_label = 0
        fire_label = 1 # Used if use_risk_categories is False

        if use_risk_categories:
            # Use tqdm to show progress for client-side mapping/filtering if complex
            # Note: The actual .map/.filter happens server-side, so tqdm only tracks
            # the *creation* of these instructions, not their execution time.
            print("  Defining fire risk categories (Low, Moderate, High)...")
            low_risk = fire_points_input.filter(ee.Filter.rangeContains('Delta T', 0, 25))\
                                       .map(lambda f: f.set(label_property, 1).set('RiskCategory', 'Low Risk'))
            moderate_risk = fire_points_input.filter(ee.Filter.rangeContains('Delta T', 25, 35))\
                                             .map(lambda f: f.set(label_property, 2).set('RiskCategory', 'Moderate Risk'))\
                                             .limit(moderate_risk_limit) # Limit moderate risk points like JS
            high_risk = fire_points_input.filter(ee.Filter.greaterThanOrEquals('Delta T', 35))\
                                         .map(lambda f: f.set(label_property, 3).set('RiskCategory', 'High Risk'))
            
            # Get counts (these trigger computation)
            print("  Fetching category counts...") 
            low_count = low_risk.size().getInfo()
            mod_count = moderate_risk.size().getInfo()
            high_count = high_risk.size().getInfo()
            print(f"  Low Risk (<25): {low_count}")
            print(f"  Moderate Risk (25-35, limited to {moderate_risk_limit}): {mod_count}")
            print(f"  High Risk (>=35): {high_count}")
            
            # Combine categorized fire points
            categorized_fires = low_risk.merge(moderate_risk).merge(high_risk)
            total_fires = categorized_fires.size().getInfo()
            print(f"  Total categorized fire points used: {total_fires}")

        else:
            # Binary classification: Label all input fires as 1
            categorized_fires = fire_points_input.map(lambda f: f.set(label_property, fire_label))
            total_fires = categorized_fires.size().getInfo()
            print(f"  Using binary classification. Total fire points: {total_fires}")

        # --- Generate Non-Fire Points ---
        print(f"Generating {num_random_points} random non-fire points...")
        random_pts = ee.FeatureCollection.randomPoints(
            region=boundary, 
            points=num_random_points, 
            seed=random_points_seed
        )
        non_fire_pts = random_pts.map(lambda feat: feat.set(label_property, no_fire_label))
        non_fire_count = non_fire_pts.size().getInfo()
        print(f"  Generated {non_fire_count} non-fire points.")

        # --- Combine Fire and Non-Fire Points ---
        training_validation_points = categorized_fires.merge(non_fire_pts)
        total_points = training_validation_points.size().getInfo()
        print(f"Total points for sampling: {total_points}")

        # --- Sample Predictor Values ---
        print("\n>>> Sending request to GEE: Sample predictor values at points...") # Log Before
        if predictor_image is None:
             print("Error: Predictor image is missing.")
             return None
        if training_validation_points is None or total_points == 0:
            print("Error: No training/validation points available for sampling.")
            return None

        # Ensure predictor image bands are float for sampling/training
        predictor_image = predictor_image.toFloat() 
        
        # Define properties to keep (only the label)
        properties_to_keep = [label_property]

        # Sample the regions
        sampled_data = predictor_image.sampleRegions(
            collection=training_validation_points,
            properties=properties_to_keep, 
            scale=predictor_image.projection().nominalScale(), # Use image native scale
            geometries=False # Don't need geometries for training
        )

        # Force execution and check size BEFORE adding random column
        print(">>> Waiting for GEE: Sampling predictor values...")
        sampled_size_check = sampled_data.size().getInfo()
        print(f"<<< GEE Response: Sampling returned {sampled_size_check} features.") # Log After

        if sampled_size_check == 0:
            print("Error: Sampling returned no features. Check predictor image, points, and permissions.")
            return None

        # Add a random column for splitting later
        print("Adding random column for train/test split...")
        sampled_data = sampled_data.randomColumn(columnName='random', seed=random_points_seed) # Explicit column name and seed
        # Verify random column exists (optional check)
        # print("Properties after adding random column:", sampled_data.first().propertyNames().getInfo())
            
        print("Training data preparation complete.")
        return sampled_data

    except ee.EEException as e: 
        print(f"GEE error during training data preparation: {e}")
        return None
    except Exception as e: 
        print(f"Unexpected error during training data preparation: {e}")
        return None

# Function to train a Random Forest classifier
def train_classifier(training_data, class_property, input_properties, num_trees=1000, variables_per_split=1, seed=42):
    """
    Trains the ee.Classifier.smileRandomForest classifier.

    Args:
        training_data: ee.FeatureCollection used for training.
        class_property: String, the name of the property containing the class label.
        input_properties: List of strings, names of the predictor properties/bands.
        num_trees: Integer, number of trees in the forest.
        variables_per_split: Integer, number of variables to consider at each split.
                             Set to None or 0 to use sqrt(num_variables).
                             JS used 1.
        seed: Integer, seed for reproducibility.

    Returns:
        A trained ee.Classifier object, or None on error.
    """
    print(f"Training Random Forest classifier with {num_trees} trees...")
    
    if training_data is None:
        print("Error: Training data is missing.")
        return None
    if not input_properties:
        print("Error: Input properties list is empty.")
        return None
        
    try:
        classifier = ee.Classifier.smileRandomForest(
            numberOfTrees=num_trees, 
            variablesPerSplit=variables_per_split, # Match JS, using 1
            # minLeafPopulation=1, # Default is 1
            # bagFraction=0.5, # Default is 0.5
            seed=seed
        ).train(
            features=training_data, 
            classProperty=class_property, 
            inputProperties=input_properties
        )
        
        print("Classifier trained successfully.")
        
        # Optional: Explain the classifier (can be computationally intensive)
        # try:
        #     explanation = classifier.explain()
        #     print("Classifier explanation:", explanation.getInfo())
        #     importance = explanation.get('importance')
        #     if importance:
        #          print("Variable Importance:")
        #          importance_info = importance.getInfo()
        #          # Ensure importance_info is a dictionary before iterating
        #          if isinstance(importance_info, dict):
        #              for prop, imp in importance_info.items():
        #                  print(f"  {prop}: {imp}")
        #          else:
        #              print("  Could not retrieve importance details in expected format.")
        # except Exception as explain_e:
        #      print(f"Could not get classifier explanation: {explain_e}")

        return classifier
        
    except ee.EEException as e:
        print(f"GEE error during classifier training: {e}")
        return None
    except Exception as e:
        print(f"Unexpected error during classifier training: {e}")
        return None

# Function to validate the Random Forest model
def validate_model(classifier, predictor_names, validation_data):
    """
    Validate the Random Forest model using test data.
    
    Args:
        classifier: Trained ee.Classifier.RandomForest
        predictor_names: List of predictor variable names
        validation_data: ee.FeatureCollection with validation data
    
    Returns:
        Dictionary with validation metrics
    """
    print("Validating Random Forest model...")
    
    # Apply the classifier to the validation data
    validated = validation_data.classify(classifier)
    
    # Get the confusion matrix
    confusion_matrix = validated.errorMatrix('RiskNumeric', 'classification')
    
    # Calculate accuracy and Kappa statistics
    accuracy = confusion_matrix.accuracy().getInfo()
    kappa = confusion_matrix.kappa().getInfo()
    
    print(f"Validation accuracy: {accuracy:.4f}")
    print(f"Kappa statistic: {kappa:.4f}")
    
    # Get confusion matrix as a list of lists
    cm_data = confusion_matrix.getInfo()
    
    # Return validation metrics
    return {
        'accuracy': accuracy,
        'kappa': kappa,
        'confusion_matrix': cm_data
    }

# Function to create a fire vulnerability map
def create_vulnerability_map(classifier, predictor_variables, predictor_names, pa):
    """
    Create a fire vulnerability map using the trained classifier.
    
    Args:
        classifier: Trained ee.Classifier.RandomForest
        predictor_variables: ee.Image with predictor variables
        predictor_names: List of predictor variable names
        pa: ee.Geometry representing the study area
    
    Returns:
        ee.Image with fire vulnerability classification
    """
    print("Creating fire vulnerability map...")
    
    # Select only the predictor bands used for training
    predictors = predictor_variables.select(predictor_names)
    
    # Apply the classifier to create a vulnerability map
    vulnerability_map = predictors.classify(classifier).clip(pa)
    
    return vulnerability_map

# Function to export the vulnerability map to Google Drive
def export_vulnerability_map(vulnerability_map, pa, filename='fire_vulnerability_map', folder='ForestFireGoa'):
    """
    Export the fire vulnerability map to Google Drive.
    
    Args:
        vulnerability_map: ee.Image with fire vulnerability classification
        pa: ee.Geometry representing the study area
        filename: Output filename
        folder: Google Drive folder name
    """
    print(f"Exporting fire vulnerability map to Google Drive folder '{folder}'...")
    
    # Start the export task
    task = ee.batch.Export.image.toDrive({
        'image': vulnerability_map,
        'description': filename,
        'folder': folder,
        'fileNamePrefix': filename,
        'region': pa,
        'scale': 30,
        'maxPixels': 1e13
    })
    
    task.start()
    print(f"Export task started with ID: {task.id}")

# Function to export the vulnerability map to an Earth Engine asset
def export_vulnerability_map_to_asset(vulnerability_map, pa, asset_name):
    """
    Export the fire vulnerability map to an Earth Engine asset.
    
    Args:
        vulnerability_map: ee.Image with fire vulnerability classification
        pa: ee.Geometry representing the study area
        asset_name: Name of the Earth Engine asset
    """
    print(f"Exporting fire vulnerability map to Earth Engine asset: {asset_name}")
    
    # Start the export task
    task = ee.batch.Export.image.toAsset({
        'image': vulnerability_map,
        'description': asset_name.split('/')[-1],
        'assetId': asset_name,
        'region': pa,
        'scale': 30,
        'maxPixels': 1e13
    })
    
    task.start()
    print(f"Export task started with ID: {task.id}")

# Main function to run the fire vulnerability mapping workflow
def main():
    """Main script execution for fire vulnerability mapping."""
    # --- Configuration ---
    PROJECT_ID = 'ee-crop-health-telangana'
    BOUNDARY_ASSET_ID = 'users/jonasnothnagel/pa_boundary'
    FIRE_ASSET_BASE = 'users/jonasnothnagel/' # Base for fire points (e.g., fire13_19, fire20_23)
    DEM_ASSET_ID = 'users/jonasnothnagel/dem' 
    ROAD_ASSET_ID = 'users/jonasnothnagel/roads' 
    TRENDFIRE_ASSET_BASE = 'users/jonasnothnagel/' # Base for JS TrendFire outputs (e.g., Trend2023_ndvi)
    NICFI_COLLECTION = 'projects/planet-nicfi/assets/basemaps/asia' 
    MODIS_LST_COLLECTION = 'MODIS/061/MOD11A1'
    
    # --- Option to use pre-computed resampled input --- #
    USE_PRECOMPUTED_RESAMPLED = True # Set to True to load the asset below
    PRECOMPUTED_RESAMPLED_ASSET = 'users/jonasnothnagel/inputResampled' # Asset ID of the uploaded 40-band tif
    # ---------------------------------------------------- #
    
    OUTPUT_ASSET_BASE = 'users/jonasnothnagel/'
    VULNERABILITY_MAP_ASSET = OUTPUT_ASSET_BASE + 'FireVulnerability_py'
    # Asset ID for the resampled image *generated* by this script (if USE_PRECOMPUTED_RESAMPLED=False)
    GENERATED_RESAMPLED_FEATURES_ASSET = OUTPUT_ASSET_BASE + 'FireVulnerability_InputsGenerated_py' 
    TRAINING_POINTS_ASSET = OUTPUT_ASSET_BASE + 'FireVulnerability_TrainingPoints_py'
    
    GRID_SCALE = 30
    GRID_CRS = 'EPSG:3857'
    NUM_RANDOM_POINTS = 90 # Match JS
    RF_TREES = 1000 # Match JS
    RF_VAR_SPLIT = 1 # Match JS
    TRAIN_TEST_SPLIT_RATIO = 0.8
    CLASS_LABEL = 'RiskNumeric'
    # Use Delta T categories (True) or binary fire/no-fire (False)
    # JS used only high risk (Delta T >= 35) vs non-fire (0). Setting this to True enables categorization,
    # but prepare_training_data currently merges all categories or uses binary.
    # To exactly match JS sampling, logic inside prepare_training_data needs adjustment.
    USE_RISK_CATEGORIES = True 
    RANDOM_SEED = 42 # For reproducibility in random points and RF training

    # --- Initialization ---
    initialize_ee(PROJECT_ID)
    boundary = get_boundary(BOUNDARY_ASSET_ID)
    if not boundary:
        return

    # --- Processing Steps ---
    print("\n--- Starting Fire Vulnerability Workflow ---")
    
    input_resampled = None # This will hold the final predictor image for sampling/classification
    fire_pts_processed = None # To hold loaded fire points
    
    # --- Mode 1: Load Precomputed Resampled Input --- #
    if USE_PRECOMPUTED_RESAMPLED:
        print(f"---> Mode: Using Precomputed Input <---")
        print(f"Attempting to load pre-computed resampled input: {PRECOMPUTED_RESAMPLED_ASSET}")
        try:
            input_resampled = ee.Image(PRECOMPUTED_RESAMPLED_ASSET).toFloat() # Ensure float
            # Verify it has bands
            if not input_resampled.bandNames().getInfo():
                 raise ee.EEException(f"Loaded asset {PRECOMPUTED_RESAMPLED_ASSET} has no bands.")
            print("Successfully loaded pre-computed resampled input.")
            print("Skipping intermediate data loading, combining, and resampling steps.")
            # Load fire points separately as they are needed for training data prep
            fire_pts_processed = load_fire_points(FIRE_ASSET_BASE, boundary)
            if fire_pts_processed is None:
                 print("Error: Failed to load fire points even when using precomputed input. Exiting.")
                 return
        except ee.EEException as e:
            print(f"Error loading pre-computed asset '{PRECOMPUTED_RESAMPLED_ASSET}': {e}")
            print("Exiting. Cannot proceed without input features.")
            return # Exit if loading precomputed fails
        except Exception as e:
             print(f"Unexpected error loading pre-computed asset: {e}")
             print("Exiting. Cannot proceed without input features.")
             return # Exit if loading precomputed fails
             
    # --- Mode 2: Generate Inputs from Components --- #
    if input_resampled is None: # This block runs if USE_PRECOMPUTED_RESAMPLED is False or loading failed
        print("---> Mode: Generating Input Features from Components <--- ")
        # 1. Load Fire Points (needed regardless for training data)
        fire_pts_processed = load_fire_points(FIRE_ASSET_BASE, boundary)
        if fire_pts_processed is None:
            print("Exiting: Failed to load fire points.")
            return
            
        # 2. Process NICFI (Optional)
        nicfi_trends = process_nicfi(boundary)
        if nicfi_trends is None:
             print("Warning: Failed to process NICFI trends. Proceeding without them.")
        
        # 3. Process MODIS LST (Optional)
        modis_trends = process_modis_lst(boundary)
        if modis_trends is None:
             print("Warning: Failed to process MODIS LST trends. Proceeding without them.")
             
        # 4. Load DEM Slope
        dem_slope = load_dem_slope(DEM_ASSET_ID, boundary)
        if dem_slope is None:
             print("Warning: Failed to process DEM slope. Proceeding without it.")
             
        # 5. Load Roads
        roads_img = load_roads(ROAD_ASSET_ID, boundary) # Corrected arguments
        if roads_img is None:
             print("Warning: Failed to process Roads layer. Proceeding without it.")
             
        # 6. Load TrendFire Outputs (Core Trends)
        trendfire_img = load_trendfire_outputs(TRENDFIRE_ASSET_BASE, boundary)
        if trendfire_img is None:
            print("Exiting: Failed to load TrendFire outputs.")
            return
            
        # 7. Combine Available Features
        print("Combining available features...")
        combined_features = combine_features(trendfire_img, nicfi_trends, modis_trends, dem_slope, roads_img)
        if combined_features is None or not combined_features.bandNames().getInfo():
             print("Exiting: Failed to combine input features or resulting image has no bands.")
             return
        print(f"Combined features image created with bands: {combined_features.bandNames().getInfo()}")
             
        # 8. Resample & Export Generated Intermediate Input
        print(f"Resampling combined features to {GRID_SCALE}m scale, CRS {GRID_CRS}...")
        input_resampled = resample_features(combined_features, GRID_SCALE, GRID_CRS, boundary)
        if input_resampled is None or not input_resampled.bandNames().getInfo():
            print("Exiting: Failed to resample features or resulting image has no bands.")
            return
        print("Resampling complete.")
            
        # Export the generated resampled image for comparison/debugging
        print(f"Exporting the generated resampled input features to asset: {GENERATED_RESAMPLED_FEATURES_ASSET}")
        export_asset(input_resampled, GENERATED_RESAMPLED_FEATURES_ASSET, boundary, "Generated_Inputs_FV_py", GRID_SCALE, GRID_CRS)
    
    # --- Classification Workflow (Uses the selected 'input_resampled') --- #
    if input_resampled is None or fire_pts_processed is None:
         print("Error: Cannot proceed to classification. Missing resampled input image or fire points.")
         return
         
    print(f"\nProceeding with classification using predictor image with bands: {input_resampled.bandNames().getInfo()}")
             
    # 9. Prepare Training Data
    training_data_fc = prepare_training_data(
        fire_points_input=fire_pts_processed, 
        boundary=boundary, 
        predictor_image=input_resampled, 
        label_property=CLASS_LABEL, 
        num_random_points=NUM_RANDOM_POINTS, 
        random_points_seed=RANDOM_SEED,
        use_risk_categories=USE_RISK_CATEGORIES
    )
    if training_data_fc is None:
        print("Exiting: Failed to prepare training data.")
        return
        
    # Export the full sampled dataset before splitting (useful for external analysis)
    export_table(training_data_fc, TRAINING_POINTS_ASSET, "Training_Validation_Points_FV_py")
    
    # 10. Split Data
    print(f"Splitting data using ratio: {TRAIN_TEST_SPLIT_RATIO:.1f} train / {1-TRAIN_TEST_SPLIT_RATIO:.1f} validation...")
    # Filter based on the random column added in prepare_training_data
    training_sample = training_data_fc.filter(ee.Filter.lt('random', TRAIN_TEST_SPLIT_RATIO))
    validation_sample = training_data_fc.filter(ee.Filter.gte('random', TRAIN_TEST_SPLIT_RATIO))
    train_count = training_sample.size().getInfo()
    val_count = validation_sample.size().getInfo()
    print(f"  Training samples: {train_count}, Validation samples: {val_count}")
    if train_count == 0 or val_count == 0:
        print("Error: Training or validation sample is empty after splitting. Check data preparation.")
        return
    
    # Get predictor names from the resampled image (used as input properties for RF)
    predictor_names = input_resampled.bandNames()
    if not predictor_names.getInfo():
         print("Error: Could not retrieve band names from the input predictor image.")
         return
         
    # 11. Train Classifier
    print("\n>>> Sending request to GEE: Train Random Forest classifier...") # Log Before
    classifier = train_classifier(
        training_data=training_sample, 
        class_property=CLASS_LABEL, 
        input_properties=predictor_names, # Pass Python list
        num_trees=RF_TREES, 
        variables_per_split=RF_VAR_SPLIT, 
        seed=RANDOM_SEED
    )
    # Note: train() itself can take time, but the result (trained classifier object)
    # is usually returned quickly to the client unless explain() is called.
    # The heavy lifting happens when classifier.classify() is called.
    print("<<< GEE Response: Classifier training definition complete.") # Log After

    if classifier is None:
        print("Exiting: Classifier training failed.")
        return
        
    # 12. Evaluate Classifier
    print("\n>>> Sending request to GEE: Evaluate classifier...") # Log Before
    evaluation_results = evaluate_classifier(classifier, validation_sample, CLASS_LABEL)
    print("<<< GEE Response: Evaluation complete.") # Log After
    
    # evaluation_results is a dictionary containing accuracy, kappa, matrix etc.
    if evaluation_results is None:
        print("Warning: Classifier evaluation step failed or produced no results.")
    
    # 13. Classify Image
    print("\n>>> Sending request to GEE: Classify the full predictor image...") # Log Before
    # Ensure the image being classified has float type bands (required by some classifiers)
    vulnerability_map = input_resampled.toFloat().classify(classifier)
    # Force execution by getting some info (optional, but confirms the step is done)
    # print(">>> Waiting for GEE: Classification...")
    # _ = vulnerability_map.bandNames().getInfo() # Or use another light .getInfo() call
    print("<<< GEE Response: Classification definition complete.") # Log After
    
    # 14. Export Final Map
    print(f"\nExporting final Fire Vulnerability Map to asset:")
    print(f"  {VULNERABILITY_MAP_ASSET}")
    export_asset(
        image=vulnerability_map.toInt8(), # Classifiers output integer labels, cast to save space
        asset_id=VULNERABILITY_MAP_ASSET, 
        region=boundary, 
        description="Fire_Vulnerability_Map_py", 
        scale=GRID_SCALE, 
        crs=GRID_CRS
    )
    
    print("\n--- Fire Vulnerability Workflow Completed --- ")
    print("Monitor GEE tasks dashboard for export completion.")

# --- Fire Vulnerability Specific Functions ---

def load_fire_points(asset_base, boundary):
    """
    Loads historical fire point assets, merges them, and filters out 2023 data.
    
    Args:
        asset_base: String, the base path for fire point assets (e.g., 'users/username/').
        boundary: ee.Geometry object for optional filtering (though JS didn't explicitly filter here).
        
    Returns:
        An ee.FeatureCollection containing fire points from 2013-2022, or None on error.
    """
    print("Loading and preparing historical fire points...")
    fire_asset_13_19 = asset_base + 'fire13_19'
    fire_asset_20_23 = asset_base + 'fire20_23'
    
    try:
        print(f"  Loading {fire_asset_13_19}...")
        fire13_19 = ee.FeatureCollection(fire_asset_13_19)
        print(f"  Loading {fire_asset_20_23}...")
        fire20_23 = ee.FeatureCollection(fire_asset_20_23)
        
        # Merge the collections
        fire13_23 = fire13_19.merge(fire20_23)
        print(f"  Merged collections. Total points (2013-2023): {fire13_23.size().getInfo()}")
        
        # Filter out points from 2023 based on the 'acq_date' property
        # Assumes 'acq_date' format allows simple string filtering like 'YYYY-MM-DD'
        fire13_22 = fire13_23.filter(ee.Filter.stringContains('acq_date', '-2023').Not())
        # Alternative filter if acq_date is ee.Date object:
        # fire13_22 = fire13_23.filter(ee.Filter.date('2013-01-01', '2022-12-31')) 
        
        print(f"  Filtered out 2023 points. Total points (2013-2022): {fire13_22.size().getInfo()}")
        
        # Optional: Filter by boundary if needed, although JS didn't seem to do it here
        # fire13_22 = fire13_22.filterBounds(boundary)
        # print(f"  Filtered by boundary. Points within boundary: {fire13_22.size().getInfo()}")
        
        return fire13_22
        
    except ee.EEException as e:
        print(f"Error loading or processing fire point assets: {e}")
        print(f"Please ensure assets '{fire_asset_13_19}' and '{fire_asset_20_23}' exist and properties like 'acq_date' are correct.")
        return None
    except Exception as e:
        print(f"An unexpected error occurred loading fire points: {e}")
        return None

def process_nicfi(boundary, start_date='2016-03-01', end_date='2024-02-28', collection_id='projects/planet-nicfi/assets/basemaps/asia'):
    """
    Loads Planet-NICFI data, calculates NDVI, and computes trends for R, G, B, N, NDVI bands.
    
    Args:
        boundary: ee.Geometry for filtering and clipping.
        start_date: Start date for filtering NICFI data.
        end_date: End date for filtering NICFI data.
        collection_id: GEE asset ID for the NICFI collection.
        
    Returns:
        ee.Image containing slope and intercept bands for R, G, B, N, NDVI, or None on error.
    """
    print(f"Processing NICFI trends ({start_date} to {end_date})...")
    
    try:
        # Load the collection
        nicfi_col = ee.ImageCollection(collection_id) \
            .map(lambda img: img.clip(boundary)) \
            .filterDate(start_date, end_date)
            
        # Function to calculate and add NDVI band
        def addNDVI(image):
            # Ensure bands are correctly named ('N', 'R', 'G', 'B')
            # May need to rename depending on the specific NICFI collection version
            # Example: .select(['B4', 'B3', 'B2', 'B1'], ['N', 'R', 'G', 'B']) if needed
            ndvi = image.normalizedDifference(['N', 'R']).rename('NDVI')
            return image.addBands(ndvi)

        nicfi_col = nicfi_col.map(addNDVI)

        # Add time band (years since 2016-01-01, as per JS)
        ref_date = ee.Date('2016-01-01')
        def addTime(image):
            img_date = image.date()
            years = img_date.difference(ref_date, 'year')
            return image.addBands(ee.Image(years).rename('time').float())

        nicfi_col_time = nicfi_col.map(addTime)

        # Indices for trend calculation
        indices = ['R', 'G', 'B', 'N', 'NDVI']
        combined_trends = ee.Image().toFloat() # Start empty

        # Calculate trends for each index
        for index in indices:
            print(f"  Calculating NICFI trend for: {index}")
            stacked = nicfi_col_time.select(['time', index])
            regression = stacked.reduce(ee.Reducer.linearFit()) # Outputs 'scale' and 'offset'
            slope = regression.select('scale').rename(f'{index}_Slope')
            intercept = regression.select('offset').rename(f'{index}_Intercept')
            combined_trends = combined_trends.addBands(slope).addBands(intercept)
        
        print("NICFI trends calculated successfully.")
        return combined_trends

    except ee.EEException as e:
        print(f"GEE error during NICFI processing: {e}")
        print(f"Please ensure you have access to the collection: {collection_id}")
        return None
    except Exception as e:
        print(f"Unexpected error during NICFI processing: {e}")
        return None

def process_modis_lst(boundary, start_date='2013-01-01', end_date='2023-02-28', collection_id='MODIS/061/MOD11A1'):
    """
    Loads MODIS LST data, masks clouds, converts to Celsius, and computes trends for Day and Night LST.
    
    Args:
        boundary: ee.Geometry for filtering and clipping.
        start_date: Start date for filtering MODIS data.
        end_date: End date for filtering MODIS data.
        collection_id: GEE asset ID for the MODIS LST collection.
        
    Returns:
        ee.Image containing slope and intercept bands for LST_Day and LST_Night, or None on error.
    """
    print(f"Processing MODIS LST trends ({start_date} to {end_date})...")

    try:
        # Function to mask clouds/poor quality pixels and convert units
        def scale_and_mask_lst(image):
            qc = image.select('QC_Day')
            # Keep only pixels with good quality (Bits 0-1: 00 = Good)
            mask = qc.bitwiseAnd(3).eq(0) 
            
            lst_day_c = image.select('LST_Day_1km') \
                             .updateMask(mask) \
                             .multiply(0.02) \
                             .subtract(273.15) \
                             .rename('LST_Day')
                             
            lst_night_c = image.select('LST_Night_1km') \
                               .updateMask(mask) \
                               .multiply(0.02) \
                               .subtract(273.15) \
                               .rename('LST_Night')
                               
            return image.addBands(lst_day_c).addBands(lst_night_c).select(['LST_Day', 'LST_Night'])

        # Load, filter, and process the collection
        modis_lst_col = ee.ImageCollection(collection_id) \
            .filterDate(start_date, end_date) \
            .filterBounds(boundary) \
            .map(scale_and_mask_lst)

        # Add time band
        ref_date = ee.Date(start_date)
        def addTime(image):
            img_date = image.date()
            years = img_date.difference(ref_date, 'year')
            return image.addBands(ee.Image(years).rename('time').float())

        modis_col_time = modis_lst_col.map(addTime)

        # Calculate trends for Day and Night LST
        trends_list = []
        for band in ['LST_Day', 'LST_Night']:
             print(f"  Calculating MODIS trend for: {band}")
             stacked = modis_col_time.select(['time', band])
             regression = stacked.reduce(ee.Reducer.linearFit()) # scale, offset
             slope = regression.select('scale').rename(f'{band}_Slope')
             intercept = regression.select('offset').rename(f'{band}_Intercept')
             trends_list.append(slope.addBands(intercept))

        # Combine Day and Night trends
        if not trends_list:
            print("Error: Could not calculate any MODIS LST trends.")
            return None
            
        combined_trends = ee.ImageCollection(trends_list).toBands()
        # Rename bands from 0_LST_Day_Slope etc. to LST_Day_Slope etc.
        original_band_names = [img.bandNames() for img in trends_list]
        flat_band_names = [item for sublist in original_band_names for item in sublist.getInfo()]
        combined_trends = combined_trends.rename(flat_band_names)

        print("MODIS LST trends calculated successfully.")
        return combined_trends.clip(boundary) # Clip final result

    except ee.EEException as e:
        print(f"GEE error during MODIS LST processing: {e}")
        return None
    except Exception as e:
        print(f"Unexpected error during MODIS LST processing: {e}")
        return None

def load_dem_slope(asset_id, boundary):
    """Loads DEM asset, calculates slope, and clips to boundary."""
    print(f"Loading DEM ({asset_id}) and calculating slope...")
    try:
        dem = ee.Image(asset_id)
        slope = ee.Terrain.slope(dem)
        print("Slope calculated successfully.")
        return slope.clip(boundary) # Return only the slope band, clipped
    except ee.EEException as e:
        print(f"GEE error loading DEM or calculating slope: {e}")
        return None
    except Exception as e:
        print(f"Unexpected error during DEM/Slope processing: {e}")
        return None

def load_roads(asset_id, boundary):
    """Loads the road layer asset and clips to boundary."""
    print(f"Loading roads layer ({asset_id})...")
    try:
        roads_img = ee.Image(asset_id)
        print("Roads layer loaded successfully.")
        return roads_img.clip(boundary)
    except ee.EEException as e:
        print(f"GEE error loading roads asset: {e}")
        return None
    except Exception as e:
        print(f"Unexpected error during roads processing: {e}")
        return None

def load_trendfire_outputs(asset_base, boundary):
    """
    Loads the individual trend assets generated by TrendFire.js and merges them.
    Assumes assets are under the user's GEE path (asset_base).

    Args:
        asset_base: String, the base path for trend assets (e.g., 'users/username/').
        boundary: ee.Geometry object for clipping.

    Returns:
        ee.Image containing all merged trend bands (Slope and Intercept), or None on error.
    """
    print(f"Loading TrendFire outputs from base path: {asset_base}...")

    # Define asset names/parts based on TrendFire.js export structure
    landsat_indices = ['ndvi', 'evi', 'mirbi', 'ndfi', 'bsi', 'ndmi', 'nbr', 'nbr2', 'msavi', 'smi', 'ST_B10']
    # JS exported these as 'users/sravanthi_earthAnalytics/Trend2023_' + indexName
    # Assuming user ran JS and exported to their own folder:
    landsat_asset_prefix = 'Trend2023_' 
    rain_asset_id_part = 'Trend2022_rain_new'
    sm_asset_id_part = 'Trend2023_SM_new'
    # JS asset was 'Trend2024_RH_new'
    rh_asset_id_part = 'Trend2024_RH_new' 

    all_trend_bands = ee.Image().toFloat() # Start with an empty image
    asset_id = '' # For error message context

    try:
        # Load Landsat index trends
        for index in landsat_indices:
            asset_id = f"{asset_base}{landsat_asset_prefix}{index}"
            print(f"  Loading Landsat trend: {asset_id}")
            trend_img = ee.Image(asset_id)
            all_trend_bands = all_trend_bands.addBands(trend_img)

        # Load Rain trend
        asset_id = f"{asset_base}{rain_asset_id_part}"
        print(f"  Loading Rain trend: {asset_id}")
        rain_trend = ee.Image(asset_id)
        all_trend_bands = all_trend_bands.addBands(rain_trend)

        # Load SM trend
        asset_id = f"{asset_base}{sm_asset_id_part}"
        print(f"  Loading SM trend: {asset_id}")
        sm_trend = ee.Image(asset_id)
        # Check expected bands: sm_surface_Slope, sm_surface_Intercept
        # Rename if necessary, e.g., if JS exported with different names
        # sm_trend = sm_trend.rename(['sm_surface_Slope', 'sm_surface_Intercept'])
        all_trend_bands = all_trend_bands.addBands(sm_trend)

        # Load RH trend
        asset_id = f"{asset_base}{rh_asset_id_part}"
        print(f"  Loading RH trend: {asset_id}")
        rh_trend = ee.Image(asset_id)
        all_trend_bands = all_trend_bands.addBands(rh_trend)

        print("Successfully loaded and merged TrendFire outputs.")
        # Clip the final merged image to the boundary
        return all_trend_bands.clip(boundary)

    except ee.EEException as e:
        print(f"Error loading GEE asset during TrendFire output loading: {e}")
        print(f"Please ensure asset '{asset_id}' exists and you have permissions.")
        return None
    except Exception as e:
        print(f"An unexpected error occurred loading TrendFire outputs: {e}")
        return None

def combine_features(trendfire_img, nicfi_trends, modis_trends, dem_slope, roads_img):
    """
    Combines all input feature bands into a single ee.Image.
    Handles potential None inputs from previous steps.

    Args:
        trendfire_img: ee.Image with TrendFire trend bands.
        nicfi_trends: ee.Image with NICFI trend bands.
        modis_trends: ee.Image with MODIS LST trend bands.
        dem_slope: ee.Image with the DEM slope band.
        roads_img: ee.Image with the roads layer.

    Returns:
        A single ee.Image containing all combined bands, or None if essential inputs are missing.
    """
    print("Combining input features...")
    
    # Check for essential inputs (can adjust based on importance)
    if trendfire_img is None:
        print("Error: TrendFire outputs are missing. Cannot combine features.")
        return None
    # Add checks for other inputs if they are considered essential
    # if nicfi_trends is None: ... etc.

    # Start with the primary trend image
    combined = trendfire_img.toFloat() # Ensure consistent type

    # Add other layers conditionally
    if nicfi_trends is not None:
        print("  Adding NICFI trends...")
        combined = combined.addBands(nicfi_trends.toFloat())
    else:
        print("  Warning: Skipping NICFI trends (input was None).")

    if modis_trends is not None:
        print("  Adding MODIS LST trends...")
        combined = combined.addBands(modis_trends.toFloat())
    else:
        print("  Warning: Skipping MODIS LST trends (input was None).")

    if dem_slope is not None:
        print("  Adding DEM slope...")
        # Ensure DEM slope band has a unique name if needed (e.g., if 'slope' already exists)
        # dem_slope = dem_slope.select('slope').rename('dem_slope') 
        combined = combined.addBands(dem_slope.toFloat()) # Assumes dem_slope is single band 'slope'
    else:
        print("  Warning: Skipping DEM slope (input was None).")
        
    if roads_img is not None:
        print("  Adding Roads layer...")
        # Ensure roads image band has a unique name (e.g., rename if it's just 'b1')
        # roads_img = roads_img.rename('roads')
        combined = combined.addBands(roads_img.toFloat())
    else:
        print("  Warning: Skipping Roads layer (input was None).")

    print("Features combined successfully.")
    # print("Combined bands:", combined.bandNames().getInfo()) # Optional: Verify bands
    return combined

def resample_features(feature_image, scale, crs, boundary):
    """
    Resamples the combined feature image to a target CRS and scale.
    Uses reduceResolution followed by reproject, matching the JS approach.
    
    Args:
        feature_image: ee.Image with all combined feature bands.
        scale: Target resolution in meters (e.g., 30).
        crs: Target CRS string (e.g., 'EPSG:3857').
        boundary: ee.Geometry for clipping the final output.
        
    Returns:
        A resampled ee.Image or None on error.
    """
    print(f"Resampling features to {scale}m scale, CRS {crs}...")
    if feature_image is None:
        print("Error: Input feature_image is None. Cannot resample.")
        return None
        
    try:
        # Define the target projection
        target_projection = ee.Projection(crs).atScale(scale)
        
        # Resample using reduceResolution followed by reproject
        resampled_image = feature_image \
            .reduceResolution(reducer=ee.Reducer.mean(), maxPixels=1024) \
            .reproject(crs=target_projection)
            
        print("Resampling completed successfully.")
        return resampled_image.clip(boundary)
        
    except ee.EEException as e:
        print(f"GEE error during resampling: {e}")
        return None
    except Exception as e:
        print(f"Unexpected error during resampling: {e}")
        return None

def evaluate_classifier(classifier, validation_data, class_property):
    """
    Evaluates the classifier on validation data using an error matrix.

    Args:
        classifier: Trained ee.Classifier object.
        validation_data: ee.FeatureCollection used for validation.
        class_property: String, the name of the property containing the true class label.

    Returns:
        A dictionary containing accuracy, kappa, and the confusion matrix, or None on error.
    """
    print("Evaluating classifier...")
    if classifier is None:
        print("Error: Classifier object is missing.")
        return None
    if validation_data is None:
        print("Error: Validation data is missing.")
        return None
        
    try:
        # Classify the validation data
        validated = validation_data.classify(classifier)
        
        # Calculate the error matrix
        # Arguments are (actual, predicted)
        error_matrix = validated.errorMatrix(class_property, 'classification')
        
        # Extract metrics
        accuracy = error_matrix.accuracy()
        kappa = error_matrix.kappa()
        confusion = error_matrix.array() # Get the matrix as an ee.Array
        
        # Get results - use getInfo() to bring server-side objects to client
        accuracy_val = accuracy.getInfo()
        kappa_val = kappa.getInfo()
        confusion_val = confusion.getInfo()
        
        print(f"  Validation Overall Accuracy: {accuracy_val:.4f}")
        print(f"  Validation Kappa Coefficient: {kappa_val:.4f}")
        print("  Validation Confusion Matrix:")
        # Pretty print the matrix (assuming simple 2x2 or small matrix)
        if confusion_val:
            for row in confusion_val:
                print(f"    {row}")
        else:
            print("    Could not retrieve confusion matrix details.")
            
        return {
            'accuracy': accuracy_val,
            'kappa': kappa_val,
            'confusion_matrix': confusion_val
        }

    except ee.EEException as e:
        print(f"GEE error during classifier evaluation: {e}")
        return None
    except Exception as e:
        print(f"Unexpected error during classifier evaluation: {e}")
        return None

# --- Reusable Helper Functions (Consider placing in a shared utils.py) ---

def initialize_ee(project_id='ee-crop-health-telangana'): # Use the verified project ID
    """Initializes the Earth Engine API."""
    try:
        # ee.Authenticate() # Often needed only once per session/machine
        ee.Initialize(project=project_id, opt_url='https://earthengine-highvolume.googleapis.com')
        print(f"Earth Engine API initialized successfully for project: {project_id}")
    except ee.EEException as e:
        print(f"Error initializing Earth Engine API for project {project_id}: {e}")
        print("Please ensure:")
        print("  1. You have authenticated (e.g., via `gcloud auth application-default login`).")
        print(f"  2. The project ID '{project_id}' is correct.")
        print("  3. The Earth Engine API is enabled for this project in Google Cloud Console.")
        print("  4. The account used for authentication has access to this project.")
        raise # Re-raise the exception to stop execution if initialization fails
    except Exception as e:
        # Catch other potential errors during initialization
        print(f"An unexpected error occurred during Earth Engine initialization: {e}")
        raise

def get_boundary(asset_id='users/jonasnothnagel/pa_boundary'):
    """
    Get the boundary geometry from the GEE asset.
    Returns an ee.Geometry object.
    """
    try:
        print(f"Loading boundary from GEE asset: {asset_id}...")
        fc = ee.FeatureCollection(asset_id)
        if fc.size().getInfo() == 0:
            print(f"Error: Asset {asset_id} is empty or inaccessible.")
            return None
        boundary_geom = fc.geometry().dissolve()
        print("Successfully loaded boundary from GEE asset.")
        return boundary_geom
    except ee.EEException as e:
        print(f"Error loading boundary asset {asset_id}: {e}")
        return None
    except Exception as e:
        print(f"An unexpected error occurred loading the boundary asset: {e}")
        return None

def export_asset(image, asset_id, region, description=None, scale=30, crs='EPSG:3857'):
    """Exports an image to a GEE asset, deleting if exists."""
    if description is None:
        description = asset_id.split('/')[-1]
    # Check and delete existing asset
    try:
        # Use ee.data.deleteAsset for deleting
        asset_info = ee.data.getInfo(asset_id) # Check existence first
        if asset_info:
             print(f"Deleting existing asset: {asset_id}...")
             ee.data.deleteAsset(asset_id)
             print(f"Deleted existing asset: {asset_id}")
    except ee.EEException as e:
        # Catch specific exception for asset not found
        if 'not found' in str(e).lower():
            print(f"Asset {asset_id} not found. Proceeding...")
        else:
            # Re-raise other GEE errors during check/delete attempt
            print(f"GEE error checking/deleting asset {asset_id}: {e}. Proceeding attempt...")
    except Exception as e:
        # Catch unexpected errors during check/delete
        print(f"Unexpected error checking/deleting asset {asset_id}: {e}. Proceeding attempt...")
        
    print(f"Starting export task to Asset: {description} (ID: {asset_id})")
    task = ee.batch.Export.image.toAsset(
        image=image.toFloat(),
        description=description,
        assetId=asset_id,
        region=region,
        scale=scale,
        crs=crs,
        maxPixels=1e13
    )
    task.start()
    print(f"Task started (id: {task.id}). Check GEE Tasks.")
    return task

def export_table(collection, asset_id, description=None):
    """Exports a FeatureCollection to a GEE asset."""
    if description is None:
        description = asset_id.split('/')[-1]
    # Check and delete existing asset
    try:
        asset_info = ee.data.getInfo(asset_id) # Check existence first
        if asset_info:
            print(f"Deleting existing asset table: {asset_id}...")
            ee.data.deleteAsset(asset_id)
            print(f"Deleted existing asset table: {asset_id}")
    except ee.EEException as e:
        if 'not found' in str(e).lower():
             print(f"Asset table {asset_id} not found. Proceeding...")
        else:
            print(f"GEE error checking/deleting asset table {asset_id}: {e}. Proceeding...")
    except Exception as e:
        print(f"Unexpected error checking/deleting asset table {asset_id}: {e}. Proceeding...")
        
    print(f"Starting export task to Table Asset: {description} (ID: {asset_id})")
    task = ee.batch.Export.table.toAsset(
        collection=collection,
        description=description,
        assetId=asset_id,
    )
    task.start()
    print(f"Table export task started (id: {task.id}). Check GEE Tasks.")
    return task

if __name__ == "__main__":
    main()