Algorithme pour trouver un ensemble minimal de codes postaux pour couvrir (presque) complètement un pays

__author__ = "Uli Köhler"
__copyright__ = "Copyright 2025, Uli Köhler"
__license__ = "CC0 1.0 Universal (CC0 1.0) Public Domain Dedication"

import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point
from shapely.ops import unary_union
import matplotlib.pyplot as plt
from tqdm.auto import tqdm
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.feature import ShapelyFeature
from sklearn.neighbors import KDTree
import cartopy.io.shapereader as shpreader

# Chargeurs PLZ et pays intégrés (autonomes)
def load_plz_data(plz_file_path: str) -> pd.DataFrame:
    """Charger les données tabulaires de codes postaux depuis un chemin de fichier (parquet ou CSV)."""
    if plz_file_path.endswith('.parquet') or plz_file_path.endswith('.parq'):
        return pd.read_parquet(plz_file_path)
    if plz_file_path.endswith('.csv') or plz_file_path.endswith('.txt'):
        return pd.read_csv(plz_file_path)
    # Essayer les lecteurs courants comme solution de secours
    try:
        return pd.read_parquet(plz_file_path)
    except Exception:
        try:
            return pd.read_csv(plz_file_path)
        except Exception as exc:
            raise RuntimeError(f'Impossible de lire les données PLZ depuis {plz_file_path} : {exc}') from exc


def create_plz_points(plz_df: pd.DataFrame) -> gpd.GeoDataFrame:
    """Créer un GeoDataFrame de points PLZ à partir d'un DataFrame.
    Attend longitude/latitude dans des colonnes nommées parmi ('lon','lng','longitude') et ('lat','latitude').
    Garantit également l'existence d'une colonne 'PLZ' en copiant les colonnes courantes de codes postaux ou en utilisant l'index comme solution de secours.
    """
    df = plz_df.copy()
    # recherche tolérante de colonnes pour les coordonnées
    lon_cols = [c for c in df.columns if c.lower() in ('lon', 'lng', 'longitude')]
    lat_cols = [c for c in df.columns if c.lower() in ('lat', 'latitude')]
    if not lon_cols or not lat_cols:
        raise ValueError('Le DataFrame PLZ doit contenir des colonnes lon et lat (par ex. lon, lat)')
    lon_col = lon_cols[0]
    lat_col = lat_cols[0]
    df = df.dropna(subset=[lon_col, lat_col]).copy()
    geometry = [Point(xy) for xy in zip(df[lon_col].astype(float), df[lat_col].astype(float))]
    gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")
    # Conserver des noms de colonnes conventionnels pour le code en aval
    if 'lon' not in gdf.columns:
        gdf['lon'] = gdf[lon_col].astype(float)
    if 'lat' not in gdf.columns:
        gdf['lat'] = gdf[lat_col].astype(float)

    # Garantir l'existence d'une colonne PLZ (code postal). Essayer d'abord les noms courants, sinon utiliser l'index comme chaîne
    plz_cols = [c for c in gdf.columns if c.lower() in ('plz', 'postal_code', 'postal', 'postleitzahl')]
    if plz_cols:
        gdf['PLZ'] = gdf[plz_cols[0]].astype(str)
    else:
        # S'il y a une colonne numérique nommée 'zip' ou 'postcode', la choisir
        zip_cols = [c for c in gdf.columns if c.lower() in ('zip', 'postcode')]
        if zip_cols:
            gdf['PLZ'] = gdf[zip_cols[0]].astype(str)
        else:
            # solution de secours : créer PLZ à partir de l'index pour avoir un identifiant stable
            gdf = gdf.reset_index(drop=False)
            gdf['PLZ'] = gdf['index'].astype(str)
            gdf.drop(columns=['index'], inplace=True)

    return gdf


def load_germany_from_natural_earth() -> object:
    """Lire les pays Natural Earth 10m et retourner la géométrie de l'Allemagne (shapely)"""
    # Utiliser le chemin shapereader de cartopy et geopandas pour charger le shapefile 10m admin_0_countries
    shp_path = shpreader.natural_earth(resolution='10m', category='cultural', name='admin_0_countries')
    countries = gpd.read_file(shp_path)
    # essayer les noms de colonnes ISO courants pour correspondre à l'Allemagne
    iso_cols = [c for c in ('ISO_A3', 'ISO3', 'adm0_a3', 'ADM0_A3', 'iso_a3') if c in countries.columns]
    if iso_cols:
        germany = countries[countries[iso_cols[0]] == 'DEU']
    else:
        # solution de secours : essayer la correspondance par NAME ou SOVEREIGNT
        name_cols = [c for c in ('NAME_EN', 'NAME', 'SOVEREIGNT', 'ADMIN') if c in countries.columns]
        if name_cols:
            germany = countries[countries[name_cols[0]].str.contains('Germany', case=False, na=False)]
        else:
            raise RuntimeError('Impossible de déterminer la colonne ISO ou nom dans les pays Natural Earth')
    if germany.empty:
        raise RuntimeError('Géométrie de l\'Allemagne introuvable dans les pays Natural Earth')
    return germany.iloc[0].geometry


def filter_plz_in_germany(plz_gdf, germany_geometry):
    """Filtrer les points PLZ qui sont en Allemagne"""
    plz_gdf = plz_gdf.to_crs("EPSG:3857")  # Web Mercator pour les calculs de distance
    germany_geometry = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:3857")[0]

    # Créer une zone tampon autour de l'Allemagne pour inclure les PLZ près de la frontière
    buffered_germany = germany_geometry.buffer(50000)  # zone tampon de 50 km

    # Filtrer les PLZ dans l'Allemagne avec zone tampon
    plz_in_germany = plz_gdf[plz_gdf.intersects(buffered_germany)]
    return plz_in_germany


def calculate_coverage(selected_plz, germany_geometry, radius_km=30):
    """Calculer le pourcentage de l'Allemagne couvert par les PLZ sélectionnés"""
    radius_m = radius_km * 1000
    coverage = unary_union(selected_plz.geometry.buffer(radius_m))
    germany_geometry = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:3857")[0]
    covered_area = germany_geometry.intersection(coverage).area
    total_area = germany_geometry.area
    return covered_area / total_area * 100


def _generate_grid_points_in_polygon(polygon_3857, step_m: float) -> np.ndarray:
    """Générer une grille régulière de points (x, y) dans le polygone donné (EPSG:3857)."""
    xmin, ymin, xmax, ymax = polygon_3857.bounds
    xs = np.arange(xmin, xmax + step_m, step_m)
    ys = np.arange(ymin, ymax + step_m, step_m)
    points = []
    # Une boucle simple suffit pour ~20k-60k points
    for x in xs:
        for y in ys:
            p = Point(x, y)
            if polygon_3857.contains(p):
                points.append((x, y))
    if not points:
        return np.empty((0, 2), dtype=float)
    return np.array(points, dtype=float)


def select_plz_max_coverage(
    plz_gdf: gpd.GeoDataFrame,
    germany_geometry,
    radius_km: float = 30,
    coverage_target: float = 0.99,
    sample_size: int = 300,
    grid_step_km: float = 5.0,
    exclusion_factor: float = 0.75,
    max_steps: int = 10000,
    patience: int = 500,
    validate_every: int = 50,
    random_state: int = 42,
    verbose: bool = True
) -> gpd.GeoDataFrame:
    """
    Algorithme unique pour sélectionner des points PLZ maximisant la zone couverte (approx).

    Utilise une grille régulière de points à l'intérieur de l'Allemagne (EPSG:3857) comme proxy de zone.
    Glouton randomisé : à chaque étape, échantillonner un sous-ensemble des candidats restants et choisir
    celui qui couvre le plus de points de grille actuellement non couverts dans le rayon.

    Améliorations :
    - Vérifications exactes périodiques de la couverture (unary_union) pour éviter un arrêt prématuré.
    - Rayon d'exclusion pour supprimer les sélections groupées : après avoir choisi un centre,
      éliminer les PLZ restants dans (exclusion_factor * radius_km).

    Retourne : GeoDataFrame (EPSG:3857) des points PLZ sélectionnés.
    """
    rng = np.random.default_rng(random_state)

    # Garantir un CRS projeté pour distance/superficie et préparer la géométrie
    work_gdf = plz_gdf.to_crs("EPSG:3857").copy()
    germany_3857 = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:3857")[0]

    radius_m = float(radius_km) * 1000.0
    exclusion_m = float(exclusion_factor) * radius_m
    grid_step_m = float(grid_step_km) * 1000.0

    # 1) Construire une grille de points à l'intérieur de l'Allemagne comme proxy de couverture
    if verbose:
        print(f"Construction d'une grille de {grid_step_km:.1f} km sur l'Allemagne pour l'approximation de couverture...")
    grid_xy = _generate_grid_points_in_polygon(germany_3857, step_m=grid_step_m)
    if grid_xy.shape[0] == 0:
        raise RuntimeError("Échec de la génération de points de grille à l'intérieur du polygone de l'Allemagne.")
    if verbose:
        print(f"Points de grille : {grid_xy.shape[0]:,}")

    # 2) KDTree sur la grille pour des requêtes de rayon rapides
    tree_grid = KDTree(grid_xy, leaf_size=40, metric='euclidean')

    # 3) Précalculer les indices de voisinage de grille pour chaque candidat PLZ
    plz_xy = np.column_stack([work_gdf.geometry.x.values, work_gdf.geometry.y.values])
    if verbose:
        print("Précalcul des voisinages de couverture candidat-vers-grille...")
    neighbor_lists = tree_grid.query_radius(plz_xy, r=radius_m, return_distance=False)

    # 3b) KDTree sur les candidats PLZ pour le filtrage d'exclusion
    tree_plz = KDTree(plz_xy, leaf_size=40, metric='euclidean')

    # 4) Boucle de sélection gloutonne randomisée sur le masque de grille non couvert
    uncovered = np.ones(grid_xy.shape[0], dtype=bool)
    selected_mask = np.zeros(len(work_gdf), dtype=bool)
    active_mask = np.ones(len(work_gdf), dtype=bool)  # candidats encore disponibles
    no_improve_steps = 0

    pbar = tqdm(total=max_steps, disable=not verbose, desc="Sélection des centres (couverture-max)")
    for step in range(max_steps):
        pbar.update(1)

        # Vérification exacte périodique de la couverture
        if step % validate_every == 0 and selected_mask.any():
            selected_tmp = work_gdf.loc[selected_mask]
            exact_cov = calculate_coverage(selected_tmp, germany_geometry, radius_km) / 100.0
            if step % (validate_every * 2) == 0:
                pbar.set_postfix_str(f"exact~{exact_cov*100:.1f}%")
            if exact_cov >= coverage_target:
                break

        # Fraction de couverture approximative sur la grille
        covered_frac = 1.0 - (uncovered.sum() / uncovered.size)
        if covered_frac >= coverage_target and selected_mask.any():
            selected_tmp = work_gdf.loc[selected_mask]
            exact_cov = calculate_coverage(selected_tmp, germany_geometry, radius_km) / 100.0
            if exact_cov >= coverage_target:
                break

        remaining_idx = np.flatnonzero(active_mask)
        if remaining_idx.size == 0:
            break

        # Échantillonner des candidats depuis l'ensemble actif
        k = int(min(sample_size, remaining_idx.size))
        sample = rng.choice(remaining_idx, size=k, replace=False) if k > 0 else []

        # Évaluer le gain marginal sur la grille
        best_gain = 0
        best_idx = None
        for idx in sample:
            neigh = neighbor_lists[idx]
            if neigh.size == 0:
                continue
            gain = int(uncovered[neigh].sum())
            if gain > best_gain:
                best_gain = gain
                best_idx = idx

        if best_idx is None or best_gain == 0:
            no_improve_steps += 1
            if no_improve_steps >= patience:
                break
            else:
                continue

        # Accepter le meilleur candidat
        selected_mask[best_idx] = True
        # Marquer les points de grille comme couverts
        uncovered[neighbor_lists[best_idx]] = False
        # Exclure les PLZ à proximité pour réduire le regroupement
        if exclusion_m > 0:
            near = tree_plz.query_radius(plz_xy[[best_idx]], r=exclusion_m, return_distance=False)[0]
            active_mask[near] = False
        # Retirer également l'index choisi
        active_mask[best_idx] = False
        no_improve_steps = 0

    pbar.close()

    selected = work_gdf.loc[selected_mask].copy()
    selected = selected.set_crs("EPSG:3857")
    return selected


def visualize_coverage(selected_plz, germany_geometry, radius_km=30):
    """Visualiser la couverture à l'aide de Cartopy"""
    # Convertir en WGS84 pour la visualisation
    selected_plz = selected_plz.to_crs("EPSG:4326")
    germany_geometry = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:4326")[0]

    # Créer la figure avec la projection Cartopy
    fig = plt.figure(figsize=(12, 12))
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.EqualEarth())

    # Ajouter les entités cartographiques
    ax.add_feature(cfeature.LAND, facecolor='lightgray')
    ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
    ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
    ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5)

    # Créer un ShapelyFeature pour l'Allemagne
    germany_feature = ShapelyFeature([germany_geometry], ccrs.PlateCarree(), facecolor='none', edgecolor='red', linewidth=1)
    ax.add_feature(germany_feature)

    # Créer et tracer les cercles (TOUS les PLZ maintenant)
    radius_m = radius_km * 1000
    buffers_3857 = gpd.GeoSeries(selected_plz.geometry).to_crs("EPSG:3857").buffer(radius_m)
    buffers_4326 = gpd.GeoSeries(buffers_3857, crs="EPSG:3857").to_crs("EPSG:4326")

    # Tracer chaque zone tampon individuellement pour que les chevauchements s'affichent comme des transparences superposées
    for geom in buffers_4326:
        ax.add_geometries([geom], crs=ccrs.PlateCarree(), facecolor='blue', edgecolor='none', alpha=0.2, zorder=2)

    # Tracer les points PLZ sélectionnés par-dessus
    ax.scatter(selected_plz.geometry.x, selected_plz.geometry.y, color='red', s=5, transform=ccrs.PlateCarree(), label='PLZ sélectionnés', zorder=3)

    # Définir l'étendue sur l'Allemagne avec un peu de marge
    ax.set_extent([4, 16, 47, 56], crs=ccrs.PlateCarree())

    # Ajouter des lignes de grille
    gl = ax.gridlines(draw_labels=True, linestyle='--')
    gl.top_labels = False
    gl.right_labels = False

    plt.title(f"Couverture PLZ de l'Allemagne avec {len(selected_plz)} centres (rayon de {radius_km} km)")
    plt.legend(loc='upper right')
    plt.show()


def run_analysis(plz_file_path, radius_km=30, coverage_target=0.99):
    """Fonction principale pour l'analyse utilisant un seul algorithme de couverture maximale"""
    # Charger la géométrie de l'Allemagne depuis Natural Earth
    print("Chargement de la géométrie de l'Allemagne depuis Natural Earth...")
    germany_geometry = load_germany_from_natural_earth()

    # Charger et traiter les données PLZ
    print("Chargement des données PLZ...")
    plz_df = load_plz_data(plz_file_path)
    print(f"{len(plz_df)} codes postaux allemands chargés")

    plz_gdf = create_plz_points(plz_df)

    print("Filtrage des PLZ en Allemagne...")
    plz_in_germany = filter_plz_in_germany(plz_gdf, germany_geometry)
    print(f"{len(plz_in_germany)} PLZ trouvés dans ou près de l'Allemagne")

    # Sélectionner les centres PLZ avec l'algorithme unique
    print("Sélection des centres PLZ avec l'algorithme de couverture maximale...")
    selected_plz = select_plz_max_coverage(
        plz_in_germany, germany_geometry, radius_km=radius_km, coverage_target=coverage_target, sample_size=300, grid_step_km=2.0, exclusion_factor=0.75, max_steps=10000, patience=500, validate_every=50, random_state=42, verbose=True
    )

    # Calculer la couverture (union exacte, unique)
    coverage_percent = calculate_coverage(selected_plz, germany_geometry, radius_km)
    print(f"\n{len(selected_plz)} centres PLZ sélectionnés couvrant {coverage_percent:.1f}% de l'Allemagne")

    # Visualiser la couverture
    visualize_coverage(selected_plz, germany_geometry=germany_geometry, radius_km=radius_km)

    # Retourner les résultats défensivement : n'inclure que les colonnes qui existent
    sel = selected_plz.to_crs("EPSG:4326")
    desired = ['PLZ', 'lat', 'lon', 'geometry']
    available = [c for c in desired if c in sel.columns or c == 'geometry']
    # Garantir la présence de la géométrie
    if 'geometry' not in available:
        available.append('geometry')

    return {
        'selected_plz': sel[available],
        'coverage_percent': coverage_percent,
        'total_plz_used': len(selected_plz)
    }

results = run_analysis(
    plz_file_path="postleitzahlen.parquet",
    radius_km=35,
    coverage_target=0.995
)

# Cellule : trouver les 50 PLZ sélectionnés les plus proches d'un PLZ donné (exemple : 63110)
plz_code = '10176'

# Garantir que les résultats sélectionnés sont disponibles ; sinon, exécuter l'analyse une fois (garde légère)
if 'results' not in globals():
    print('`results` introuvable dans l\'espace de noms du notebook — exécution de `run_analysis` (cela peut prendre du temps)')
    results = run_analysis(plz_file_path="postleitzahlen.parquet", radius_km=35, coverage_target=0.995)

selected = results['selected_plz']

# Charger la table PLZ complète pour localiser les coordonnées du PLZ de référence
plz_df = load_plz_data("postleitzahlen.parquet")
plz_gdf = create_plz_points(plz_df)

# Essayer une correspondance exacte sur PLZ, sinon une correspondance par préfixe
ref = plz_gdf[plz_gdf['PLZ'].astype(str) == str(plz_code)]
if ref.empty:
    ref = plz_gdf[plz_gdf['PLZ'].astype(str).str.startswith(str(plz_code))]

if ref.empty:
    raise ValueError(f"PLZ de référence {plz_code} introuvable dans le jeu de données PLZ")

ref_point = ref.iloc[0].geometry

# Reprojeter les deux dans un CRS projeté pour des distances métriques (EPSG:3857)
selected_3857 = selected.to_crs("EPSG:3857").copy()
ref_3857 = gpd.GeoSeries([ref_point], crs="EPSG:4326").to_crs("EPSG:3857")[0]

# Calculer les distances (mètres)
selected_3857['distance_m'] = selected_3857.geometry.distance(ref_3857)

# Si le PLZ de référence lui-même est dans l'ensemble sélectionné, l'exclure pour obtenir les autres centres les plus proches
if 'PLZ' in selected_3857.columns:
    selected_for_sort = selected_3857[selected_3857['PLZ'].astype(str) != str(plz_code)]
else:
    selected_for_sort = selected_3857

# Prendre les 50 plus proches
nearest = selected_for_sort.nsmallest(50, 'distance_m').copy()

# Reconvertir en coordonnées géographiques et préparer un DataFrame convivial
nearest_geo = nearest.to_crs("EPSG:4326")
nearest_geo['lat'] = nearest_geo.geometry.y
nearest_geo['lon'] = nearest_geo.geometry.x
nearest_geo['distance_km'] = (nearest_geo['distance_m'] / 1000.0).round(3)

# Choisir les colonnes à afficher : privilégier les champs humains s'ils existent
prefer_columns = ['PLZ', 'name', 'place', 'Ortsteil', 'stadtteil', 'place_name', 'lat', 'lon', 'distance_km']
available = [c for c in prefer_columns if c in nearest_geo.columns]
# si aucune des colonnes de nom préférées n'existe, afficher toutes les colonnes non-géométriques
if not available:
    available = [c for c in nearest_geo.columns if c != 'geometry']

nearest_50_df = nearest_geo[available].reset_index(drop=True)

print(f"{len(nearest_50_df)} PLZ sélectionnés les plus proches de {plz_code} (top {min(50, len(nearest_50_df))}) :")
nearest_50_df.head(50)  # afficher le DataFrame