batiments.py

import os
import argparse
import requests
import geopandas as gpd
import pandas as pd
import folium
from pyproj import Transformer
from shapely.geometry import Point
from folium import LayerControl


def fetch_buildings(x: float, y: float, radius: float):
    """Fetch building features from Geoportail WFS within a bounding box."""
    minx, miny = x - radius, y - radius
    maxx, maxy = x + radius, y + radius

    url = "https://data.geopf.fr/wfs/ows"
    params = {
        "service": "WFS",
        "version": "1.1.0",
        "request": "GetFeature",
        "typename": "BDTOPO_V3:batiment",
        "outputFormat": "application/json",
        "srsname": "EPSG:3857",
        "bbox": f"{minx},{miny},{maxx},{maxy},EPSG:3857",
    }

    response = requests.get(url, params=params)
    response.raise_for_status()
    return response.json()


def process_buildings(geojson: dict) -> gpd.GeoDataFrame:
    """Convert GeoJSON to GeoDataFrame and clean attributes."""
    gdf = gpd.GeoDataFrame.from_features(geojson["features"], crs="EPSG:3857")

    if "date_d_apparition" in gdf.columns:
        gdf["annee_construction"] = pd.to_datetime(
            gdf["date_d_apparition"], errors="coerce", format="%Y-%m-%dZ"
        ).dt.year.astype("Int64")

    return gdf


def get_building_gdf(lat: float, lon: float, radius: float) -> gpd.GeoDataFrame:
    """Fetch and process buildings around a point."""
    transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
    x, y = transformer.transform(lon, lat)

    buildings_data = fetch_buildings(x, y, radius)
    gdf = process_buildings(buildings_data)

    gdf = gdf.to_crs("EPSG:2154")

    gdf_polygon = gpd.GeoDataFrame(geometry=[Point(lon, lat)], crs="EPSG:4326").to_crs(
        "EPSG:2154"
    )
    gdf_polygon["geometry"] = gdf_polygon["geometry"].buffer(radius)

    len_before = len(gdf)
    gdf = gpd.sjoin(gdf, gdf_polygon, how="inner", predicate="within")
    gdf = gdf.reset_index(drop=True).drop(columns=["index_right"])
    len_after = len(gdf)
    print(f"Removed {len_before - len_after} buildings outside the polygon.")

    return gdf


def get_building_map(gdf: gpd.GeoDataFrame, column: str):
    """Generate a folium map from a GeoDataFrame colored by a given column."""
    tiles = folium.TileLayer(
        tiles="https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}",
        attr="Google",
        name="Google Satellite",
        overlay=False,
        control=True,
        max_zoom=25,
    )

    m = gdf.explore(
        column=column,
        cmap="jet",
        tiles="OpenStreetMap",
        popup=False,
        name="GeoDataframe",
    )
    tiles.add_to(m)
    LayerControl().add_to(m)
    return m


def main():
    parser = argparse.ArgumentParser(
        description="Fetch building data from Geoportail and generate maps."
    )
    parser.add_argument(
        "--lat", type=float, required=True, help="Latitude in EPSG:4326"
    )
    parser.add_argument(
        "--lon", type=float, required=True, help="Longitude in EPSG:4326"
    )
    parser.add_argument("--radius", type=float, required=True, help="Radius in meters")
    parser.add_argument("--output", type=str, required=True, help="Output folder path")

    args = parser.parse_args()
    lat, lon, radius, output = args.lat, args.lon, args.radius, args.output

    os.makedirs(output, exist_ok=True)

    print("Fetching building data...")
    gdf = get_building_gdf(lat, lon, radius)

    gpkg_path = os.path.join(output, "batiments.gpkg")
    print(f"Saving GeoPackage to {gpkg_path}")
    gdf.to_file(gpkg_path, driver="GPKG")

    if "hauteur" in gdf.columns:
        print("Generating height map...")
        m1 = get_building_map(gdf, column="hauteur")
        m1.save(os.path.join(output, "buildings_hauteur.html"))

    else:
        print("⚠️ No 'hauteur' column found in data, skipping height map.")

    if "annee_construction" in gdf.columns:
        print("Generating construction year map...")
        m2 = get_building_map(gdf, column="annee_construction")
        m2.save(os.path.join(output, "buildings_annee.html"))
    else:
        print("⚠️ No 'annee_construction' column found in data, skipping year map.")

    print("✅ All done!")


if __name__ == "__main__":
    main()