osm-labo/counting_osm_objects/get_poly.py

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

"""
Script pour récupérer le polygone d'une commune française à partir de son code INSEE.

Ce script:
1. Demande un code INSEE
2. Interroge l'API Overpass Turbo pour obtenir les limites administratives
3. Extrait le polygone de la commune
4. Sauvegarde le polygone dans un fichier

Usage:
    python get_poly.py [code_insee]

Si le code INSEE n'est pas fourni en argument, le script le demandera interactivement.
"""

import sys
import os
import json
import urllib.parse
import urllib.request
import argparse


def get_insee_code():
    """
    Récupère le code INSEE soit depuis les arguments de ligne de commande,
    soit en demandant à l'utilisateur.

    Returns:
        str: Le code INSEE de la commune
    """
    parser = argparse.ArgumentParser(
        description="Récupère le polygone d'une commune à partir de son code INSEE"
    )
    parser.add_argument("insee", nargs="?", help="Code INSEE de la commune")
    args = parser.parse_args()

    if args.insee:
        return args.insee

    # Si le code INSEE n'est pas fourni en argument, le demander
    return input("Entrez le code INSEE de la commune: ")


def query_overpass_api(insee_code):
    """
    Interroge l'API Overpass pour obtenir les limites administratives d'une commune.

    Args:
        insee_code (str): Le code INSEE de la commune

    Returns:
        dict: Les données GeoJSON de la commune
    """
    print(f"Récupération des limites administratives pour la commune {insee_code}...")

    # Construire la requête Overpass QL pour obtenir la relation administrative
    query = f"""
    [out:json][timeout:60];
    (
      relation["boundary"="administrative"]["admin_level"="8"]["ref:INSEE"="{insee_code}"];
      way(r);
      node(w);
    );
    out geom;
    """

    # Encoder la requête pour l'URL
    encoded_query = urllib.parse.quote(query)

    # Construire l'URL de l'API Overpass
    url = f"https://overpass-api.de/api/interpreter?data={encoded_query}"

    try:
        # Envoyer la requête à l'API
        print("Envoi de la requête à Overpass API...")
        with urllib.request.urlopen(url) as response:
            data = json.loads(response.read().decode("utf-8"))

        # Afficher des informations sur la réponse (version réduite pour production)
        print(
            f"Réponse reçue de l'API Overpass. Nombre d'éléments: {len(data.get('elements', []))}"
        )

        return data
    except Exception as e:
        print(f"Erreur lors de la requête à l'API Overpass: {e}")
        raise RuntimeError(f"Erreur lors de la requête à l'API Overpass: {e}")


def extract_polygon(data):
    """
    Extrait le polygone des données GeoJSON.

    Args:
        data (dict): Les données GeoJSON de la commune

    Returns:
        list: Liste des coordonnées du polygone
    """
    print("Extraction du polygone des données...")

    # Vérifier si des éléments ont été trouvés
    if not data.get("elements"):
        print("Aucune limite administrative trouvée pour ce code INSEE.")
        raise ValueError("Aucune limite administrative trouvée pour ce code INSEE.")

    try:
        # Collecter tous les nœuds (points) avec leurs coordonnées
        nodes = {}
        for element in data["elements"]:
            if element["type"] == "node":
                nodes[element["id"]] = (element["lon"], element["lat"])

        # Trouver les ways qui forment le contour de la commune
        ways = []
        for element in data["elements"]:
            if element["type"] == "way":
                ways.append(element)

        # Si aucun way n'est trouvé, essayer d'extraire directement les coordonnées des nœuds
        if not ways and nodes:
            print("Aucun way trouvé. Utilisation directe des nœuds...")
            polygon = list(nodes.values())
            return polygon

        # Trouver la relation administrative
        relation = None
        for element in data["elements"]:
            if (
                element["type"] == "relation"
                and element.get("tags", {}).get("boundary") == "administrative"
            ):
                relation = element
                break

        if not relation:
            print("Aucune relation administrative trouvée.")
            # Si nous avons des ways, nous pouvons essayer de les utiliser directement
            if ways:
                print("Tentative d'utilisation directe des ways...")
                # Prendre le premier way comme contour
                way = ways[0]
                polygon = []
                for node_id in way.get("nodes", []):
                    if node_id in nodes:
                        polygon.append(nodes[node_id])
                return polygon
            raise ValueError(
                "Impossible de trouver une relation administrative ou des ways"
            )

        # Extraire les ways qui forment le contour extérieur de la relation
        outer_ways = []
        for member in relation.get("members", []):
            if member.get("role") == "outer" and member.get("type") == "way":
                # Trouver le way correspondant
                for way in ways:
                    if way["id"] == member["ref"]:
                        outer_ways.append(way)
                        break

        # Si aucun way extérieur n'est trouvé, utiliser tous les ways
        if not outer_ways:
            print("Aucun way extérieur trouvé. Utilisation de tous les ways...")
            outer_ways = ways

        # Construire le polygone à partir des ways extérieurs
        polygon = []
        for way in outer_ways:
            for node_id in way.get("nodes", []):
                if node_id in nodes:
                    polygon.append(nodes[node_id])

        if not polygon:
            raise ValueError("Impossible d'extraire le polygone de la relation")

        print(f"Polygone extrait avec {len(polygon)} points.")
        return polygon
    except Exception as e:
        print(f"Erreur lors de l'extraction du polygone: {e}")
        raise RuntimeError(f"Erreur lors de l'extraction du polygone: {e}")


def save_polygon_to_file(polygon, insee_code):
    """
    Sauvegarde le polygone dans un fichier.

    Args:
        polygon (list): Liste des coordonnées du polygone
        insee_code (str): Le code INSEE de la commune

    Returns:
        str: Le chemin du fichier créé

    Raises:
        ValueError: Si le polygone est vide ou invalide
        IOError: Si une erreur survient lors de l'écriture du fichier
    """
    if not polygon:
        raise ValueError("Le polygone est vide")

    try:
        # Créer le répertoire de sortie s'il n'existe pas
        output_dir = os.path.join(
            os.path.dirname(os.path.abspath(__file__)), "polygons"
        )
        os.makedirs(output_dir, exist_ok=True)

        # Définir le nom du fichier de sortie
        output_file = os.path.join(output_dir, f"commune_{insee_code}.poly")

        print(f"Sauvegarde du polygone dans le fichier {output_file}...")

        # Écrire le polygone dans le fichier au format .poly (format utilisé par Osmosis)
        with open(output_file, "w") as f:
            f.write(f"commune_{insee_code}\n")
            f.write("1\n")  # Numéro de section

            # Écrire les coordonnées
            for i, (lon, lat) in enumerate(polygon):
                f.write(f"    {lon:.7f}    {lat:.7f}\n")

            # Fermer le polygone en répétant le premier point
            if len(polygon) > 1 and polygon[0] != polygon[-1]:
                lon, lat = polygon[0]
                f.write(f"    {lon:.7f}    {lat:.7f}\n")

            f.write("END\n")
            f.write("END\n")

        print(f"Polygone sauvegardé avec succès dans {output_file}")
        return output_file
    except IOError as e:
        print(f"Erreur lors de l'écriture du fichier: {e}")
        raise  # Re-raise the IOError
    except Exception as e:
        print(f"Erreur inattendue lors de la sauvegarde du polygone: {e}")
        raise RuntimeError(f"Erreur inattendue lors de la sauvegarde du polygone: {e}")


def main():
    """
    Fonction principale du script.
    """
    try:
        # Récupérer le code INSEE
        insee_code = get_insee_code()

        # Vérifier que le code INSEE est valide (format numérique ou alphanumérique pour les DOM-TOM)
        if not insee_code:
            raise ValueError("Le code INSEE ne peut pas être vide")

        if not insee_code.isalnum() or len(insee_code) not in [5, 3]:
            raise ValueError(
                "Code INSEE invalide. Il doit être composé de 5 chiffres (ou 3 pour certains territoires)."
            )

        # Interroger l'API Overpass
        data = query_overpass_api(insee_code)

        # Extraire le polygone
        polygon = extract_polygon(data)

        # Sauvegarder le polygone dans un fichier
        output_file = save_polygon_to_file(polygon, insee_code)

        print(
            f"Terminé. Le polygone de la commune {insee_code} a été sauvegardé dans {output_file}"
        )
        return 0  # Succès
    except ValueError as e:
        print(f"Erreur de validation: {e}")
        return 1  # Erreur
    except KeyboardInterrupt:
        print("\nOpération annulée par l'utilisateur.")
        return 1  # Erreur
    except Exception as e:
        print(f"Erreur inattendue: {e}")
        return 1  # Erreur


if __name__ == "__main__":
    sys.exit(main())