utils.py

#!/usr/bin/env python3
"""Geo utilities for viewing maps and merging data."""

import argparse
import difflib
from pathlib import Path

import yaml

# Load config
CONFIG_PATH = Path(__file__).parent / "data_config.yaml"

def load_config():
    """Load configuration from YAML file."""
    if not CONFIG_PATH.exists():
        raise FileNotFoundError(f"Config file not found: {CONFIG_PATH}")
    with open(CONFIG_PATH) as f:
        return yaml.safe_load(f)

CONFIG = load_config()

# Paths (from config)
_paths = CONFIG["paths"]
ADM0_GPKG = _paths["admin0_gpkg"]
ADM1_GPKG = _paths["admin1_gpkg"]
DEFAULT_CSV = _paths["admin1_csv"]
CITIES_SHP = _paths["cities_shp"]
DEFAULT_OUTPUT_DIR = _paths["output_dir"]

# Default level from config (admin0 or admin1)
DEFAULT_GPKG = ADM1_GPKG if CONFIG["default_level"] == "admin1" else ADM0_GPKG

# Column settings (from config)
_cols = CONFIG["columns"]
NAME_COLUMNS = _cols["name_search"]
LABEL_COLUMNS = _cols["labels"]
MERGE_KEYS = _cols["merge_keys"]
WB_REGION_COL = _cols.get("wb_region_col", "WB_REGION")
WB_REGIONS = _cols.get("wb_regions", {})

# Cities settings (from config)
_cities_cfg = CONFIG["cities"]
CITIES_NAME_COL = _cities_cfg["name_col"]
CITIES_POP_COL = _cities_cfg["population_col"]
CITIES_COUNTRY_COL = _cities_cfg["country_col"]
CITIES_TOP_N = _cities_cfg["top_n_default"]


def make_output_path(output_dir: str, gpkg_path: str, filter_term: str = None, suffix: str = "", ext: str = "html") -> Path:
    """Generate output path: {output_dir}/{source}_{filter}_{suffix}.{ext}"""
    out = Path(output_dir)
    out.mkdir(parents=True, exist_ok=True)
    source = Path(gpkg_path).stem.replace(" ", "_").replace("-", "_")[:20]
    parts = [source]
    if filter_term:
        parts.append(filter_term.lower().replace(" ", "_"))
    if suffix:
        parts.append(suffix)
    return out / f"{'_'.join(parts)}.{ext}"


def fuzzy_filter(gdf, search_terms, search_columns: list = None):
    """Filter GeoDataFrame by fuzzy matching search terms across name columns.

    Args:
        gdf: GeoDataFrame to filter
        search_terms: Single term (str) or multiple terms (list) to match
        search_columns: Columns to search in (default: NAME_COLUMNS from config)

    Returns:
        Filtered GeoDataFrame containing features matching any of the search terms
    """
    if not search_terms:
        return gdf

    # Normalize to list
    if isinstance(search_terms, str):
        search_terms = [search_terms]

    cols = search_columns or [c for c in NAME_COLUMNS if c in gdf.columns]
    if not cols:
        print(f"Warning: No name columns found. Available: {gdf.columns.tolist()}")
        return gdf

    # Build lookup: lowercase value -> (original value, column, row indices)
    lookup = {}
    for col in cols:
        for idx, val in gdf[col].dropna().items():
            key = str(val).lower()
            if key not in lookup:
                lookup[key] = (val, col, [])
            lookup[key][2].append(idx)

    # Match each search term and collect indices
    all_indices = set()
    for term in search_terms:
        matches = difflib.get_close_matches(term.lower(), lookup.keys(), n=1, cutoff=0.4)
        if not matches:
            print(f"No match found for '{term}'")
            continue

        matched_key = matches[0]
        matched_val, matched_col, indices = lookup[matched_key]
        score = difflib.SequenceMatcher(None, term.lower(), matched_key).ratio()
        print(f"Matched '{term}' → '{matched_val}' (column: {matched_col}, score: {score:.2f})")
        all_indices.update(indices)

    if not all_indices:
        print("No matches found. Returning full dataset.")
        return gdf

    return gdf.loc[list(all_indices)]


def filter_by_region(gdf, region_codes: list, adm0_gpkg: str = ADM0_GPKG):
    """Filter GeoDataFrame by World Bank region codes.

    Args:
        gdf: GeoDataFrame to filter (admin1 level)
        region_codes: List of WB region codes (e.g., ['AFR', 'SAR']) or full names
        adm0_gpkg: Path to Admin0 GPKG containing WB_REGION column

    Returns:
        Filtered GeoDataFrame containing only features in specified regions
    """
    if not region_codes:
        return gdf

    # Normalize region codes (accept both codes and full names)
    normalized_codes = []
    for code in region_codes:
        code_upper = code.upper()
        if code_upper in WB_REGIONS:
            normalized_codes.append(code_upper)
        else:
            # Try matching by name
            for wb_code, wb_name in WB_REGIONS.items():
                if code.lower() in wb_name.lower():
                    normalized_codes.append(wb_code)
                    break
            else:
                print(f"Warning: Unknown region '{code}'. Valid: {', '.join(WB_REGIONS.keys())}")

    if not normalized_codes:
        print("No valid regions found. Returning full dataset.")
        return gdf

    print(f"Filtering by regions: {', '.join(normalized_codes)}")

    # Load Admin0 to get country-region mapping
    import geopandas as gpd
    if not Path(adm0_gpkg).exists():
        print(f"Warning: Admin0 file not found: {adm0_gpkg}")
        return gdf

    adm0 = gpd.read_file(adm0_gpkg)
    if WB_REGION_COL not in adm0.columns:
        print(f"Warning: {WB_REGION_COL} column not found in Admin0")
        return gdf

    # Get countries in the specified regions
    countries_in_regions = adm0[adm0[WB_REGION_COL].isin(normalized_codes)]

    # Find matching column between admin0 and admin1
    # Try common columns: NAM_0, SOVEREIGN, ISO_A3
    match_cols = [("NAM_0", "NAM_0"), ("SOVEREIGN", "SOVEREIGN"), ("ISO_A3", "ISO_A3")]

    for adm0_col, adm1_col in match_cols:
        if adm0_col in countries_in_regions.columns and adm1_col in gdf.columns:
            country_values = countries_in_regions[adm0_col].dropna().unique()
            filtered = gdf[gdf[adm1_col].isin(country_values)]
            print(f"  Matched {len(filtered)} features via {adm1_col}")
            return filtered

    print("Warning: Could not find matching column between Admin0 and Admin1")
    return gdf


def auto_merge(gdf, csv_path: str = DEFAULT_CSV):
    """Auto-merge CSV data to GeoDataFrame if possible."""
    import pandas as pd
    
    if not Path(csv_path).exists():
        return gdf
    
    df = pd.read_csv(csv_path)
    gpkg_cols = set(gdf.columns)
    csv_cols = set(df.columns)
    common = gpkg_cols & csv_cols
    
    if not common:
        return gdf
    
    # Pick merge key
    merge_key = next((p for p in MERGE_KEYS if p in common), sorted(common)[0])
    
    new_cols = [c for c in df.columns if c not in gpkg_cols or c == merge_key]
    df_subset = df[new_cols].drop_duplicates(subset=[merge_key])
    
    merged = gdf.merge(df_subset, on=merge_key, how="left")
    added_cols = [c for c in merged.columns if c not in gpkg_cols]
    if added_cols:
        print(f"Auto-merged {len(added_cols)} columns from CSV (key: {merge_key})")
        print(f"  Added columns: {', '.join(added_cols)}")
    return merged


def get_label_col(gdf):
    """Get best available label column."""
    return next((c for c in LABEL_COLUMNS if c in gdf.columns), None)


def load_cities(gdf, top_n: int = None):
    """Load top N cities within the filtered regions, sorted by population.

    Args:
        gdf: GeoDataFrame of filtered regions (cities must fall within these)
        top_n: Number of cities to return

    Returns None if cities data not available.
    """
    if not Path(CITIES_SHP).exists():
        print(f"Cities file not found: {CITIES_SHP}")
        return None

    import geopandas as gpd

    print(f"Loading cities: {CITIES_SHP}")
    cities = gpd.read_file(CITIES_SHP)
    print(f"  Loaded {len(cities)} cities")

    # Convert cities to match the reference CRS
    if cities.crs != gdf.crs:
        cities = cities.to_crs(gdf.crs)

    # Spatial join: keep only cities within the filtered regions
    cities_in_regions = gpd.sjoin(cities, gdf[["geometry"]], predicate="within", how="inner")
    print(f"  Cities in filtered regions: {len(cities_in_regions)}")

    if len(cities_in_regions) == 0:
        return None

    # Sort by population and take top N
    if CITIES_POP_COL in cities_in_regions.columns:
        cities_in_regions = cities_in_regions.sort_values(CITIES_POP_COL, ascending=False)

    n = top_n or CITIES_TOP_N
    cities_in_regions = cities_in_regions.head(n)
    print(f"  Selected top {len(cities_in_regions)} cities")

    # List city names
    if CITIES_NAME_COL in cities_in_regions.columns:
        for _, row in cities_in_regions.iterrows():
            name = row[CITIES_NAME_COL]
            pop = row.get(CITIES_POP_COL, 0) if CITIES_POP_COL in cities_in_regions.columns else 0
            print(f"    - {name}: pop={pop:,.0f}")

    return cities_in_regions


def plot_static(gpkg_path: str = DEFAULT_GPKG, column: str = None, filter_terms: list = None, regions: list = None, output_dir: str = DEFAULT_OUTPUT_DIR, show: bool = False, cities: bool = True, labels: bool = False):
    """Plot static map using matplotlib with auto-coloring and labels.

    Args:
        gpkg_path: Path to GPKG file
        column: Column to color by
        filter_terms: List of country/region names to filter (fuzzy match)
        regions: List of WB region codes to filter (e.g., ['AFR', 'SAR'])
        output_dir: Output directory
        show: Show plot interactively
        cities: Include cities layer
        labels: Show admin1 region labels
    """
    import geopandas as gpd
    import matplotlib.pyplot as plt
    import pandas as pd
    import numpy as np

    print(f"Loading: {gpkg_path}")
    gdf = gpd.read_file(gpkg_path)
    print(f"Loaded {len(gdf)} features, {len(gdf.columns)} columns")
    print(f"  Columns: {', '.join(gdf.columns)}")

    gdf = auto_merge(gdf)

    # Apply filters (region filter first, then fuzzy filter)
    if regions:
        gdf = filter_by_region(gdf, regions)
    if filter_terms:
        gdf = fuzzy_filter(gdf, filter_terms)

    # Determine filter label for output file naming
    filter_label = None
    if regions:
        filter_label = "_".join(regions)
    if filter_terms:
        filter_label = (filter_label + "_" if filter_label else "") + "_".join(filter_terms)
    
    fig, ax = plt.subplots(1, 1, figsize=(14, 10))
    
    # Get bounds with padding for context
    b = gdf.total_bounds  # [minx, miny, maxx, maxy]
    pad_x = (b[2] - b[0]) * 0.1
    pad_y = (b[3] - b[1]) * 0.1
    
    # Define view bounds (what will actually be visible)
    view_bounds = (b[0] - pad_x, b[1] - pad_y, b[2] + pad_x, b[3] + pad_y)
    
    # Load Admin 0 for context
    adm0 = None
    if Path(ADM0_GPKG).exists():
        adm0 = gpd.read_file(ADM0_GPKG)
        if filter_label:
            # Select countries intersecting view bounds (cx preserves GeoDataFrame)
            adm0_view = adm0.cx[view_bounds[0]:view_bounds[2], view_bounds[1]:view_bounds[3]]
            # Draw neighboring countries in light gray first
            if len(adm0_view) > 0:
                adm0_view.plot(ax=ax, facecolor="#f0f0f0", edgecolor="#cccccc", linewidth=0.5)
    
    # Draw regions with unique colors per region (not by country)
    gdf = gdf.reset_index(drop=True)
    gdf["_color_idx"] = range(len(gdf))  # unique color per region
    gdf.plot(column="_color_idx", ax=ax, legend=False, cmap="Set3", 
             edgecolor="white", linewidth=0.5, alpha=0.75)
    
    # Draw Admin 0 borders ON TOP
    if adm0 is not None and filter_label:
        adm0_view = adm0.cx[view_bounds[0]:view_bounds[2], view_bounds[1]:view_bounds[3]]
        if len(adm0_view) > 0:
            adm0_view.plot(ax=ax, facecolor="none", edgecolor="black", linewidth=1.5)
            
            # Add country names only if centroid is within view
            name_col = next((c for c in ["NAM_0", "ADM0_NAME", "WB_NAME", "ADMIN"] if c in adm0_view.columns), None)
            if name_col:
                for idx, row in adm0_view.iterrows():
                    centroid = row.geometry.centroid
                    # Only label if centroid is inside view bounds
                    if (view_bounds[0] < centroid.x < view_bounds[2] and 
                        view_bounds[1] < centroid.y < view_bounds[3]):
                        name = str(row[name_col]) if pd.notna(row[name_col]) else ""
                        if name:
                            ax.text(centroid.x, centroid.y, name, fontsize=9, ha='center', va='center',
                                    fontweight='bold', color='#333333', alpha=0.6)
    
    # Add region labels when requested (with fallback columns)
    REGION_LABEL_COLS = ["NAM_1_GAUL", "NAM_1_STAT", "P_NAME_1", "NAM_1", "ADM1_NAME"]
    if labels and len(gdf) < 100:
        try:
            from adjustText import adjust_text
            texts = []
            for idx, row in gdf.iterrows():
                # Try columns in order until finding a non-empty value
                label = ""
                for col in REGION_LABEL_COLS:
                    if col in gdf.columns and pd.notna(row[col]) and str(row[col]).strip():
                        label = str(row[col])
                        break
                if label:
                    centroid = row.geometry.centroid
                    texts.append(ax.text(centroid.x, centroid.y, label, fontsize=7, ha='center',
                                       bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.7, edgecolor='none')))
            if texts:
                adjust_text(texts, ax=ax, arrowprops=dict(arrowstyle='-', color='gray', lw=0.5))
                print(f"Added {len(texts)} region labels")
        except ImportError:
            print("Note: Install 'adjustText' for non-overlapping labels")
    
    # Set focused view (same as view_bounds)
    if filter_label:
        ax.set_xlim(view_bounds[0], view_bounds[2])
        ax.set_ylim(view_bounds[1], view_bounds[3])
    
    # Add cities layer (on top of everything)
    cities_added = False
    if cities:
        city_data = load_cities(gdf)
        if city_data is not None and len(city_data) > 0:
            cities_added = True
            # Size and color based on population if available
            if CITIES_POP_COL in city_data.columns:
                pop = city_data[CITIES_POP_COL].fillna(0)
                # Normalize sizes: min 30, max 200
                pop_norm = (pop - pop.min()) / (pop.max() - pop.min() + 1)
                sizes = 30 + pop_norm * 170
                # Plot with color gradient based on population
                scatter = ax.scatter(
                    city_data.geometry.x, city_data.geometry.y,
                    c=pop, cmap='YlOrRd', s=sizes,
                    edgecolor='black', linewidth=0.5, zorder=10
                )
                # Add colorbar legend for population (in millions)
                cbar = plt.colorbar(scatter, ax=ax, shrink=0.5, pad=0.02)
                cbar.set_label('Population (M)', fontsize=9)
                cbar.ax.tick_params(labelsize=8)
                cbar.formatter = plt.FuncFormatter(lambda x, _: f'{x/1e6:.1f}')
                cbar.update_ticks()
            else:
                ax.scatter(
                    city_data.geometry.x, city_data.geometry.y,
                    c='red', s=50, edgecolor='black', linewidth=0.5, zorder=10
                )
            # Add city names
            for _, row in city_data.iterrows():
                name = row[CITIES_NAME_COL] if CITIES_NAME_COL in city_data.columns else ""
                if pd.notna(name) and str(name).strip():
                    ax.annotate(str(name), (row.geometry.x, row.geometry.y),
                                fontsize=7, ha='left', va='bottom',
                                xytext=(3, 3), textcoords='offset points',
                                fontweight='bold', color='#222222',
                                bbox=dict(boxstyle='round,pad=0.15', facecolor='white', alpha=0.8, edgecolor='none'))
    
    # Add legend
    from matplotlib.lines import Line2D
    legend_elements = []
    legend_elements.append(Line2D([0], [0], marker='s', color='w', markerfacecolor='#d4edda', 
                                   markersize=10, markeredgecolor='white', label='Admin regions'))
    if cities_added:
        legend_elements.append(Line2D([0], [0], marker='o', color='w', markerfacecolor='#ff6b35',
                                       markersize=8, markeredgecolor='black', label='Cities'))
    if filter_label:
        legend_elements.append(Line2D([0], [0], color='black', linewidth=1.5, label='Country borders'))
    ax.legend(handles=legend_elements, loc='lower left', fontsize=8, framealpha=0.9)

    ax.set_title(filter_label.replace("_", " ").title() if filter_label else "World Map", fontsize=14, fontweight='bold')
    ax.set_axis_off()
    plt.tight_layout()

    out = make_output_path(output_dir, gpkg_path, filter_label, column or "", "png")
    fig.savefig(out, dpi=150, bbox_inches="tight")
    print(f"Saved: {out}")
    
    if show:
        plt.show()
    plt.close()


def plot_interactive(gpkg_path: str = DEFAULT_GPKG, column: str = None, filter_terms: list = None, regions: list = None, output_dir: str = DEFAULT_OUTPUT_DIR, show: bool = True, cities: bool = True):
    """Plot interactive map using Folium.

    Args:
        gpkg_path: Path to GPKG file
        column: Column to color by
        filter_terms: List of country/region names to filter (fuzzy match)
        regions: List of WB region codes to filter (e.g., ['AFR', 'SAR'])
        output_dir: Output directory
        show: Open in browser
        cities: Include cities layer
    """
    import geopandas as gpd
    import folium
    import pandas as pd

    print(f"Loading: {gpkg_path}")
    gdf = gpd.read_file(gpkg_path)
    print(f"Loaded {len(gdf)} features, {len(gdf.columns)} columns")
    print(f"  Columns: {', '.join(gdf.columns)}")

    gdf = auto_merge(gdf)

    # Apply filters (region filter first, then fuzzy filter)
    if regions:
        gdf = filter_by_region(gdf, regions)
    if filter_terms:
        gdf = fuzzy_filter(gdf, filter_terms)

    # Determine filter label for output file naming
    filter_label = None
    if regions:
        filter_label = "_".join(regions)
    if filter_terms:
        filter_label = (filter_label + "_" if filter_label else "") + "_".join(filter_terms)
    
    if gdf.crs and gdf.crs.to_epsg() != 4326:
        gdf = gdf.to_crs(epsg=4326)
    
    bounds = gdf.total_bounds
    center = [(bounds[1] + bounds[3]) / 2, (bounds[0] + bounds[2]) / 2]
    pad = (bounds[2] - bounds[0]) * 0.1  # 10% padding
    
    m = folium.Map(location=center, zoom_start=5)
    
    # Add Admin 0 borders (no fill) - clipped to view when filtered
    if Path(ADM0_GPKG).exists() and filter_label:
        adm0 = gpd.read_file(ADM0_GPKG)
        if adm0.crs and adm0.crs.to_epsg() != 4326:
            adm0 = adm0.to_crs(epsg=4326)
        # Clip to view bounds
        adm0 = adm0.cx[bounds[0]-pad:bounds[2]+pad, bounds[1]-pad:bounds[3]+pad]
        if len(adm0) > 0:
            border_style = lambda x: {"fillColor": "none", "color": "black", "weight": 2, "fillOpacity": 0}
            folium.GeoJson(adm0.__geo_interface__, style_function=border_style, name="Country borders").add_to(m)
    
    # Add data layer (regions) - on top of borders with unique colors
    tooltip_cols = [c for c in ["NAM_0", "NAM_1", "NAM_1_STAT", "P_NAME_1", "SOVEREIGN", "WB_REGION"] if c in gdf.columns]
    
    # Assign colors to each region
    import matplotlib.colors as mcolors
    cmap = mcolors.ListedColormap(plt.colormaps["Set3"].colors[:len(gdf)])
    gdf = gdf.reset_index(drop=True)
    gdf["_color"] = [mcolors.to_hex(cmap(i)) for i in range(len(gdf))]
    
    def region_style(feature):
        return {
            "fillColor": feature["properties"].get("_color", "steelblue"),
            "color": "white",
            "weight": 1,
            "fillOpacity": 0.7
        }
    
    folium.GeoJson(
        gdf.__geo_interface__,
        style_function=region_style,
        tooltip=folium.GeoJsonTooltip(fields=tooltip_cols) if tooltip_cols else None,
        name="Regions"
    ).add_to(m)
    
    # Add cities layer (on top) in a FeatureGroup for layer control
    if cities:
        city_data = load_cities(gdf)
        if city_data is not None and len(city_data) > 0:
            cities_group = folium.FeatureGroup(name="Cities")

            # Get population range for sizing
            if CITIES_POP_COL in city_data.columns:
                pop = city_data[CITIES_POP_COL].fillna(0)
                pop_min, pop_max = pop.min(), pop.max()
                pop_range = pop_max - pop_min + 1

            for _, row in city_data.iterrows():
                name = row[CITIES_NAME_COL] if CITIES_NAME_COL in city_data.columns else ""
                name = str(name) if pd.notna(name) else ""
                population = row[CITIES_POP_COL] if CITIES_POP_COL in city_data.columns else 0
                population = population if pd.notna(population) else 0

                # Size based on population (radius 4-15)
                if CITIES_POP_COL in city_data.columns and pop_range > 0:
                    pop_norm = (population - pop_min) / pop_range
                    radius = 4 + pop_norm * 11
                    # Color gradient: yellow to red based on population
                    r = int(255)
                    g = int(255 * (1 - pop_norm * 0.8))
                    b = int(50 * (1 - pop_norm))
                    color = f'#{r:02x}{g:02x}{b:02x}'
                else:
                    radius = 6
                    color = 'red'

                tooltip_text = f"<b>{name}</b><br>Pop: {population:,.0f}" if population else name
                folium.CircleMarker(
                    location=[row.geometry.y, row.geometry.x],
                    radius=radius,
                    color='black',
                    weight=1,
                    fill=True,
                    fill_color=color,
                    fill_opacity=0.9,
                    tooltip=tooltip_text
                ).add_to(cities_group)

            cities_group.add_to(m)

    # Add layer control to toggle layers on/off
    folium.LayerControl().add_to(m)

    m.fit_bounds([[bounds[1]-pad, bounds[0]-pad], [bounds[3]+pad, bounds[2]+pad]])
    
    out = make_output_path(output_dir, gpkg_path, filter_label, "interactive", "html")
    m.save(str(out))
    print(f"Saved: {out}")
    
    if show:
        import webbrowser
        webbrowser.open(str(out.absolute()))
    
    return m


def merge_csv_to_gpkg(gpkg_path: str = DEFAULT_GPKG, csv_path: str = DEFAULT_CSV, output_path: str = None):
    """Merge CSV data to GPKG by auto-detecting common columns and save."""
    import geopandas as gpd
    import pandas as pd
    
    print(f"Loading GPKG: {gpkg_path}")
    gdf = gpd.read_file(gpkg_path)
    print(f"  → {len(gdf)} features, {len(gdf.columns)} columns")
    
    print(f"Loading CSV: {csv_path}")
    df = pd.read_csv(csv_path)
    print(f"  → {len(df)} rows, {len(df.columns)} columns")
    
    gpkg_cols = set(gdf.columns)
    csv_cols = set(df.columns)
    common = gpkg_cols & csv_cols
    
    if not common:
        print("Error: No common columns found for merge.")
        print(f"  GPKG columns: {sorted(gpkg_cols)}")
        print(f"  CSV columns: {sorted(csv_cols)}")
        return None
    
    print(f"Common columns: {sorted(common)}")
    
    merge_key = next((p for p in MERGE_KEYS if p in common), sorted(common)[0])
    print(f"Using merge key: {merge_key}")
    
    new_cols = [c for c in df.columns if c not in gpkg_cols or c == merge_key]
    df_subset = df[new_cols].drop_duplicates(subset=[merge_key])
    
    merged = gdf.merge(df_subset, on=merge_key, how="left")
    matched = merged[merge_key].isin(df_subset[merge_key]).sum()
    print(f"Merged: {matched}/{len(gdf)} features matched")
    
    if output_path:
        Path(output_path).parent.mkdir(parents=True, exist_ok=True)
        merged.to_file(output_path, driver="GPKG")
        print(f"Saved: {output_path}")
    else:
        out = make_output_path(DEFAULT_OUTPUT_DIR, gpkg_path, suffix="merged", ext="gpkg")
        merged.to_file(str(out), driver="GPKG")
        print(f"Saved: {out}")
    
    return merged


def simplify_geometries(gpkg_path: str, tolerance: float = 0.01, output_path: str = None):
    """Simplify geometries to reduce file size.
    
    Args:
        gpkg_path: Path to input GPKG file
        tolerance: Simplification tolerance in degrees (~0.01 = ~1km)
        output_path: Output path (default: adds _simple suffix)
    
    Returns:
        Path to simplified file
    """
    import geopandas as gpd
    
    input_path = Path(gpkg_path)
    if not input_path.exists():
        print(f"Error: File not found: {gpkg_path}")
        return None
    
    # Default output path: same name with _simple suffix
    if output_path is None:
        output_path = input_path.parent / f"{input_path.stem}_simple.gpkg"
    else:
        output_path = Path(output_path)
    
    # Get original file size
    original_size = input_path.stat().st_size / (1024 * 1024)  # MB
    
    print(f"Loading: {gpkg_path}")
    gdf = gpd.read_file(gpkg_path)
    print(f"  Features: {len(gdf)}, Columns: {len(gdf.columns)}")
    print(f"  Original size: {original_size:.2f} MB")
    
    # Simplify geometries
    print(f"Simplifying with tolerance={tolerance}...")
    gdf["geometry"] = gdf.geometry.simplify(tolerance=tolerance, preserve_topology=True)
    
    # Save
    output_path.parent.mkdir(parents=True, exist_ok=True)
    gdf.to_file(str(output_path), driver="GPKG")
    
    # Get new file size
    new_size = output_path.stat().st_size / (1024 * 1024)  # MB
    reduction = (1 - new_size / original_size) * 100
    
    print(f"Saved: {output_path}")
    print(f"  New size: {new_size:.2f} MB ({reduction:.1f}% reduction)")
    
    return output_path


def simplify_all_boundaries(tolerance: float = 0.01):
    """Simplify both Admin 0 and Admin 1 boundary files."""
    results = []
    for gpkg in [ADM0_GPKG, ADM1_GPKG]:
        if Path(gpkg).exists():
            result = simplify_geometries(gpkg, tolerance)
            if result:
                results.append(result)
    return results


def main():
    parser = argparse.ArgumentParser(description="Geo utilities for maps and data merging")
    subparsers = parser.add_subparsers(dest="command", help="Commands")
    
    # Common filter arguments (added to static and interactive)
    region_help = f"Filter by WB region: {', '.join(WB_REGIONS.keys())} (can specify multiple)"
    filter_help = "Filter by country/region name (can specify multiple, e.g., -f kenya tanzania)"

    # Static map
    p_static = subparsers.add_parser("static", help="Plot static map")
    p_static.add_argument("gpkg", nargs="?", default=DEFAULT_GPKG, help="Path to GPKG file")
    p_static.add_argument("--column", "-c", help="Column to color by")
    p_static.add_argument("--region", "-r", dest="regions", nargs="+", help=region_help)
    p_static.add_argument("--filter", "-f", dest="filter_terms", nargs="+", help=filter_help)
    p_static.add_argument("--output-dir", "-o", default=DEFAULT_OUTPUT_DIR, help="Output directory")
    p_static.add_argument("--show", "-s", action="store_true", help="Show plot (default: save only)")
    p_static.add_argument("--no-cities", action="store_true", help="Don't show cities layer")
    p_static.add_argument("--labels", "-l", action="store_true", help="Show admin1 region labels")

    # Interactive map
    p_inter = subparsers.add_parser("interactive", help="Plot interactive map")
    p_inter.add_argument("gpkg", nargs="?", default=DEFAULT_GPKG, help="Path to GPKG file")
    p_inter.add_argument("--column", "-c", help="Column to color by")
    p_inter.add_argument("--region", "-r", dest="regions", nargs="+", help=region_help)
    p_inter.add_argument("--filter", "-f", dest="filter_terms", nargs="+", help=filter_help)
    p_inter.add_argument("--output-dir", "-o", default=DEFAULT_OUTPUT_DIR, help="Output directory")
    p_inter.add_argument("--no-show", "-n", action="store_true", help="Don't open browser")
    p_inter.add_argument("--no-cities", action="store_true", help="Don't show cities layer")
    
    # Merge
    p_merge = subparsers.add_parser("merge", help="Merge CSV data to GPKG")
    p_merge.add_argument("gpkg", nargs="?", default=DEFAULT_GPKG, help="Path to GPKG file")
    p_merge.add_argument("csv", nargs="?", default=DEFAULT_CSV, help="Path to CSV file")
    p_merge.add_argument("--output", "-o", help="Output GPKG path")
    
    # Simplify
    p_simplify = subparsers.add_parser("simplify", help="Simplify geometries to reduce file size")
    p_simplify.add_argument("gpkg", nargs="?", help="Path to GPKG file (default: simplify all boundaries)")
    p_simplify.add_argument("--tolerance", "-t", type=float, default=0.01, help="Simplification tolerance (~0.01 = ~1km)")
    p_simplify.add_argument("--output", "-o", help="Output GPKG path")
    
    args = parser.parse_args()

    if args.command == "static":
        plot_static(args.gpkg, args.column, args.filter_terms, args.regions, args.output_dir, args.show, cities=not args.no_cities, labels=args.labels)
    elif args.command == "interactive":
        plot_interactive(args.gpkg, args.column, args.filter_terms, args.regions, args.output_dir, not args.no_show, cities=not args.no_cities)
    elif args.command == "merge":
        merge_csv_to_gpkg(args.gpkg, args.csv, args.output)
    elif args.command == "simplify":
        if args.gpkg:
            simplify_geometries(args.gpkg, args.tolerance, args.output)
        else:
            simplify_all_boundaries(args.tolerance)
    else:
        parser.print_help()


if __name__ == "__main__":
    main()