EEMT Calculation Workflows¶

#!/usr/bin/env python3
"""
Traditional EEMT Calculation
Based on Rasmussen et al. (2005, 2014)
"""

import numpy as np
import rasterio
import pandas as pd
from datetime import datetime

def calculate_pet_hamon(temp_mean, temp_max, temp_min, daylight_hours):
    """
    Calculate potential evapotranspiration using Hamon's equation

    Parameters:
    temp_mean: mean monthly temperature [°C]
    temp_max, temp_min: daily temperature extremes [°C]
    daylight_hours: day length [hours]

    Returns:
    PET in mm/month
    """

    # Saturated vapor pressure (Tetens equation)
    es = 0.6108 * np.exp(17.27 * temp_mean / (temp_mean + 237.3))  # kPa

    # Hamon PET equation
    pet_daily = 0.55 * (daylight_hours / 12) * (es / (temp_mean + 273.15)) * 25.4

    # Convert to monthly (approximate)
    days_in_month = 30.4  # Average
    pet_monthly = pet_daily * days_in_month

    return pet_monthly

def calculate_npp_lieth(temperature, precipitation, pet):
    """
    Calculate NPP using Lieth (1975) temperature-based approach

    Parameters:
    temperature: mean monthly temperature [°C]
    precipitation: monthly precipitation [mm]
    pet: potential evapotranspiration [mm]

    Returns:
    NPP in kg/m²/yr
    """

    # Only calculate NPP for months with water surplus
    npp_monthly = np.zeros_like(temperature)

    for i, (temp, precip, evap) in enumerate(zip(temperature, precipitation, pet)):
        if precip > evap and temp > 0:  # Growing conditions
            # Lieth equation: NPP = 3000[1 - exp(1.315 - 0.119T)]^-1
            npp_monthly[i] = 3000 * (1 - np.exp(1.315 - 0.119 * temp))**(-1)
            npp_monthly[i] *= (precip - evap) / precip  # Scale by water availability
        else:
            npp_monthly[i] = 0

    # Convert g/m²/month to kg/m²/yr  
    npp_annual = np.sum(npp_monthly) / 1000  # g to kg conversion

    return npp_annual

def calculate_eemt_traditional(temperature_data, precipitation_data, daylight_data):
    """
    Calculate traditional EEMT for each pixel

    Parameters:
    temperature_data: dict with 'tmin', 'tmax', 'tmean' monthly arrays
    precipitation_data: monthly precipitation array [mm]
    daylight_data: monthly daylight hours array [hours]

    Returns:
    EEMT array in MJ/m²/yr
    """

    # Get array dimensions
    shape = temperature_data['tmean'].shape
    eemt_result = np.zeros(shape)

    # Process each pixel
    for i in range(shape[0]):
        for j in range(shape[1]):

            # Extract pixel time series
            temp_mean = temperature_data['tmean'][i, j, :]
            temp_max = temperature_data['tmax'][i, j, :]
            temp_min = temperature_data['tmin'][i, j, :]
            precip = precipitation_data[i, j, :]
            daylight = daylight_data[i, j, :]

            # Skip if any data is missing
            if np.any(np.isnan([temp_mean, precip])):
                eemt_result[i, j] = np.nan
                continue

            # Calculate PET
            pet = calculate_pet_hamon(temp_mean, temp_max, temp_min, daylight)

            # Calculate effective precipitation (F)
            effective_precip = np.maximum(0, precip - pet)  # mm/month

            # Convert to mass flux [kg/m²/s]
            seconds_per_month = 30.4 * 24 * 3600
            F = (effective_precip / 1000) / seconds_per_month  # kg/m²/s

            # Calculate E_PPT [W/m²]
            c_w = 4180  # J/kg/K
            delta_T = np.maximum(0, temp_mean - 0)  # °C above freezing
            E_PPT = F * c_w * delta_T

            # Calculate NPP
            npp_annual = calculate_npp_lieth(temp_mean, precip, pet)  # kg/m²/yr
            npp_flux = npp_annual / (365 * 24 * 3600)  # kg/m²/s

            # Calculate E_BIO [W/m²]
            h_BIO = 22e6  # J/kg
            E_BIO = npp_flux * h_BIO

            # Calculate EEMT [W/m²]
            eemt_flux = np.mean(E_BIO + E_PPT)  # Average over months

            # Convert to MJ/m²/yr
            eemt_result[i, j] = eemt_flux * 365 * 24 * 3600 / 1e6

    return eemt_result

# Example usage
def run_traditional_workflow(dem_file, climate_dir, output_file):
    """Complete traditional EEMT workflow"""

    # Load climate data (assumes NetCDF format)
    import xarray as xr

    # Load DAYMET data
    tmin = xr.open_dataset(f"{climate_dir}/tmin_monthly.nc")
    tmax = xr.open_dataset(f"{climate_dir}/tmax_monthly.nc") 
    precip = xr.open_dataset(f"{climate_dir}/prcp_monthly.nc")

    # Calculate mean temperature
    tmean = (tmin + tmax) / 2

    # Calculate daylight hours (simplified)
    # This should use actual solar geometry calculations
    daylight_hours = np.full_like(tmean, 12.0)  # Placeholder

    # Prepare data arrays
    temp_data = {
        'tmean': tmean.values,
        'tmax': tmax.values, 
        'tmin': tmin.values
    }

    # Calculate EEMT
    eemt = calculate_eemt_traditional(temp_data, precip.values, daylight_hours)

    # Save results
    with rasterio.open(dem_file) as dem_src:
        profile = dem_src.profile.copy()
        profile.update(dtype='float32', count=1)

        with rasterio.open(output_file, 'w', **profile) as dst:
            dst.write(eemt.astype(np.float32), 1)

    print(f"Traditional EEMT saved to: {output_file}")

# Run workflow
# run_traditional_workflow('dem.tif', 'climate_data/', 'eemt_traditional.tif')

#!/usr/bin/env python3
"""
Topographic EEMT Calculation
Based on Rasmussen et al. (2014) methodology
"""

import numpy as np
import rasterio
import subprocess
import tempfile
import os
from pathlib import Path

class TopographicEEMT:
    """Calculate EEMT with topographic controls"""

    def __init__(self, dem_path, climate_dir, output_dir):
        self.dem_path = Path(dem_path)
        self.climate_dir = Path(climate_dir)
        self.output_dir = Path(output_dir)
        self.output_dir.mkdir(parents=True, exist_ok=True)

    def calculate_solar_radiation(self):
        """Calculate topographically-modified solar radiation"""

        # Run annual solar radiation calculation
        from grass_solar_calculator import GrassSolarCalculator

        solar_calc = GrassSolarCalculator(
            str(self.dem_path), 
            str(self.output_dir / 'solar')
        )

        # Calculate for full year
        solar_calc.calculate_annual_solar()
        solar_calc.calculate_monthly_summaries()

        print("✓ Solar radiation calculation completed")

    def calculate_topographic_temperature(self, base_temp, lapse_rate=-6.5):
        """
        Calculate topographically-modified temperature
        Following Eq. 6 from Rasmussen et al. (2014)
        """

        with rasterio.open(self.dem_path) as dem_src:
            elevation = dem_src.read(1)
            profile = dem_src.profile

        # Load solar radiation ratio (S_topo / S_flat)
        solar_ratio_file = self.output_dir / 'solar' / 'solar_ratio_annual.tif'

        if solar_ratio_file.exists():
            with rasterio.open(solar_ratio_file) as src:
                solar_ratio = src.read(1)
        else:
            print("Warning: Solar ratio not found, using elevation only")
            solar_ratio = np.ones_like(elevation)

        # Calculate temperature modification
        # T_i = T_b - T_lapse * (z_i - z_b)/1000 + C * (S_i - 1/S_i) * (1 - LAI_i/LAI_max)

        base_elevation = np.nanmin(elevation)
        elevation_diff = (elevation - base_elevation) / 1000  # km

        # Temperature lapse rate effect
        temp_lapse_effect = lapse_rate * elevation_diff

        # Solar radiation effect (simplified - no LAI for TOPO method)
        solar_effect = 2.0 * (solar_ratio - 1/solar_ratio)  # C=2 constant

        # Modified temperature
        temp_modified = base_temp - temp_lapse_effect + solar_effect

        return temp_modified, profile

    def calculate_mcwi(self):
        """
        Calculate Mass Conservative Wetness Index
        Following Rasmussen et al. (2014) Eq. 9-10
        """

        # Use GRASS to calculate flow accumulation and slope
        temp_location = tempfile.mkdtemp()

        grass_commands = f"""
# Import DEM
r.in.gdal input={self.dem_path} output=dem

# Calculate flow accumulation using D-infinity
r.terraflow elevation=dem filled=dem_filled direction=flow_dir \\
            swatershed=watersheds accumulation=flow_accum tci=twi

# Calculate slope in degrees
r.slope.aspect elevation=dem slope=slope_deg

# Calculate traditional wetness index  
r.mapcalc "wetness_index = log(flow_accum / tan(slope_deg * 3.14159/180))"

# Calculate mass conservative wetness index (MCWI)
# Normalize by mean wetness index to conserve mass
r.univar wetness_index
"""

        # Execute GRASS commands and calculate MCWI
        # (Implementation details for MCWI calculation)

        return self.output_dir / 'mcwi.tif'

    def calculate_effective_precipitation(self, precipitation, temperature):
        """
        Calculate effective precipitation with topographic redistribution
        Using Penman-Monteith and Budyko curve approach
        """

        # Load MCWI for water redistribution
        with rasterio.open(self.output_dir / 'mcwi.tif') as src:
            mcwi = src.read(1)

        # Calculate PET using Penman-Monteith (simplified)
        # This is a placeholder - full implementation needs wind, humidity, radiation
        pet = self.calculate_penman_monteith(temperature, precipitation)

        # Calculate AET using Budyko curve
        aridity_index = pet / precipitation
        w = 2.63  # Empirical constant

        # Zhang-Budyko equation
        aet_ratio = (1 + w * aridity_index) / (1 + w * aridity_index + 1/aridity_index)
        aet = precipitation * aet_ratio

        # Effective precipitation
        effective_precip = precipitation - aet

        # Redistribute using MCWI
        effective_precip_redistributed = effective_precip * mcwi

        return effective_precip_redistributed

    def calculate_npp_topographic(self, elevation, aspect, slope):
        """
        Calculate NPP with topographic controls
        Following Eq. 11 from Rasmussen et al. (2014)
        """

        # Calculate northness
        aspect_rad = np.deg2rad(aspect)
        slope_rad = np.deg2rad(slope)
        northness = np.cos(aspect_rad) * np.sin(slope_rad)

        # Empirical relationship from Whittaker and Niering (1975)
        # NPP = 0.39z + 346n - 187 [g/m²/yr]
        npp = 0.39 * elevation + 346 * northness - 187

        # Set minimum NPP
        npp = np.maximum(npp, 100)  # g/m²/yr minimum

        # Convert to kg/m²/yr
        npp_kg = npp / 1000

        return npp_kg

    def run_complete_workflow(self):
        """Execute complete topographic EEMT workflow"""

        print("Starting Topographic EEMT Calculation...")

        # Step 1: Calculate solar radiation
        print("1. Calculating solar radiation...")
        self.calculate_solar_radiation()

        # Step 2: Calculate MCWI
        print("2. Calculating mass conservative wetness index...")
        mcwi_file = self.calculate_mcwi()

        # Step 3: Load climate data and DEM
        print("3. Loading input data...")
        with rasterio.open(self.dem_path) as src:
            elevation = src.read(1)
            profile = src.profile

        # Load climate data (implementation depends on data format)
        # This is a placeholder for actual climate data loading
        climate_data = self.load_climate_data()

        # Step 4: Calculate topographically modified temperature
        print("4. Calculating topographic temperature modification...")
        temp_modified, _ = self.calculate_topographic_temperature(
            climate_data['temperature'], 
            lapse_rate=-6.5
        )

        # Step 5: Calculate effective precipitation with redistribution
        print("5. Calculating effective precipitation...")
        effective_precip = self.calculate_effective_precipitation(
            climate_data['precipitation'],
            temp_modified
        )

        # Step 6: Calculate NPP with topographic effects
        print("6. Calculating topographic NPP...")

        # Load slope and aspect from DEM
        slope, aspect = self.calculate_slope_aspect()
        npp = self.calculate_npp_topographic(elevation, aspect, slope)

        # Step 7: Calculate EEMT components
        print("7. Calculating EEMT components...")

        # E_BIO calculation
        h_bio = 22e6  # J/kg
        npp_flux = npp / (365 * 24 * 3600)  # kg/m²/s
        e_bio = npp_flux * h_bio  # W/m²

        # E_PPT calculation  
        c_w = 4180  # J/kg/K
        delta_t = np.maximum(0, temp_modified - 273.15)  # K above freezing
        precip_flux = effective_precip / (30.4 * 24 * 3600)  # kg/m²/s (monthly avg)
        e_ppt = precip_flux * c_w * delta_t  # W/m²

        # Calculate EEMT
        eemt_flux = e_bio + e_ppt  # W/m²
        eemt_annual = eemt_flux * 365 * 24 * 3600 / 1e6  # MJ/m²/yr

        # Step 8: Save results
        print("8. Saving results...")

        # Save EEMT
        with rasterio.open(self.output_dir / 'eemt_topographic.tif', 'w', **profile) as dst:
            dst.write(eemt_annual.astype(np.float32), 1)

        # Save components
        with rasterio.open(self.output_dir / 'e_bio_topographic.tif', 'w', **profile) as dst:
            dst.write((e_bio * 365 * 24 * 3600 / 1e6).astype(np.float32), 1)

        with rasterio.open(self.output_dir / 'e_ppt_topographic.tif', 'w', **profile) as dst:
            dst.write((e_ppt * 365 * 24 * 3600 / 1e6).astype(np.float32), 1)

        print(f"✓ Topographic EEMT calculation completed")
        print(f"Results saved to: {self.output_dir}")

        return eemt_annual

# Usage example
if __name__ == '__main__':

    # Initialize calculator
    calculator = TopographicEEMT(
        dem_path='data/elevation/dem.tif',
        climate_dir='data/climate/',
        output_dir='results/eemt_topographic/'
    )

    # Run complete workflow
    eemt_result = calculator.run_complete_workflow()

    # Print summary statistics
    print(f"\nEEMT Summary Statistics:")
    print(f"  Mean: {np.nanmean(eemt_result):.2f} MJ/m²/yr")
    print(f"  Min:  {np.nanmin(eemt_result):.2f} MJ/m²/yr") 
    print(f"  Max:  {np.nanmax(eemt_result):.2f} MJ/m²/yr")
    print(f"  Std:  {np.nanstd(eemt_result):.2f} MJ/m²/yr")

#!/usr/bin/env python3
"""
Vegetation-Enhanced EEMT Calculation  
Based on Rasmussen et al. (2014) EEMT_TOPO-VEG approach
"""

import numpy as np
import rasterio
from scipy import ndimage

class VegetationEEMT(TopographicEEMT):
    """EEMT calculation with full vegetation integration"""

    def __init__(self, dem_path, climate_dir, output_dir, vegetation_data=None):
        super().__init__(dem_path, climate_dir, output_dir)
        self.vegetation_data = vegetation_data

    def calculate_lai_from_ndvi(self, ndvi_file):
        """
        Calculate Leaf Area Index from NDVI
        Using Qi et al. (2000) polynomial for semiarid regions
        """

        with rasterio.open(ndvi_file) as src:
            ndvi = src.read(1)
            profile = src.profile

        # Qi et al. (2000) polynomial: LAI = ax³ + bx² + cx + d
        a, b, c, d = 18.99, -15.24, 6.124, -0.352
        lai = a * ndvi**3 + b * ndvi**2 + c * ndvi + d

        # Constrain LAI to reasonable range
        lai = np.clip(lai, 0, 10)

        # Save LAI
        lai_file = self.output_dir / 'lai.tif'
        with rasterio.open(lai_file, 'w', **profile) as dst:
            dst.write(lai.astype(np.float32), 1)

        return lai, lai_file

    def calculate_canopy_height_from_lidar(self, lidar_file):
        """
        Extract canopy height from LiDAR data
        """

        # This would process LiDAR point clouds or canopy height models
        # For now, return placeholder based on LAI
        with rasterio.open(self.output_dir / 'lai.tif') as src:
            lai = src.read(1)
            profile = src.profile

        # Estimate canopy height from LAI (simplified relationship)
        # In practice, use actual LiDAR processing
        canopy_height = lai * 2.5  # Rough approximation

        # Save canopy height
        height_file = self.output_dir / 'canopy_height.tif'
        with rasterio.open(height_file, 'w', **profile) as dst:
            dst.write(canopy_height.astype(np.float32), 1)

        return canopy_height, height_file

    def calculate_npp_vegetation(self, canopy_height):
        """
        Calculate NPP from canopy height
        Following Eq. 12 from Rasmussen et al. (2014)
        """

        # Polynomial relationship: NPP = 196 + 36h - 0.61h² - 12.09h³
        h = canopy_height
        npp = 196 + 36*h - 0.61*h**2 - 12.09*h**3

        # Set minimum NPP
        npp = np.maximum(npp, 100)  # g/m²/yr

        # Convert to kg/m²/yr
        npp_kg = npp / 1000

        return npp_kg

    def calculate_surface_resistance(self, lai):
        """
        Calculate surface resistance from LAI
        Following Schulze et al. (1994) and Kelliher et al. (1995)
        """

        # Maximum leaf stomatal conductance 
        g_max = 0.008  # m/s

        # Bulk surface conductance from LAI
        # Polynomial fit to literature data
        g_surface = g_max * (1 - np.exp(-0.5 * lai))

        # Surface resistance (inverse of conductance)
        r_surface = 1 / np.maximum(g_surface, 1e-6)  # Avoid division by zero

        # Constrain to reasonable range
        r_surface = np.clip(r_surface, 38, 1000)  # s/m

        return r_surface

    def calculate_aet_penman_monteith(self, temperature, humidity, wind_speed, 
                                    net_radiation, lai):
        """
        Calculate actual evapotranspiration using full Penman-Monteith
        Including surface and aerodynamic resistance
        """

        # Calculate surface resistance
        r_surface = self.calculate_surface_resistance(lai)

        # Calculate aerodynamic resistance (simplified)
        # In practice, use canopy height and wind profile
        canopy_height = lai * 2.0  # Rough estimate
        r_aero = 208 / np.maximum(wind_speed, 0.1) * np.log(2.0 / (0.1 * canopy_height))
        r_aero = np.clip(r_aero, 10, 500)  # s/m

        # Psychrometric constant
        gamma = 0.665  # kPa/°C

        # Slope of saturation vapor pressure curve
        delta = 4098 * (0.6108 * np.exp(17.27 * temperature / (temperature + 237.3))) / (temperature + 237.3)**2

        # Vapor pressure deficit
        es = 0.6108 * np.exp(17.27 * temperature / (temperature + 237.3))
        ea = humidity * es / 100  # Assuming humidity is relative humidity %
        vpd = es - ea

        # Penman-Monteith equation
        numerator = delta * net_radiation + gamma * 900 * vpd / (temperature + 273) / r_aero
        denominator = delta + gamma * (1 + r_surface / r_aero)

        aet = numerator / denominator  # mm/day

        return aet

    def run_vegetation_workflow(self, ndvi_file=None, lidar_file=None):
        """Execute complete vegetation EEMT workflow"""

        print("Starting Vegetation EEMT Calculation...")

        # Step 1: Calculate solar radiation (inherited)
        print("1. Calculating solar radiation...")
        self.calculate_solar_radiation()

        # Step 2: Process vegetation data
        print("2. Processing vegetation data...")

        if ndvi_file:
            lai, lai_file = self.calculate_lai_from_ndvi(ndvi_file)
        else:
            print("Warning: No NDVI data provided, using default LAI")
            lai = np.ones((100, 100)) * 2.0  # Placeholder

        if lidar_file:
            canopy_height, height_file = self.calculate_canopy_height_from_lidar(lidar_file)
        else:
            canopy_height, height_file = self.calculate_canopy_height_from_lidar(None)

        # Step 3: Calculate vegetation-modified NPP
        print("3. Calculating vegetation NPP...")
        npp = self.calculate_npp_vegetation(canopy_height)

        # Step 4: Calculate AET with vegetation effects
        print("4. Calculating AET with vegetation controls...")

        # Load climate data
        climate_data = self.load_climate_data()

        # Calculate AET using Penman-Monteith with vegetation resistance
        aet = self.calculate_aet_penman_monteith(
            climate_data['temperature'],
            climate_data['humidity'], 
            climate_data['wind_speed'],
            climate_data['net_radiation'],
            lai
        )

        # Step 5: Calculate effective precipitation
        effective_precip = climate_data['precipitation'] - aet

        # Step 6: Apply topographic redistribution
        with rasterio.open(self.output_dir / 'mcwi.tif') as src:
            mcwi = src.read(1)

        effective_precip_redistributed = effective_precip * mcwi

        # Step 7: Calculate EEMT
        print("5. Calculating final EEMT...")

        # Load DEM for output profile
        with rasterio.open(self.dem_path) as src:
            profile = src.profile

        # Energy calculations
        h_bio = 22e6  # J/kg
        c_w = 4180   # J/kg/K

        # Convert fluxes to W/m²
        npp_flux = npp / (365 * 24 * 3600)  # kg/m²/s
        precip_flux = effective_precip_redistributed / (30.4 * 24 * 3600)  # kg/m²/s

        # Energy components
        e_bio = npp_flux * h_bio
        e_ppt = precip_flux * c_w * np.maximum(0, climate_data['temperature'] - 273.15)

        # Total EEMT
        eemt_flux = e_bio + e_ppt
        eemt_annual = eemt_flux * 365 * 24 * 3600 / 1e6  # MJ/m²/yr

        # Step 8: Save results
        print("6. Saving results...")

        outputs = {
            'eemt_vegetation.tif': eemt_annual,
            'e_bio_vegetation.tif': e_bio * 365 * 24 * 3600 / 1e6,
            'e_ppt_vegetation.tif': e_ppt * 365 * 24 * 3600 / 1e6,
            'npp_vegetation.tif': npp,
            'lai.tif': lai,
            'canopy_height.tif': canopy_height
        }

        for filename, data in outputs.items():
            output_path = self.output_dir / filename
            with rasterio.open(output_path, 'w', **profile) as dst:
                dst.write(data.astype(np.float32), 1)

        print(f"✓ Vegetation EEMT calculation completed")
        print(f"Results saved to: {self.output_dir}")

        return eemt_annual

# Command line interface
if __name__ == '__main__':
    import argparse

    parser = argparse.ArgumentParser(description='Vegetation EEMT Calculator')
    parser.add_argument('dem', help='Input DEM file')
    parser.add_argument('--climate-dir', required=True, help='Climate data directory')
    parser.add_argument('--output-dir', required=True, help='Output directory')
    parser.add_argument('--ndvi', help='NDVI raster file for LAI calculation')
    parser.add_argument('--lidar', help='LiDAR file for canopy height')

    args = parser.parse_args()

    # Run vegetation EEMT calculation
    calculator = VegetationEEMT(args.dem, args.climate_dir, args.output_dir)
    result = calculator.run_vegetation_workflow(args.ndvi, args.lidar)

    print(f"Vegetation EEMT range: {np.nanmin(result):.1f} - {np.nanmax(result):.1f} MJ/m²/yr")

#!/usr/bin/env python3
"""
Complete EEMT calculation pipeline integrating all three approaches
Enhanced version of eemt/run-workflow
"""

import argparse
import sys
from pathlib import Path
import numpy as np
import rasterio

def run_complete_eemt_pipeline(dem_file, output_dir, climate_dir, 
                             start_year=2015, end_year=2020,
                             vegetation_data=None, validation_data=None):
    """
    Run complete EEMT pipeline with all three calculation methods

    Parameters:
    dem_file: Path to elevation data
    output_dir: Output directory for results
    climate_dir: Directory containing climate data
    start_year, end_year: Analysis time period
    vegetation_data: Dict with 'ndvi' and 'lidar' file paths
    validation_data: Dict with validation datasets (soil depth, biomass, etc.)
    """

    print("=== EEMT Complete Pipeline ===")

    output_dir = Path(output_dir) 
    output_dir.mkdir(parents=True, exist_ok=True)

    results = {}

    # Method 1: Traditional EEMT
    print("\n1. Calculating Traditional EEMT...")

    try:
        eemt_trad = run_traditional_workflow(
            dem_file, 
            climate_dir, 
            output_dir / 'eemt_traditional.tif'
        )
        results['traditional'] = eemt_trad
        print("✓ Traditional EEMT completed")
    except Exception as e:
        print(f"✗ Traditional EEMT failed: {e}")

    # Method 2: Topographic EEMT  
    print("\n2. Calculating Topographic EEMT...")

    try:
        topo_calculator = TopographicEEMT(dem_file, climate_dir, output_dir / 'topographic')
        eemt_topo = topo_calculator.run_complete_workflow()
        results['topographic'] = eemt_topo
        print("✓ Topographic EEMT completed")
    except Exception as e:
        print(f"✗ Topographic EEMT failed: {e}")

    # Method 3: Vegetation EEMT
    print("\n3. Calculating Vegetation EEMT...")

    try:
        veg_calculator = VegetationEEMT(
            dem_file, climate_dir, output_dir / 'vegetation',
            vegetation_data
        )

        ndvi_file = vegetation_data.get('ndvi') if vegetation_data else None
        lidar_file = vegetation_data.get('lidar') if vegetation_data else None

        eemt_veg = veg_calculator.run_vegetation_workflow(ndvi_file, lidar_file)
        results['vegetation'] = eemt_veg
        print("✓ Vegetation EEMT completed")
    except Exception as e:
        print(f"✗ Vegetation EEMT failed: {e}")

    # Comparison Analysis
    print("\n4. Generating Comparison Analysis...")

    if len(results) > 1:
        generate_comparison_analysis(results, output_dir / 'comparison')

    # Validation (if data provided)
    if validation_data:
        print("\n5. Running Validation Analysis...")
        run_validation_analysis(results, validation_data, output_dir / 'validation')

    print(f"\n=== Pipeline Complete ===")
    print(f"Results saved to: {output_dir}")

    return results

def generate_comparison_analysis(results, output_dir):
    """Generate comparison plots and statistics between EEMT methods"""

    output_dir.mkdir(parents=True, exist_ok=True)

    # Load first result to get spatial structure
    first_key = list(results.keys())[0]

    if isinstance(results[first_key], str):
        # Results are file paths
        comparison_data = {}
        for method, filepath in results.items():
            with rasterio.open(filepath) as src:
                comparison_data[method] = src.read(1)
                profile = src.profile
    else:
        # Results are arrays
        comparison_data = results

    # Calculate difference maps
    if 'traditional' in comparison_data and 'topographic' in comparison_data:
        diff_topo_trad = comparison_data['topographic'] - comparison_data['traditional']

        with rasterio.open(output_dir / 'difference_topo_minus_trad.tif', 'w', **profile) as dst:
            dst.write(diff_topo_trad.astype(np.float32), 1)

    if 'vegetation' in comparison_data and 'topographic' in comparison_data:
        diff_veg_topo = comparison_data['vegetation'] - comparison_data['topographic'] 

        with rasterio.open(output_dir / 'difference_veg_minus_topo.tif', 'w', **profile) as dst:
            dst.write(diff_veg_topo.astype(np.float32), 1)

    # Summary statistics
    stats_file = output_dir / 'comparison_statistics.txt'
    with open(stats_file, 'w') as f:
        f.write("EEMT Method Comparison Statistics\\n")
        f.write("=" * 40 + "\\n\\n")

        for method, data in comparison_data.items():
            f.write(f"{method.upper()} EEMT:\\n")
            f.write(f"  Mean: {np.nanmean(data):.2f} MJ/m²/yr\\n")
            f.write(f"  Std:  {np.nanstd(data):.2f} MJ/m²/yr\\n")
            f.write(f"  Min:  {np.nanmin(data):.2f} MJ/m²/yr\\n")
            f.write(f"  Max:  {np.nanmax(data):.2f} MJ/m²/yr\\n\\n")

    print(f"✓ Comparison analysis saved to: {output_dir}")

def main():
    """Main command line interface"""

    parser = argparse.ArgumentParser(description='Complete EEMT Pipeline')
    parser.add_argument('dem', help='Input DEM file path')
    parser.add_argument('--output', '-o', required=True, help='Output directory')
    parser.add_argument('--climate', '-c', required=True, help='Climate data directory')
    parser.add_argument('--start-year', type=int, default=2015, help='Start year')
    parser.add_argument('--end-year', type=int, default=2020, help='End year')
    parser.add_argument('--ndvi', help='NDVI file for vegetation analysis')
    parser.add_argument('--lidar', help='LiDAR file for canopy height')
    parser.add_argument('--validation-soil', help='Soil depth data for validation')
    parser.add_argument('--validation-biomass', help='Biomass data for validation')

    args = parser.parse_args()

    # Prepare vegetation data
    vegetation_data = {}
    if args.ndvi:
        vegetation_data['ndvi'] = args.ndvi
    if args.lidar:
        vegetation_data['lidar'] = args.lidar

    # Prepare validation data  
    validation_data = {}
    if args.validation_soil:
        validation_data['soil_depth'] = args.validation_soil
    if args.validation_biomass:
        validation_data['biomass'] = args.validation_biomass

    # Run complete pipeline
    results = run_complete_eemt_pipeline(
        args.dem,
        args.output,
        args.climate,
        args.start_year,
        args.end_year,
        vegetation_data if vegetation_data else None,
        validation_data if validation_data else None
    )

    # Success summary
    success_count = len(results)
    print(f"\\nPipeline completed with {success_count}/3 methods successful")

    return success_count > 0

if __name__ == '__main__':
    success = main()
    sys.exit(0 if success else 1)

#!/usr/bin/env python3
"""
Automated data acquisition for EEMT calculations
Integrates with public data sources
"""

from data_sources.elevation import download_3dep, download_opentopo
from data_sources.climate import download_daymet_spatial, download_prism
from data_sources.satellite import download_landsat_ndvi

def setup_eemt_project(study_area, years, project_dir):
    """
    Automated setup of complete EEMT project with public data

    Parameters:
    study_area: [west, south, east, north] bounding box
    years: [start_year, end_year] 
    project_dir: output directory
    """

    project_dir = Path(project_dir)
    project_dir.mkdir(parents=True, exist_ok=True)

    print(f"Setting up EEMT project for {study_area}")
    print(f"Time period: {years[0]}-{years[1]}")
    print(f"Project directory: {project_dir}")

    # 1. Download elevation data
    print("\\n1. Downloading elevation data...")
    dem_file = download_3dep(study_area, resolution='10m')
    # Fallback to global data if US data unavailable  
    if not dem_file:
        dem_file = download_opentopo(study_area, 'SRTMGL1')

    # 2. Download climate data
    print("\\n2. Downloading climate data...")
    climate_files = download_daymet_spatial(
        study_area, range(years[0], years[1]+1), 
        ['tmin', 'tmax', 'prcp', 'vp']
    )

    # 3. Download vegetation data
    print("\\n3. Downloading vegetation data...")
    ndvi_file = download_landsat_ndvi(study_area, years[0])

    # 4. Set up analysis directories
    print("\\n4. Setting up analysis structure...")

    analysis_config = {
        'dem_file': dem_file,
        'climate_dir': project_dir / 'climate',
        'vegetation_data': {'ndvi': ndvi_file},
        'output_dir': project_dir / 'results'
    }

    # Save configuration
    import json
    with open(project_dir / 'eemt_config.json', 'w') as f:
        json.dump(analysis_config, f, indent=2, default=str)

    print("\\n✓ Project setup completed!")
    print(f"Configuration saved to: {project_dir}/eemt_config.json")
    print(f"Run analysis with: python run_complete_eemt.py {project_dir}/eemt_config.json")

    return analysis_config

# Example usage
if __name__ == '__main__':

    # Arizona Sky Islands study area
    bbox = [-111.0, 32.0, -110.0, 32.5]
    years = [2015, 2020]

    config = setup_eemt_project(bbox, years, 'arizona_eemt_project')