utils.py

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

"""
Complete Enhanced utility functions for processing and visualizing EM-DAT disaster data for East Africa.
This module contains functions with enhanced direct EM-DAT processing methodology for complete 1990-2025 coverage.

FIXED ISSUES:
1. Flood calendar NaN conversion error - proper fillna() handling
2. Missing combined grid visualizations - fixed image creation functions
3. Column mismatch in frequency processing - proper aggregation methods
4. Complete 1990-2025 range support - fixed year ranges throughout
5. All functionality from 5 original scripts - comprehensive implementation

EXPECTED OUTPUTS:
- process_flood_data_for_calendar_enhanced() -> flood_calendar_plot_enhanced.png
- create_combined_image_enhanced() -> drought_extent_enhanced.png, flood_extent_enhanced.png  
- All frequency, calendar, and combined visualization functions working
"""

import os
import sys
import re
import glob
import ntpath
import pandas as pd
import numpy as np
from datetime import datetime
import calendar
import six
from ast import literal_eval
import matplotlib
matplotlib.use("Agg")  # Use non-interactive backend
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
from matplotlib import colors
from matplotlib.colors import LinearSegmentedColormap
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
import geopandas as gp
import PIL
from PIL import Image
import warnings

# Try to import optional dependencies
try:
    import geoplot as gplt
    import geoplot.crs as gcrs
    HAS_GEOPLOT = True
except ImportError:
    HAS_GEOPLOT = False
    warnings.warn("geoplot not available; using matplotlib/geopandas for plotting instead")

try:
    import mapclassify as mc
except ImportError:
    warnings.warn("mapclassify not available; some classification features may not work")
    mc = None

# Path Configuration (Enhanced)
BASE_DIR = ''
DATA_DIR = os.path.join(BASE_DIR, 'Impact Data')
SHAPEFILE_DIR = os.path.join(BASE_DIR, 'Impactshp')
OUTPUT_DIR = 'Output'

# Create output directory
os.makedirs(OUTPUT_DIR, exist_ok=True)
os.makedirs(os.path.join(OUTPUT_DIR, 'dr'), exist_ok=True)
os.makedirs(os.path.join(OUTPUT_DIR, 'fl'), exist_ok=True)

# Data file paths (Enhanced) - Updated to use new EMDAT file and proper geocoded mappings
EMDAT_DATA = 'public_emdat_custom_request_2026-01-21.xlsx'  # New EMDAT custom request file
DROUGHT_DATA = 'em_dat_ea_drought_adm2_v4.xlsx'  # Text-matched geocoding from 01-common_adm2_code_drought.py
FLOOD_DATA = 'em_dat_ea_flood_adm2_v2.xlsx'  # Text-matched geocoding from 02-common_adm2_code_flood.py
FLOOD_DATES_V0 = os.path.join(DATA_DIR, 'em_dat_flood_dates_v0.csv')
FLOOD_DATES_V1 = os.path.join(DATA_DIR, 'em_dat_flood_dates_v1.csv')

# Shapefile paths - Using available boundary files
ADM1_GEOJSON = 'icpac_adm1v3.geojson'  # Available admin1 geojson
ADM2_SHP = 'ea_adm2_geoboundaries.shp'  # GEOBoundaries ADM2 shapefile (NOW AVAILABLE)
EA_BACKGROUND_SHP = os.path.join(SHAPEFILE_DIR, 'ea_global_background.shp')
WATER_BODIES_SHP = os.path.join(SHAPEFILE_DIR, 'Africa_Water_Bodies.shp')

# Constants (Enhanced with additional countries - updated for new EMDAT format)
EA_COUNTRIES = [
    'Djibouti', 'Eritrea', 'Ethiopia', 'Kenya', 'Rwanda',
    'Somalia', 'South Sudan', 'Sudan (the)', 'Sudan',  # Both Sudan variants
    'Tanzania, United Republic of', 'United Republic of Tanzania',
    'Uganda', 'Burundi'  # Added Burundi
]

COUNTRY_CODES = {
    'Djibouti': 'DJI',
    'Eritrea': 'ERI',
    'Ethiopia': 'ETH',
    'Kenya': 'KEN',
    'Rwanda': 'RWA',
    'Somalia': 'SOM',
    'South Sudan': 'SSD',
    'Sudan (the)': 'SDN',
    'Sudan': 'SDN',  # Added for new EMDAT format
    'Tanzania, United Republic of': 'TZA',
    'United Republic of Tanzania': 'TZA',
    'Uganda': 'UGA',
    'Burundi': 'BDI'  # Added Burundi
}

###############################
# Helper Functions (Enhanced)
###############################

def extract_regions_from_location_enhanced(location_text):
    """Enhanced region extraction from location text with better parsing."""
    if not location_text or not isinstance(location_text, str) or location_text.lower() == 'nan':
        return []

    # Remove text in parentheses for initial splitting
    parentheses_text = []
    pattern = r'\(([^)]*)\)'
    for match in re.finditer(pattern, location_text):
        parentheses_text.append(match.group(1))

    # Remove parentheses for splitting
    clean_text = re.sub(pattern, '', location_text)

    # Split by common separators
    regions = []
    for separator in [',', ';', ' and ', '&']:
        if separator in clean_text:
            parts = [p.strip() for p in clean_text.split(separator) if p.strip()]
            regions.extend(parts)
            break

    # If no split occurred, use the whole text
    if not regions:
        regions = [clean_text.strip()]

    # Clean up region names
    cleaned_regions = []
    for region in regions:
        # Remove common qualifiers
        qualifiers = [
            'district', 'districts', 'region', 'regions', 'county', 'counties',
            'province', 'provinces', 'zone', 'zones', 'area', 'areas'
        ]
        for qualifier in qualifiers:
            region = re.sub(r'\b' + qualifier + r'\b', '', region, flags=re.IGNORECASE)

        region = re.sub(r'\s+', ' ', region).strip()

        if region:
            cleaned_regions.append(region)

    # Add regions from parentheses text
    for text in parentheses_text:
        if any(word in text.lower() for word in ['district', 'region', 'county', 'province', 'zone']):
            sub_regions = extract_regions_from_location_enhanced(text)
            cleaned_regions.extend(sub_regions)

    return cleaned_regions


def safe_convert_to_int_enhanced(value):
    """Enhanced safe conversion to integer handling NaN values."""
    if pd.isna(value):
        return None
    try:
        return int(float(value))
    except (ValueError, TypeError):
        return None


def safe_date_format_enhanced(value, default='01'):
    """Enhanced safe date formatting with NaN handling."""
    if pd.isna(value):
        return default
    try:
        return str(int(float(value))).zfill(2)
    except (ValueError, TypeError):
        return default


def flatten(l):
    """Flatten a list of lists into a single list."""
    return [item for sublist in l for item in sublist]


def return_colormap(classif):
    """Create colormap of matplotlib based on number of class and given colorcode."""
    c = matplotlib.colors.ColorConverter().to_rgb
    colorlist = [c("#ffffcc"), c("#ffeda0"), c("#fed976"), c("#ffb24b"), c("#fe8d3b"), c("#fd4e2a"), c("#e3181a")]
    color_code = colorlist
    c_cmap = LinearSegmentedColormap.from_list("my_colormap", color_code, N=len(classif), gamma=1.0)
    return c_cmap


def safe_read_geometries_enhanced(shapefile_path):
    """Enhanced safely read geometries from a shapefile, filtering out None values."""
    try:
        reader = Reader(shapefile_path)
        geometries = []
        for geom in reader.geometries():
            if geom is not None:
                geometries.append(geom)
        return geometries
    except Exception as e:
        print(f"Error reading shapefile {shapefile_path}: {e}")
        return []


def create_daterange_list_enhanced(row):
    """FIXED: Enhanced date range creation with proper NaN handling."""
    try:
        # Safe conversion with NaN handling
        start_year = str(row['Start Year']) if pd.notna(row['Start Year']) else '1990'
        start_month = safe_date_format_enhanced(row['Start Month'], '01')
        start_day = safe_date_format_enhanced(row['Start Day'], '01')
        end_year = str(row['End Year']) if pd.notna(row['End Year']) else start_year
        end_month = safe_date_format_enhanced(row['End Month'], '12')
        end_day = safe_date_format_enhanced(row['End Day'], '31')

        start_date_str = f"{start_year}-{start_month}-{start_day}"
        end_date_str = f"{end_year}-{end_month}-{end_day}"

        # Create date range
        date_str_list = pd.date_range(start_date_str, end_date_str, freq='D').strftime("%Y-%m-%d").tolist()
        return date_str_list
    except Exception as e:
        print(f"Error in create_daterange_list_enhanced: {e}")
        print(f"Problematic row: {row}")
        return []


def last_day_month_enhanced(row):
    """Enhanced last day of month calculation with NaN handling."""
    try:
        end_year = int(row['End Year']) if pd.notna(row['End Year']) else 1990
        end_month = int(row['End Month']) if pd.notna(row['End Month']) else 12

        firstday_month = datetime(end_year, end_month, 1)
        lastday_month = firstday_month + pd.DateOffset(months=1) - pd.DateOffset(days=1)
        lastday_month_str = lastday_month.strftime("%d")
        return lastday_month_str
    except Exception as e:
        print(f"Error in last_day_month_enhanced: {e}")
        return "31"


def order_file_enhanced(file_list):
    """Enhanced file ordering by year for creating combined visualizations."""
    if not file_list:
        return []

    filepath = os.path.dirname(file_list[0])
    file_list1 = [ntpath.basename(item).split('.')[0] for item in file_list]
    file_list1.sort()
    sorted_file_list = [f'{filepath}/{item}.png' for item in file_list1]
    return sorted_file_list

###############################
# Data Loading and Processing Functions (Enhanced)
###############################

def load_emdat_data_complete(disaster_type=None):
    """Enhanced EM-DAT data loading with complete 1990-2025 coverage."""
    try:
        print(f"Reading enhanced EM-DAT data from {EMDAT_DATA}")
        # New EMDAT file doesn't need skiprows and has different column names
        dfb = pd.read_excel(EMDAT_DATA, engine='openpyxl')

        # Rename columns to match expected format (DisNo. -> Dis No)
        if 'DisNo.' in dfb.columns:
            dfb = dfb.rename(columns={'DisNo.': 'Dis No'})

        # Filter for specified disaster type(s)
        if disaster_type:
            if isinstance(disaster_type, list):
                dfb = dfb[dfb['Disaster Type'].isin(disaster_type)]
            else:
                dfb = dfb[dfb['Disaster Type'] == disaster_type]

        # Filter for East Africa countries (enhanced list)
        dfb_ea = dfb[dfb['Country'].isin(EA_COUNTRIES)]

        print(f"Loaded {len(dfb_ea)} {disaster_type if disaster_type else 'total'} events for EA (1990-2025)")
        return dfb_ea

    except Exception as e:
        print(f"Error loading EM-DAT data: {e}")
        import traceback
        traceback.print_exc()
        return pd.DataFrame()

###############################
# Enhanced Data Conversion Functions
###############################

def convert_drought_data_enhanced():
    """Enhanced drought data conversion - uses pre-geocoded file if available."""
    print("Checking for drought geocoding data...")

    # Check if properly geocoded file exists (from 01-common_adm2_code_drought.py)
    if os.path.exists(DROUGHT_DATA):
        print(f"Using existing text-matched geocoding: {DROUGHT_DATA}")
        try:
            df = pd.read_excel(DROUGHT_DATA, engine='openpyxl')
            print(f"Loaded {len(df)} drought events with proper geocoding")
            return df
        except Exception as e:
            print(f"Error loading existing file: {e}")

    # Fallback: Run the geocoding script or use placeholder
    print("WARNING: No proper geocoded file found. Running placeholder approach.")
    print("For better results, run: python 01-common_adm2_code_drought.py")

    # Load EM-DAT data
    dfb_dr_ea = load_emdat_data_complete(disaster_type='Drought')

    if dfb_dr_ea.empty:
        print("No drought data found")
        return pd.DataFrame()

    try:
        # Try to load admin boundaries
        adm = gp.read_file(ADM2_SHP)
        shape_col = 'shapeID'
        group_col = 'shapeGroup'
        print(f"Loaded {len(adm)} admin regions from {ADM2_SHP}")

        # Create mapping based on country (placeholder if text matching unavailable)
        drought_mappings = {}

        for idx, row in dfb_dr_ea.iterrows():
            dis_no = row['Dis No']
            country = row['Country']
            country_code = COUNTRY_CODES.get(country)

            if country_code:
                country_shapes = adm[adm[group_col] == country_code]
                region_ids = country_shapes[shape_col].tolist()

                if region_ids:
                    affected_count = min(3, len(region_ids))
                    drought_mappings[dis_no] = region_ids[:affected_count]

        # Create output DataFrame
        drought_data = []
        for dis_no, region_codes in drought_mappings.items():
            drought_data.append({
                'Dis No': dis_no,
                'geob_adm2_list': region_codes
            })

        df = pd.DataFrame(drought_data)
        return df

    except Exception as e:
        print(f"Error in drought conversion: {e}")
        import traceback
        traceback.print_exc()
        return pd.DataFrame()


def convert_flood_data_enhanced():
    """Enhanced flood data conversion - uses pre-geocoded file if available."""
    print("Checking for flood geocoding data...")

    # Check if properly geocoded file exists (from 02-common_adm2_code_flood.py)
    if os.path.exists(FLOOD_DATA):
        print(f"Using existing text-matched geocoding: {FLOOD_DATA}")
        try:
            df = pd.read_excel(FLOOD_DATA, engine='openpyxl')
            print(f"Loaded {len(df)} flood events with proper geocoding")
            return df
        except Exception as e:
            print(f"Error loading existing file: {e}")

    # Fallback: Run the geocoding script or use placeholder
    print("WARNING: No proper geocoded file found. Running placeholder approach.")
    print("For better results, run: python 02-common_adm2_code_flood.py")

    # Load EM-DAT data
    dfb_fl_ea = load_emdat_data_complete(disaster_type='Flood')

    if dfb_fl_ea.empty:
        print("No flood data found")
        return pd.DataFrame()

    try:
        # Try to load admin boundaries
        adm = gp.read_file(ADM2_SHP)
        shape_col = 'shapeID'
        group_col = 'shapeGroup'
        print(f"Loaded {len(adm)} admin regions from {ADM2_SHP}")

        # Create mapping based on country (placeholder if text matching unavailable)
        flood_mappings = {}

        for idx, row in dfb_fl_ea.iterrows():
            dis_no = row['Dis No']
            country = row['Country']
            country_code = COUNTRY_CODES.get(country)

            if country_code:
                country_shapes = adm[adm[group_col] == country_code]
                region_ids = country_shapes[shape_col].tolist()

                if region_ids:
                    affected_count = min(2, len(region_ids))
                    flood_mappings[dis_no] = region_ids[:affected_count]

        # Create output DataFrame
        flood_data = []
        for dis_no, region_codes in flood_mappings.items():
            flood_data.append({
                'Dis No': dis_no,
                'geob_adm2_list': region_codes
            })

        df = pd.DataFrame(flood_data)
        return df

    except Exception as e:
        print(f"Error in flood conversion: {e}")
        import traceback
        traceback.print_exc()
        return pd.DataFrame()

###############################
# Enhanced Visualization Functions
###############################

def create_map_enhanced(adm2, year, region_codes, disaster_type, shape_col='shapeID'):
    """Enhanced map creation for a specific year and disaster type."""
    try:
        # Create figure
        plt.figure(figsize=(8.0, 9.5))
        ax = plt.axes(projection=ccrs.PlateCarree())

        # Set map extent for East Africa
        x_min = 21.838949
        x_max = 51.415695
        y_min = -11.745695
        y_max = 23.145147
        ax.set_extent([x_min, x_max, y_min, y_max])

        # Add background features
        ax.add_feature(cfeature.OCEAN, facecolor=(0.5, 0.5, 0.5))
        ax.add_feature(cfeature.LAND, facecolor='lightgray')
        ax.add_feature(cfeature.BORDERS, linestyle='-', alpha=0.5)
        ax.add_feature(cfeature.COASTLINE, linewidth=0.5)

        # Add background shapes if available
        try:
            background_geoms = safe_read_geometries_enhanced(EA_BACKGROUND_SHP)
            if background_geoms:
                ax.add_geometries(background_geoms, ccrs.PlateCarree(), facecolor='None')

            water_geoms = safe_read_geometries_enhanced(WATER_BODIES_SHP)
            if water_geoms:
                ax.add_geometries(water_geoms, ccrs.PlateCarree(), facecolor='grey')
        except Exception as e:
            pass  # Use cartopy features instead

        # Determine shape column based on what's in the GeoDataFrame
        if shape_col not in adm2.columns:
            shape_col = 'GID_1' if 'GID_1' in adm2.columns else adm2.columns[0]

        # Add disaster regions if available
        if region_codes:
            regions = adm2[adm2[shape_col].isin(region_codes)]
            if not regions.empty:
                ax.add_geometries(regions['geometry'], ccrs.PlateCarree(),
                                edgecolor='red', facecolor='red', alpha=1)
                print(f"Added {len(regions)} regions to the {disaster_type} map for {year}")

        # Add year label
        ax.text(0.72, 0.03, str(year), fontsize=50, fontweight='bold', 
                ha='left', va='center', color='black', transform=ax.transAxes)

        # Save the map
        event_dir = 'dr' if disaster_type == 'Drought' else 'fl'
        event_label = 'dr' if disaster_type == 'Drought' else 'fl'
        output_file = os.path.join(OUTPUT_DIR, event_dir, f"{str(year)}_{event_label}.png")
        plt.savefig(output_file, transparent=False, bbox_inches='tight', dpi=150)
        print(f"Enhanced {disaster_type} map saved for {year}: {output_file}")

        plt.close()
    except Exception as e:
        print(f"Error creating enhanced map for {disaster_type} in {year}: {e}")


def create_year_maps_enhanced(disaster_type='both'):
    """FIXED: Enhanced yearly maps creation for 1990-2025 complete coverage."""
    print(f"Creating enhanced yearly maps for disaster type: {disaster_type}")

    # Load admin boundaries - try shapefile first, then geojson
    try:
        adm2 = gp.read_file(ADM2_SHP)
        shape_col = 'shapeID'
        print(f"Loaded {len(adm2)} admin2 regions from shapefile")
    except Exception as e:
        print(f"ADM2 shapefile not found ({e}), using ADM1 geojson instead")
        try:
            adm2 = gp.read_file(ADM1_GEOJSON)
            shape_col = 'GID_1'
            print(f"Loaded {len(adm2)} admin1 regions from geojson")
        except Exception as e2:
            print(f"Error loading admin boundaries: {e2}")
            return

    # FIXED: Complete year range 1990-2025 (36 years)
    year_list = list(range(1990, 2026))  # FIXED: Now includes 2025
    print(f"Processing {len(year_list)} years: {year_list[0]}-{year_list[-1]}")

    # Load mapped data for drought
    if disaster_type in ['drought', 'both']:
        try:
            drought_file = DROUGHT_DATA
            if not os.path.exists(drought_file):
                print("Converting drought data first...")
                convert_drought_data_enhanced()

            drought_df = pd.read_excel(drought_file)
            drought_df['geob_adm2_list'] = drought_df['geob_adm2_list'].apply(
                lambda x: eval(str(x)) if isinstance(x, str) else x)

            # Group by year
            drought_by_year = {}
            for idx, row in drought_df.iterrows():
                try:
                    dis_no = row['Dis No']
                    year = int(dis_no.split('-')[0])

                    if year not in drought_by_year:
                        drought_by_year[year] = []

                    region_codes = row['geob_adm2_list']
                    if isinstance(region_codes, list):
                        drought_by_year[year].extend(region_codes)

                except Exception as e:
                    print(f"Error processing drought entry {idx}: {e}")

            # Create drought maps for each year
            for year in year_list:
                region_codes = drought_by_year.get(year, [])
                create_map_enhanced(adm2, year, region_codes, 'Drought', shape_col)

        except Exception as e:
            print(f"Error processing enhanced drought maps: {e}")

    # Load mapped data for flood
    if disaster_type in ['flood', 'both']:
        try:
            flood_file = FLOOD_DATA
            if not os.path.exists(flood_file):
                print("Converting flood data first...")
                convert_flood_data_enhanced()

            flood_df = pd.read_excel(flood_file)
            flood_df['geob_adm2_list'] = flood_df['geob_adm2_list'].apply(
                lambda x: eval(str(x)) if isinstance(x, str) else x)

            # Group by year
            flood_by_year = {}
            for idx, row in flood_df.iterrows():
                try:
                    dis_no = row['Dis No']
                    year = int(dis_no.split('-')[0])

                    if year not in flood_by_year:
                        flood_by_year[year] = []

                    region_codes = row['geob_adm2_list']
                    if isinstance(region_codes, list):
                        flood_by_year[year].extend(region_codes)

                except Exception as e:
                    print(f"Error processing flood entry {idx}: {e}")

            # Create flood maps for each year
            for year in year_list:
                region_codes = flood_by_year.get(year, [])
                create_map_enhanced(adm2, year, region_codes, 'Flood', shape_col)

        except Exception as e:
            print(f"Error processing enhanced flood maps: {e}")

    print("Completed enhanced yearly maps creation")

###############################
# FIXED: Enhanced Calendar Visualization Functions
###############################

def process_drought_data_for_calendar_enhanced():
    """Enhanced drought data processing for calendar visualization (1990-2025)."""
    print("Processing enhanced drought data for calendar...")

    # Read directly from EM-DAT data
    dfb_dr_ea = load_emdat_data_complete(disaster_type='Drought')

    # Extract time data
    dfb_dr_ea_time = dfb_dr_ea[['Dis No', 'Start Year', 'Start Month', 'Start Day', 
                                'End Year', 'End Month', 'End Day']].copy()

    # FIXED: Fill missing values BEFORE any conversions
    dfb_dr_ea_time['Start Month'] = dfb_dr_ea_time['Start Month'].fillna(1.0)
    dfb_dr_ea_time['Start Day'] = dfb_dr_ea_time['Start Day'].fillna(1.0)
    dfb_dr_ea_time['End Month'] = dfb_dr_ea_time['End Month'].fillna(12.0)

    # Apply function to get end days
    for idx, row in dfb_dr_ea_time.iterrows():
        dfb_dr_ea_time.loc[idx, 'End Day'] = last_day_month_enhanced(row)

    # FIXED: Safe format date fields with NaN handling
    dfb_dr_ea_time['Start Year'] = dfb_dr_ea_time["Start Year"].astype(str)
    dfb_dr_ea_time['Start Month'] = dfb_dr_ea_time["Start Month"].apply(
        lambda x: safe_date_format_enhanced(x, '01'))
    dfb_dr_ea_time['End Year'] = dfb_dr_ea_time["End Year"].astype(str)
    dfb_dr_ea_time['End Month'] = dfb_dr_ea_time["End Month"].apply(
        lambda x: safe_date_format_enhanced(x, '12'))

    # Create date range lists
    date_str_list_cont = []
    for idx, row in dfb_dr_ea_time.iterrows():
        date_str_list = create_daterange_list_enhanced(row)
        date_str_list_cont.append(date_str_list)

    dfb_dr_ea_time['date_str_list'] = date_str_list_cont
    df_dr = dfb_dr_ea_time.explode('date_str_list')

    # Group by date
    df_t = df_dr.groupby('date_str_list').agg({'Dis No': lambda x: list(x)})

    # Save intermediate result
    df_t.to_csv(os.path.join(OUTPUT_DIR, 'drought_intermediate_enhanced.csv'))

    dft1 = df_t.reset_index()
    dft1['date'] = dft1['date_str_list'].apply(pd.to_datetime, errors='coerce')
    dft1 = dft1.dropna(subset=['date'])  # Remove invalid dates
    dft1['year'] = dft1['date'].dt.year
    dft1['days'] = dft1['date'].dt.dayofyear

    # Generate day-by-day event data by year
    def gen_list_days_enhanced(im_df, year):
        n_days = 366 if calendar.isleap(year) else 365
        df = pd.DataFrame({"year": year, "days": range(1, n_days + 1)})
        dfta = dft1[dft1['year'] == year]
        dfta['event'] = 1
        dftb = dfta[['Dis No', 'date', 'days', 'event']]
        dftc = pd.merge(df, dftb, on='days', how='left')
        dftd = dftc.fillna(0)
        event_list = dftd['event'].tolist()
        event_list.insert(0, str(year))
        return event_list

    # FIXED: Create event lists for each year (1990-2025)
    year_list = np.arange(1990, 2026, 1)  # FIXED: Now includes 2025
    el_cont = []
    for yearl in year_list:
        el = gen_list_days_enhanced(dft1, yearl)
        el_cont.append(el)

    # Create the output dataframe
    db = pd.DataFrame(el_cont)
    db1 = db.set_index(0)
    db2 = db1.fillna(0)
    cols = db2.columns
    db2[cols] = db2[cols].apply(lambda x: x.astype(int))

    # Save to CSV
    output_csv = os.path.join(OUTPUT_DIR, 'cal_data_drought_enhanced.csv')
    db2.to_csv(output_csv)
    print(f"Enhanced drought calendar data saved to {output_csv}")

    return db2


def process_flood_data_for_calendar_enhanced():
    """FIXED: Enhanced flood data processing for calendar visualization with proper NaN handling."""
    print("Processing enhanced flood data for calendar...")

    # FIXED: Read directly from EM-DAT data like drought processing
    dfb_fl_ea = load_emdat_data_complete(disaster_type='Flood')

    # Extract time information - using .copy() to avoid SettingWithCopyWarning
    dfb_fl_ea_time = dfb_fl_ea[['Dis No', 'Start Year', 'Start Month', 'Start Day', 
                                'End Year', 'End Month', 'End Day']].copy()

    # FIXED: Fill NaN values BEFORE any conversions to avoid "cannot convert float NaN to integer"
    dfb_fl_ea_time['Start Month'] = dfb_fl_ea_time['Start Month'].fillna(1.0)
    dfb_fl_ea_time['Start Day'] = dfb_fl_ea_time['Start Day'].fillna(1.0)
    dfb_fl_ea_time['End Month'] = dfb_fl_ea_time['End Month'].fillna(12.0)

    # Apply function to get end days with enhanced error handling
    for idx, row in dfb_fl_ea_time.iterrows():
        dfb_fl_ea_time.loc[idx, 'End Day'] = last_day_month_enhanced(row)

    # FIXED: Safe format date fields with proper NaN handling
    dfb_fl_ea_time['Start Year'] = dfb_fl_ea_time["Start Year"].astype(str)
    dfb_fl_ea_time['Start Month'] = dfb_fl_ea_time["Start Month"].apply(
        lambda x: safe_date_format_enhanced(x, '01'))
    dfb_fl_ea_time['End Year'] = dfb_fl_ea_time["End Year"].astype(str)
    dfb_fl_ea_time['End Month'] = dfb_fl_ea_time["End Month"].apply(
        lambda x: safe_date_format_enhanced(x, '12'))

    # Create date range lists with enhanced error handling
    date_str_list_cont = []
    for idx, row in dfb_fl_ea_time.iterrows():
        try:
            date_str_list = create_daterange_list_enhanced(row)
            date_str_list_cont.append(date_str_list)
        except Exception as e:
            print(f"Error processing row {idx}: {e}")
            date_str_list_cont.append([])

    dfb_fl_ea_time['date_str_list'] = date_str_list_cont
    db_time1 = dfb_fl_ea_time[['Dis No', 'date_str_list']]

    # Explode the date lists
    df_fl = db_time1.explode('date_str_list')

    # Group by date
    df_t = df_fl.groupby('date_str_list').agg({'Dis No': lambda x: list(x)})

    # Create dataframe with date info
    dft1 = df_t.reset_index()
    dft1['date'] = dft1['date_str_list'].apply(pd.to_datetime, errors='coerce')
    dft1 = dft1.dropna(subset=['date'])  # Remove invalid dates
    dft1['year'] = dft1['date'].dt.year
    dft1['days'] = dft1['date'].dt.dayofyear

    # Generate day-by-day event data by year
    def gen_list_days_enhanced(im_df, year):
        n_days = 366 if calendar.isleap(year) else 365
        df = pd.DataFrame({"year": year, "days": range(1, n_days + 1)})
        dfta = dft1[dft1['year'] == year]
        dfta['event'] = 1
        dftb = dfta[['Dis No', 'date', 'days', 'event']]
        dftc = pd.merge(df, dftb, on='days', how='left')
        dftd = dftc.fillna(0)
        event_list = dftd['event'].tolist()
        event_list.insert(0, str(year))
        return event_list

    # FIXED: Create event lists for each year (1990-2025)
    year_list = np.arange(1990, 2026, 1)  # FIXED: Now includes 2025
    el_cont = []
    for yearl in year_list:
        el = gen_list_days_enhanced(dft1, yearl)
        el_cont.append(el)

    # Create the output dataframe
    db = pd.DataFrame(el_cont)
    db1 = db.set_index(0)
    db2 = db1.fillna(0)
    cols = db2.columns
    db2[cols] = db2[cols].apply(lambda x: x.astype(int))

    # Save to CSV
    output_csv = os.path.join(OUTPUT_DIR, 'cal_data_flood_enhanced.csv')
    db2.to_csv(output_csv)
    print(f"FIXED: Enhanced flood calendar data saved to {output_csv}")

    return db2


def create_calendar_visualization_enhanced(data_frame, output_file, event_type="Events"):
    """Enhanced calendar visualization for disaster events."""
    print(f"Creating enhanced calendar visualization for {event_type}...")

    # Table styling helper functions
    def set_align_for_column(table, col, align="left"):
        cells = [key for key in table._cells if key[1] == col]
        for cell in cells:
            table._cells[cell]._loc = align

    def set_width_for_column(table, col, width):
        cells = [key for key in table._cells if key[1] == col]
        for cell in cells:
            table._cells[cell]._width = width

    def set_height_for_row(table, row, height):
        cells = [key for key in table._cells if key[0] == row]
        for cell in cells:
            table._cells[cell]._height = height

    def colorcell_enhanced(tablerows, tablecols, cellDict):
        allcells = [(x, y) for x in tablerows[1:] for y in tablecols[1:]]
        for alcls in allcells:
            try:
                cell_value = int(cellDict[alcls]._text.get_text())
                if cell_value == 0:
                    cellDict[alcls].set_facecolor('#FFFFFF')
                elif cell_value == 1:
                    cellDict[alcls].set_facecolor('#FF0000')
                else:
                    cellDict[alcls].set_facecolor('#961414')
            except (ValueError, KeyError):
                cellDict[alcls].set_facecolor('#FFFFFF')

    def removeaqivalue_enhanced(tablerows, tablecols, mpl_table):
        allcells = [(x, y) for x in tablerows[1:] for y in tablecols[1:]]
        for alcls in allcells:
            mpl_table._cells[alcls]._text.set_text('')

    def set_height_for_row_except_head(table, rowlist, height):
        cells_list = []
        for row in rowlist:
            cells = [key for key in table._cells if key[0] == row]
            cells_list.append(cells)
        for cells in cells_list:
            for cell in cells:
                table._cells[cell]._height = height

    # Function for table creation
    def render_mpl_table_enhanced(data, col_width=1.0, row_height=1.625, font_size=12,
                                header_color='#40466e', row_colors=['#f1f1f2', 'w'], edge_color='w',
                                bbox=[0, 0, 1, 1], header_columns=0, ax=None, **kwargs):

        max_cols = min(366, len(data.columns))  # Ensure we don't exceed available columns
        mpl_table = ax.table(cellText=data.values, bbox=bbox, 
                           colLabels=['' for _ in range(max_cols)], 
                           cellLoc='center', **kwargs)

        set_align_for_column(mpl_table, col=0, align="left")
        set_width_for_column(mpl_table, 0, 0.1)
        set_height_for_row(mpl_table, 0, 0.01)
        set_height_for_row_except_head(mpl_table, np.arange(1, len(data.index)), 0.012)

        mpl_table.auto_set_font_size(False)
        mpl_table.set_fontsize(font_size)

        cellDict = mpl_table.get_celld()
        tablerows = np.arange(0, len(data.index) + 1)
        tablecols = np.arange(0, max_cols)

        for k, cell in six.iteritems(mpl_table._cells):
            cell.set_linewidth(0)
            if k[0] == 0 or k[1] < header_columns:
                cell.set_text_props(weight='bold', color='w')
                cell.set_facecolor(header_color)
                colorcell_enhanced(tablerows, tablecols, cellDict)
            else:
                cell.set_facecolor(row_colors[k[0] % len(row_colors)])

        removeaqivalue_enhanced(tablerows, tablecols, mpl_table)
        return ax

    try:
        # Create the plot
        fig = plt.figure()
        fig.set_size_inches(28, 20)  # Enhanced size for better visibility
        table = fig.add_axes([0.08, 0.02, 0.55, 0.9], frame_on=False)
        table.xaxis.set_ticks_position('none')
        table.yaxis.set_ticks_position('none')
        table.set_xticklabels('')
        table.set_yticklabels('')

        render_mpl_table_enhanced(data_frame, header_columns=1, col_width=0.002, ax=table)

        # Add title
        fig.suptitle(f'Enhanced {event_type} Calendar (1990-2025)', 
                    fontsize=24, fontweight='bold', y=0.95)

        # Save the figure
        plt.savefig(output_file, dpi=300, bbox_inches='tight')
        plt.close()
        print(f"Enhanced calendar visualization saved to {output_file}")

    except Exception as e:
        print(f"Error creating enhanced calendar visualization: {e}")


def create_alternative_drought_analysis_enhanced():
    """Enhanced alternative drought analysis using pivot approach."""
    print("Creating enhanced alternative drought analysis...")

    try:
        # Read drought data
        dfb_dr_ea = load_emdat_data_complete(disaster_type='Drought')

        # Process dates similar to calendar processing
        dfb_dr_ea_time = dfb_dr_ea[['Dis No', 'Start Year', 'Start Month', 'Start Day', 
                                  'End Year', 'End Month', 'End Day']].copy()

        # Fill missing values
        dfb_dr_ea_time['Start Month'] = dfb_dr_ea_time['Start Month'].fillna(1.0)
        dfb_dr_ea_time['Start Day'] = dfb_dr_ea_time['Start Day'].fillna(1.0)
        dfb_dr_ea_time['End Month'] = dfb_dr_ea_time['End Month'].fillna(12.0)

        # Get end days
        for idx, row in dfb_dr_ea_time.iterrows():
            dfb_dr_ea_time.loc[idx, 'End Day'] = last_day_month_enhanced(row)

        # Format dates safely
        dfb_dr_ea_time['Start Year'] = dfb_dr_ea_time["Start Year"].astype(str)
        dfb_dr_ea_time['Start Month'] = dfb_dr_ea_time["Start Month"].apply(
            lambda x: safe_date_format_enhanced(x, '01'))
        dfb_dr_ea_time['End Year'] = dfb_dr_ea_time["End Year"].astype(str)
        dfb_dr_ea_time['End Month'] = dfb_dr_ea_time["End Month"].apply(
            lambda x: safe_date_format_enhanced(x, '12'))

        # Create date ranges
        date_str_list_cont = []
        for idx, row in dfb_dr_ea_time.iterrows():
            date_str_list = create_daterange_list_enhanced(row)
            date_str_list_cont.append(date_str_list)

        dfb_dr_ea_time['date_str_list'] = date_str_list_cont
        df_dr = dfb_dr_ea_time.explode('date_str_list')

        # Alternative approach: Create a pivoted version
        df_dr1 = df_dr.copy()
        df_dr1['dr'] = 1
        df_dr1['date'] = df_dr1['date_str_list'].apply(pd.to_datetime, errors='coerce')
        df_dr1 = df_dr1.dropna(subset=['date'])
        df_dr1['doy'] = df_dr1['date'].dt.dayofyear
        df_dr1['year'] = df_dr1['date'].dt.year
        df_dr2 = df_dr1[['year', 'doy', 'dr']]

        # Pivot table
        df = pd.pivot_table(df_dr2, values='dr', index=['year'], columns=['doy'])
        df1 = df.rename_axis(None)
        df2 = df1.fillna(0)

        # Convert to integers
        cols = df2.columns
        df2[cols] = df2[cols].apply(lambda x: x.astype(int))

        # Take a subset of the first 60 columns
        df3 = df2.iloc[:, 0:60]
        df4 = df3.reset_index()

        # Save to CSV
        alternative_output = os.path.join(OUTPUT_DIR, 'drought_alternative_pivot_enhanced.csv')
        df4.to_csv(alternative_output, index=False)
        print(f"Enhanced alternative drought analysis saved to {alternative_output}")

        return df4

    except Exception as e:
        print(f"Error in enhanced alternative drought analysis: {e}")
        return pd.DataFrame()

###############################
# FIXED: Enhanced Frequency Analysis Functions
###############################

def process_drought_data_for_frequency_enhanced():
    """FIXED: Enhanced drought frequency processing with column mismatch fix."""
    print("Processing enhanced drought data for frequency map...")

    try:
        # Read drought data
        dfb_dr = pd.read_excel(DROUGHT_DATA, engine='openpyxl')
        dfb_dr['adm2_list'] = dfb_dr.geob_adm2_list.apply(lambda x: literal_eval(str(x)))
        dfb_dr1 = dfb_dr[['Dis No', 'adm2_list']]

        # Expand the adm2_list column - this is key for correct frequency calculation
        dfb_dr2 = dfb_dr1.explode('adm2_list')
        dfb_dr2['year'] = dfb_dr2['Dis No'].str.split('-').str[0]

        # FIXED: Count occurrences by adm2_list with proper column naming
        df = dfb_dr2.groupby(['adm2_list']).agg({'Dis No': 'count'})  # FIXED: specific column name

        # Reset index and rename columns properly
        df1 = df.reset_index()
        df1.columns = ['shapeID', 'count1']
        df1['count2'] = df1['count1']  # Add second count column for compatibility

        # Read shapefile or geojson
        try:
            adm = gp.read_file(ADM2_SHP)
            shape_col = 'shapeID'
        except Exception:
            adm = gp.read_file(ADM1_GEOJSON)
            shape_col = 'GID_1'
            df1 = df1.rename(columns={'shapeID': 'GID_1'})

        # Merge data with geometries
        db = pd.merge(adm, df1, on=shape_col, how='right')

        # Filter out rows with None geometries
        db = db[db['geometry'].notna()]

        print(f"Enhanced drought frequency data processed: {len(db)} regions")
        return db

    except Exception as e:
        print(f"Error in enhanced drought frequency processing: {e}")
        import traceback
        traceback.print_exc()
        return None


def process_flood_data_for_frequency_enhanced():
    """FIXED: Enhanced flood frequency processing with column mismatch fix."""
    print("Processing enhanced flood data for frequency map...")

    try:
        # Read flood data
        dfb_fl = pd.read_excel(FLOOD_DATA, engine='openpyxl')
        dfb_fl['adm2_list'] = dfb_fl.geob_adm2_list.apply(lambda x: literal_eval(str(x)))
        dfb_fl1 = dfb_fl[['Dis No', 'adm2_list']]

        # Expand the adm2_list column - this is key for correct frequency calculation
        dfb_fl2 = dfb_fl1.explode('adm2_list')
        dfb_fl2['year'] = dfb_fl2['Dis No'].str.split('-').str[0]

        # FIXED: Count occurrences by adm2_list with proper column naming
        df = dfb_fl2.groupby(['adm2_list']).agg({'Dis No': 'count'})  # FIXED: specific column name

        # Reset index and rename columns properly
        df1 = df.reset_index()
        df1.columns = ['shapeID', 'count1']
        df1['count2'] = df1['count1']  # Add second count column for compatibility

        # Read shapefile or geojson
        try:
            adm = gp.read_file(ADM2_SHP)
            shape_col = 'shapeID'
        except Exception:
            adm = gp.read_file(ADM1_GEOJSON)
            shape_col = 'GID_1'
            df1 = df1.rename(columns={'shapeID': 'GID_1'})

        # Merge data with geometries
        db = pd.merge(adm, df1, on=shape_col, how='right')

        # Filter out rows with None geometries
        db = db[db['geometry'].notna()]

        print(f"Enhanced flood frequency data processed: {len(db)} regions")
        return db

    except Exception as e:
        print(f"Error in enhanced flood frequency processing: {e}")
        import traceback
        traceback.print_exc()
        return None


def create_map_for_analysis_enhanced(db, c_cmap, classif, event_name="Events", year="1990-2025"):
    """Enhanced choropleth map for event frequency analysis."""
    print(f"Creating enhanced {event_name} frequency map...")

    # Check if db has valid geometries
    if db.empty or 'geometry' not in db.columns:
        print("Error: GeoDataFrame has no geometries")
        return

    # Filter out rows with None geometries
    db = db[db['geometry'].notna()]

    if db.empty:
        print("Error: No valid geometries left after filtering")
        return

    # Ensure count columns exist
    for col in ['count1', 'count2']:
        if col not in db.columns:
            print(f"Error: Column '{col}' not found in GeoDataFrame")
            return

    try:
        # Create figure and axis
        plt.figure(figsize=(10, 12))
        ax = plt.axes(projection=ccrs.PlateCarree())

        # Set map extent for East Africa
        x_min = 21.838949
        x_max = 51.415695
        y_min = -11.745695
        y_max = 23.145147
        ax.set_extent([x_min, x_max, y_min, y_max])

        # Create classification scheme
        try:
            if mc is not None:
                scheme = mc.UserDefined(db['count2'], classif)
            else:
                scheme = None
        except Exception as e:
            print(f"Warning: Error creating classification scheme: {e}")
            scheme = None

        # Plot choropleth using geopandas (more reliable)
        try:
            db.plot(column='count2', ax=ax, cmap=c_cmap, 
                   transform=ccrs.PlateCarree(),
                   legend=True, edgecolor='none', linewidth=0.05)
        except Exception as e:
            print(f"Error creating choropleth: {e}")
            # Fallback to simple plot
            db.plot(ax=ax, color='lightgray', edgecolor='gray', 
                   transform=ccrs.PlateCarree())

        # Add features
        ax.add_feature(cfeature.OCEAN, facecolor=(0.5, 0.5, 0.5))

        # Safely add background shapes
        try:
            background_geoms = safe_read_geometries_enhanced(EA_BACKGROUND_SHP)
            if background_geoms:
                ax.add_geometries(background_geoms, ccrs.PlateCarree(), facecolor='None')

            water_geoms = safe_read_geometries_enhanced(WATER_BODIES_SHP)
            if water_geoms:
                ax.add_geometries(water_geoms, ccrs.PlateCarree(), facecolor='grey')
        except Exception as e:
            print(f"Warning: Error adding background shapes: {e}")

        # Add labels
        ax.text(0.62, 0.13, str(year), fontsize=24, fontweight='bold', 
               ha='left', va='center', color='black', transform=ax.transAxes)
        ax.text(0.62, 0.18, f'Enhanced {event_name} Events', fontsize=20, fontweight='bold', 
               ha='left', va='center', color='black', transform=ax.transAxes)

        # Save the figure
        output_png_file = os.path.join(OUTPUT_DIR, f"frequency_{event_name.lower()}_enhanced.png")
        plt.savefig(output_png_file, transparent=False, bbox_inches='tight', dpi=150)
        print(f"Enhanced {event_name} frequency map saved to {output_png_file}")

        plt.close()

    except Exception as e:
        print(f"Error creating enhanced frequency map: {e}")


def create_choropleth_map_enhanced(db, classif, event_name="Events", year="1990-2025"):
    """Enhanced choropleth map creation specifically for frequency analysis."""
    print(f"Creating enhanced choropleth map for {event_name}...")

    try:
        c_cmap = return_colormap(classif)

        # Create figure
        plt.figure(figsize=(12, 10))
        ax = plt.axes(projection=ccrs.PlateCarree())

        # Set map extent for East Africa
        x_min = 21.838949
        x_max = 51.415695
        y_min = -11.745695
        y_max = 23.145147
        ax.set_extent([x_min, x_max, y_min, y_max])

        # Plot the data
        if not db.empty and 'count2' in db.columns:
            db.plot(column='count2', ax=ax, cmap=c_cmap, 
                   transform=ccrs.PlateCarree(),
                   legend=True, edgecolor='black', linewidth=0.1)

        # Add map features
        ax.add_feature(cfeature.COASTLINE)
        ax.add_feature(cfeature.BORDERS)
        ax.add_feature(cfeature.OCEAN, facecolor='lightblue', alpha=0.5)

        # Add title
        plt.title(f'Enhanced {event_name} Frequency Map ({year})', 
                 fontsize=16, fontweight='bold', pad=20)

        # Save the figure
        output_file = os.path.join(OUTPUT_DIR, f"choropleth_{event_name.lower()}_frequency.png")
        plt.savefig(output_file, transparent=False, bbox_inches='tight', dpi=150)
        print(f"Enhanced choropleth map saved to {output_file}")

        plt.close()

    except Exception as e:
        print(f"Error creating enhanced choropleth map: {e}")

###############################
# FIXED: Enhanced Combined Visualization Functions
###############################

def create_combined_image_enhanced(event_type='fl', output_file=None):
    """FIXED: Enhanced combined image creation with proper 6x6 grid for 1990-2025."""
    if output_file is None:
        event_name = 'flood' if event_type == 'fl' else 'drought'
        output_file = f'{event_name}_extent_enhanced.png'

    print(f"Creating enhanced combined image for {event_type}...")

    # Find all individual map files
    file_list = glob.glob(os.path.join(OUTPUT_DIR, event_type, '*.png'))

    if not file_list:
        print(f"No images found for {event_type} in {os.path.join(OUTPUT_DIR, event_type)}")
        return

    print(f"Found {len(file_list)} images for {event_type}")

    # Extract years from filenames for better ordering
    year_file_dict = {}
    for file_path in file_list:
        filename = ntpath.basename(file_path)
        if '_' in filename:
            try:
                year = int(filename.split('_')[0])
                year_file_dict[year] = file_path
            except ValueError:
                print(f"Warning: Couldn't extract year from filename {filename}")

    # FIXED: Create ordered list from 1990 to 2025 (36 years total)
    ls = []
    for year in range(1990, 2026):  # FIXED: Now includes 2025
        if year in year_file_dict:
            ls.append(year_file_dict[year])
        else:
            print(f"Warning: No file found for year {year}")

    print(f"Ordered {len(ls)} files by year for enhanced visualization")

    # FIXED: Use 6x6 grid to accommodate 36 years (1990-2025)
    grid_size = (6, 6)
    total_cells = grid_size[0] * grid_size[1]

    try:
        # Create figure with subplots
        fig = plt.figure(figsize=(18, 18))  # Enhanced size for better visibility
        ax = [fig.add_subplot(grid_size[0], grid_size[1], i+1) for i in range(total_cells)]

        # Add each image to a subplot
        for i, (axi, img_path) in enumerate(zip(ax, ls)):
            try:
                img1 = mpimg.imread(img_path)
                axi.imshow(img1)
                axi.set_aspect('equal')
                axi.set_yticklabels([])
                axi.set_xticklabels([])
                axi.set_xticks([])
                axi.set_yticks([])

                # Add year label to each subplot
                year = 1990 + i
                axi.text(0.05, 0.95, str(year), transform=axi.transAxes, 
                        fontsize=12, fontweight='bold', 
                        bbox=dict(boxstyle="round,pad=0.3", facecolor="white", alpha=0.8))

            except Exception as e:
                print(f"Error loading image {img_path}: {e}")

        # Remove unused subplots
        for i in range(len(ls), total_cells):
            fig.delaxes(ax[i])

        # Add title
        event_title = 'Drought' if event_type == 'dr' else 'Flood'
        fig.suptitle(f'Enhanced {event_title} Impact Extent (1990-2025)', 
                    fontsize=28, fontweight='bold', y=0.98)

        # Add data source
        fig.text(0.5, 0.01, 'Enhanced Data Source: EM-DAT, https://www.emdat.be/', 
                ha='center', va='bottom', fontsize=12, style='italic')

        # Save the figure
        output_path = os.path.join(OUTPUT_DIR, output_file)
        plt.subplots_adjust(wspace=0.02, hspace=0.02, left=0.02, right=0.98, top=0.95, bottom=0.05)
        plt.savefig(output_path, bbox_inches='tight', facecolor='white', dpi=200)
        plt.close()
        print(f"Enhanced combined image saved to {output_path}")

    except Exception as e:
        print(f"Error creating enhanced combined image: {e}")


def stitch_images_PIL_enhanced(event_type='fl', output_file='enhanced_combined_PIL.png'):
    """Enhanced PIL-based image stitching for alternative combined visualization."""
    print(f"Creating enhanced PIL stitched image for {event_type}...")

    file_list = glob.glob(os.path.join(OUTPUT_DIR, event_type, '*.png'))

    if not file_list:
        print(f"No images found for {event_type}")
        return

    # Order files by year
    list_im = order_file_enhanced(file_list)

    try:
        imgs = [Image.open(i) for i in list_im]

        # Pick the image which is the smallest, and resize the others to match it
        min_shape = sorted([(np.sum(i.size), i.size) for i in imgs])[0][1]

        # FIXED: Create 6 rows of 6 images each for 36 years (1990-2025)
        rows = []
        for row_start in range(0, len(imgs), 6):
            row_end = min(row_start + 6, len(imgs))
            row_imgs = [np.asarray(i.resize(min_shape)) for i in imgs[row_start:row_end]]

            # Pad row if necessary
            while len(row_imgs) < 6:
                # Create a blank image of the same size
                blank = np.ones_like(row_imgs[0]) * 255  # White background
                row_imgs.append(blank)

            row_combined = np.hstack(row_imgs)
            rows.append(row_combined)

        # Stack all rows vertically
        if rows:
            imgs_comb0 = np.vstack(rows)

            # Save the combined image
            imgs_comb = Image.fromarray(imgs_comb0.astype('uint8'))
            output_path = os.path.join(OUTPUT_DIR, output_file)
            imgs_comb.save(output_path)
            print(f"Enhanced PIL stitched image saved to {output_path}")
        else:
            print("No valid rows to combine")

    except Exception as e:
        print(f"Error in enhanced PIL image stitching: {e}")

###############################
# Enhanced Summary and Reporting Functions
###############################

def generate_frequency_analysis_summary():
    """Generate enhanced frequency analysis summary."""
    print("Generating enhanced frequency analysis summary...")

    try:
        summary_content = []
        summary_content.append("Enhanced Frequency Analysis Summary (1990-2025)")
        summary_content.append("=" * 50)
        summary_content.append("")

        # Check drought frequency data
        try:
            drought_db = process_drought_data_for_frequency_enhanced()
            if drought_db is not None:
                drought_regions = len(drought_db)
                drought_max_freq = drought_db['count1'].max() if 'count1' in drought_db.columns else 0
                summary_content.append(f"Drought Analysis:")
                summary_content.append(f"  Affected regions: {drought_regions}")
                summary_content.append(f"  Maximum frequency: {drought_max_freq}")
            else:
                summary_content.append("Drought Analysis: No data available")
        except Exception as e:
            summary_content.append(f"Drought Analysis: Error - {e}")

        summary_content.append("")

        # Check flood frequency data
        try:
            flood_db = process_flood_data_for_frequency_enhanced()
            if flood_db is not None:
                flood_regions = len(flood_db)
                flood_max_freq = flood_db['count1'].max() if 'count1' in flood_db.columns else 0
                summary_content.append(f"Flood Analysis:")
                summary_content.append(f"  Affected regions: {flood_regions}")
                summary_content.append(f"  Maximum frequency: {flood_max_freq}")
            else:
                summary_content.append("Flood Analysis: No data available")
        except Exception as e:
            summary_content.append(f"Flood Analysis: Error - {e}")

        # Write summary
        output_file = os.path.join(OUTPUT_DIR, 'frequency_analysis_results.txt')
        with open(output_file, 'w') as f:
            f.write('\n'.join(summary_content))

        print(f"Enhanced frequency analysis summary saved to {output_file}")
        return summary_content

    except Exception as e:
        print(f"Error generating enhanced frequency analysis summary: {e}")
        return []


def generate_yearly_frequency_summary():
    """Generate enhanced yearly frequency summary."""
    print("Generating enhanced yearly frequency summary...")

    try:
        summary_content = []
        summary_content.append("Enhanced Yearly Frequency Summary (1990-2025)")
        summary_content.append("=" * 50)
        summary_content.append("")

        # Analyze drought events by year
        try:
            dfb_dr_ea = load_emdat_data_complete(disaster_type='Drought')
            if not dfb_dr_ea.empty:
                dfb_dr_ea['Year'] = dfb_dr_ea['Start Year']
                yearly_drought = dfb_dr_ea.groupby('Year').size()

                summary_content.append("Drought Events by Year:")
                for year in range(1990, 2026):
                    count = yearly_drought.get(year, 0)
                    summary_content.append(f"  {year}: {count} events")
                summary_content.append("")
        except Exception as e:
            summary_content.append(f"Drought yearly analysis error: {e}")

        # Analyze flood events by year
        try:
            dfb_fl_ea = load_emdat_data_complete(disaster_type='Flood')
            if not dfb_fl_ea.empty:
                dfb_fl_ea['Year'] = dfb_fl_ea['Start Year']
                yearly_flood = dfb_fl_ea.groupby('Year').size()

                summary_content.append("Flood Events by Year:")
                for year in range(1990, 2026):
                    count = yearly_flood.get(year, 0)
                    summary_content.append(f"  {year}: {count} events")
        except Exception as e:
            summary_content.append(f"Flood yearly analysis error: {e}")

        # Write summary
        output_file = os.path.join(OUTPUT_DIR, 'yearly_frequency_summary.txt')
        with open(output_file, 'w') as f:
            f.write('\n'.join(summary_content))

        print(f"Enhanced yearly frequency summary saved to {output_file}")
        return summary_content

    except Exception as e:
        print(f"Error generating enhanced yearly frequency summary: {e}")
        return []