Get Panoramax Sequences in a Bounding Box as a GeoPackage File

panoramax_sequences.py

import requests
import mercantile
import mapbox_vector_tile
from shapely.geometry import LineString
import geopandas as gpd

def transform(x_local, y_local, west, south, east, north, extent=4096):
    """
    Transforme les coordonnées locales d'une tuile MVT en coordonnées géographiques (WGS84).

    Args:
        x_local (float): Coordonnée locale x (en pixels ou unités arbitraires).
        y_local (float): Coordonnée locale y (en pixels ou unités arbitraires).
        west (float): Longitude du coin ouest de la tuile.
        south (float): Latitude du coin sud de la tuile.
        east (float): Longitude du coin est de la tuile.
        north (float): Latitude du coin nord de la tuile.
        extent (int, optional): Taille de la tuile en unités locales. Par défaut, 4096.

    Returns:
        tuple: (longitude, latitude) en WGS84.
    """
    x_geo = west + (x_local / extent) * (east - west)
    y_geo = south + (y_local / extent) * (north - south)
    return x_geo, y_geo

def get_sequence(left, bottom, right, top, output_file="sequences.gpkg"):
    """
    Télécharge une tuile MVT depuis Panoramax pour une étendue géographique donnée,
    extrait les géométries de type LineString, les transforme en coordonnées géographiques,
    et les enregistre dans un fichier GeoPackage.

    Args:
        left (float): Longitude minimale de l'étendue géographique (WGS84).
        bottom (float): Latitude minimale de l'étendue géographique (WGS84).
        right (float): Longitude maximale de l'étendue géographique (WGS84).
        top (float): Latitude maximale de l'étendue géographique (WGS84).
        output_file (str, optional): Chemin du fichier de sortie. Par défaut, "sequences.gpkg".

    Returns:
        None: Le résultat est enregistré dans un fichier GeoPackage.
    """
    # 1. Déterminer la tuile MVT couvrant l'étendue géographique
    tile = mercantile.bounding_tile(left, bottom, right, top)

    # 2. Récupérer l'étendue géographique de la tuile (en WGS84)
    bounds = mercantile.bounds(tile)
    east = bounds.east
    west = bounds.west
    north = bounds.north
    south = bounds.south

    # 3. Construire l'URL pour télécharger la tuile MVT
    x, y, z = tile.x, tile.y, tile.z
    url = f'https://panoramax.openstreetmap.fr/api/map/{z}/{x}/{y}.mvt'
    r = requests.get(url, allow_redirects=True)

    # 4. Télécharger le contenu binaire de la tuile MVT
    response = requests.get(url, allow_redirects=True)
    response.raise_for_status()  # Vérifier que la requête a réussi

    # 5. Décoder le contenu binaire en mémoire
    decoded_data = mapbox_vector_tile.decode(response.content)

    # 6. Extraire les features de la couche 'sequences'
    extent = decoded_data['sequences']['extent']
    feats = decoded_data['sequences']['features']

    # 7. Transformer les coordonnées locales en coordonnées géographiques
    lines = []
    datas = []
    for i, feat in enumerate(feats):
        
        coords = feat['geometry']['coordinates']
        data = feat['properties']
        
        coords2 = []
        # Cas 1 : Les coordonnées sont des paires de valeurs (ex: [x, y])
        if isinstance(coords[0][0], int):
            for coord in coords:
                coord2 = transform(coord[0], coord[1], west, south, east, north, extent)
                coords2.append(coord2)
        # Cas 2 : Les coordonnées sont imbriquées dans des listes (ex: [[x, y], ...])
        else:
            for groupe in coords:
                for coord in groupe:
                    coord2 = transform(coord[0], coord[1], west, south, east, north, extent)
                    coords2.append(coord2)

        # Créer une LineString à partir des coordonnées transformées
        line = LineString(coords2)
        
        lines.append(line)
        datas.append(data)

    print(f"{len(lines)} lignes créées.")

    # 8. Créer un GeoDataFrame avec les LineString
    lines_df = gpd.GeoDataFrame(data=datas, geometry=lines, crs="EPSG:4326")

    # 9. Exporter les résultats dans un fichier GeoPackage
    lines_df.to_file(output_file, driver="GPKG")
    
# Exemple d'appel de la fonction
get_sequence(5.74144, 45.15984, 5.74985, 45.16798, "sequences.gpkg")