points2prints ¶

def convert_bd_topo_call(
    input_path: Path,
    output_path: Path,
    input_output: InputOutput,
    verbose: Verbose,
):
    """Convert a BD TOPO file to Parquet format.

    Parameters
    ----------
    input_path : Path
        Path to the input BD TOPO file containing all the layers including the buildings layer (e.g., a .gpkg file).
    output_path : Path
        Path to the output parquet file.
    input_output: InputOutput
        The handler for input and output file issues.
    verbose: Verbose
        The verbosity level for logging.
    """
    with LoggingContext(verbose=verbose):
        convert_bd_topo_implementation(
            input_path=input_path,
            output_path=output_path,
            input_output=input_output,
        )

def convert_bd_topo_implementation(
    input_path: Path,
    output_path: Path,
    input_output: InputOutput,
):
    """Convert a BD TOPO file to Parquet format.

    Parameters
    ----------
    input_path : Path
        Path to the input BD TOPO file containing all the layers including the buildings layer (e.g., a .gpkg file).
    output_path : Path
        Path to the output parquet file.
    input_output: InputOutput
        The handler for input and output file issues.
    """
    message_prefix = "Converting BD TOPO to Parquet"
    input_output.handle_input(
        message_prefix=message_prefix,
        input_files=[input_path],
    )
    output_action = input_output.handle_output(
        message_prefix=message_prefix,
        behaviour=OutputBehaviour.ALL_OR_NOTHING,
        output_files=[[output_path]],
    )
    if output_action == OutputActionEnum.SKIP:
        return

    # Create a temporary directory to store intermediate files
    logging.info(f"Creating temporary directory for processing {input_path.name}.")
    with TemporaryDirectory() as temp_dir:
        temp_dir_path = Path(temp_dir)

        # Extract the buildings from the BD TOPO file
        logging.info(f"Extracting buildings from {input_path.name}.")
        sql_query = f"SELECT ST_ForcePolygonCCW({INITIAL_GEOMETRY_COLUMN_NAME}) AS {FINAL_GEOMETRY_COLUMN_NAME}, * FROM batiment"
        extracted_file_path = temp_dir_path / "extracted.gpkg"
        extraction_command = [
            "gdal",
            "vector",
            "sql",
            "-i",
            str(input_path),
            "-o",
            str(extracted_file_path),
            "--output-layer",
            "buildings",
            "--sql",
            sql_query,
        ]

        return_code = run_command_with_tqdm_logging(extraction_command, display=True)
        if return_code != 0:
            logging.error(f"Failed to extract buildings from {input_path.name}.")
        else:
            logging.info(f"Successfully extracted buildings from {input_path.name}.")

        # Add a column with the integer part of the building ID
        logging.info(f"Processing extracted data from {input_path.name}.")
        db_path = temp_dir_path / "duckdb.duckdb"
        processed_file_path = temp_dir_path / "processed.parquet"
        with DuckDBConnectionManager(db_path) as con:
            con.create_schema(SCHEMA_NAME)

            read_query = f"""
                CREATE OR REPLACE TABLE {SCHEMA_NAME}.{TABLE_NAME} AS (
                    SELECT CAST(array_slice(cleabs, 9, -1) AS INT64) AS cleabs_int, *
                    FROM '{extracted_file_path}'
                );
            """

            con.execute(read_query)

            con.export_parquet(
                schema_name=SCHEMA_NAME,
                table_name=TABLE_NAME,
                geom_col_name=FINAL_GEOMETRY_COLUMN_NAME,
                output_file=processed_file_path,
            )

        # Format the file properly with gpio
        logging.info(
            f"Formatting the processed data from {input_path.name} to Parquet."
        )
        format_command = [
            "gpio",
            "convert",
            str(processed_file_path),
            str(output_path),
        ]

        return_code = run_command_with_tqdm_logging(format_command, display=True)
        if return_code != 0:
            logging.error(
                f"Failed to format the data from {input_path.name} to Parquet."
            )
        else:
            logging.info(
                f"Successfully formatted the data from {input_path.name} to Parquet."
            )

async def collect_existing_tiles(
    tiles: List[Box2154],
    concurrency: int = DEFAULT_CONCURRENCY,
) -> Dict[str, Optional[str]]:
    """Collect the best matching WFS tile for each requested tile box.

    Parameters
    ----------
    tiles
        Requested tiles in EPSG:2154.
    concurrency
        Maximum HTTP concurrency for the client.

    Returns
    -------
    dict[str, str | None]
        Mapping from requested tile name to the matched download URL.
    """
    # The WFS endpoint enforces a strict per-second rate limit; querying
    # sequentially avoids silent partial/empty responses.
    semaphore = asyncio.Semaphore(1)
    timeout = httpx.Timeout(timeout=60.0, connect=30.0)
    limits = httpx.Limits(
        max_connections=concurrency, max_keepalive_connections=concurrency
    )
    name_to_url: Dict[str, Optional[str]] = {
        _tile_box_to_name(tile_box): None for tile_box in tiles
    }

    async with httpx.AsyncClient(
        timeout=timeout, limits=limits, headers=STAC_HEADERS
    ) as client:
        with tqdm(total=len(tiles), desc="Checking tiles", unit="tile") as pbar:
            for i, tile_box in enumerate(tiles):
                logging.debug(f"Checking existence of tile: {tile_box}")
                tile_results = await _fetch_tiles_by_bbox(client, semaphore, tile_box)
                requested_tile_name = _tile_box_to_name(tile_box)

                matched_url: Optional[str] = None
                matched_tile_name: Optional[str] = None
                best_overlap = -1
                for tile_name, tile_url in tile_results:
                    candidate_box = _nw_coordinates_to_tile_box(tile_name)
                    if candidate_box is None:
                        continue

                    overlap = _intersection_area(tile_box, candidate_box) / 1000000
                    # We only accept a match if the overlap is greater than 90% of the requested tile area
                    if overlap > best_overlap and overlap > 0.9:
                        best_overlap = overlap
                        matched_url = tile_url
                        matched_tile_name = tile_name

                if matched_url is None and tile_results:
                    logging.debug(
                        f"No overlapping tile match for {requested_tile_name}; available tiles: {[tile_name for tile_name, _ in tile_results]}"
                    )
                elif (
                    matched_tile_name is not None
                    and matched_tile_name != requested_tile_name
                ):
                    logging.debug(
                        f"Requested {requested_tile_name} matched WFS tile {matched_tile_name} with overlap {best_overlap}"
                    )

                name_to_url[requested_tile_name] = matched_url

                pbar.update(1)
                if i < len(tiles) - 1:
                    await asyncio.sleep(WFS_MIN_QUERY_INTERVAL_SECONDS)

    return name_to_url

def download_lidar_hd_call(
    xmin: int,
    xmax: int,
    ymin: int,
    ymax: int,
    output_path_template: Path,
    input_output: InputOutput,
    verbose: Verbose,
    concurrency: Optional[int],
):
    """Download LiDAR HD data for a specified bounding box.

    Parameters
    ----------
    xmin : int
        Minimum X coordinate of the requested bounding box in EPSG:2154.
    xmax : int
        Maximum X coordinate of the requested bounding box in EPSG:2154.
    ymin : int
        Minimum Y coordinate of the requested bounding box in EPSG:2154.
    ymax : int
        Maximum Y coordinate of the requested bounding box in EPSG:2154.
    output_path_template : Path
        Path to save the downloaded files. The path can contain the values {xmin}, {ymin}, {xmax}, {ymax}, {file_name} which will be replaced with the corresponding values. The values also have their kilometre equivalents {xmin_km}, {ymin_km}, {xmax_km}, {ymax_km}.
    input_output: InputOutput
        The handler for input and output file issues.
    verbose: Verbose
        The verbosity level for logging.
    concurrency: Optional[int]
        Maximum number of concurrent download workers.
    """
    with LoggingContext(verbose=verbose):
        download_lidar_hd_data_implementation(
            xmin=xmin,
            xmax=xmax,
            ymin=ymin,
            ymax=ymax,
            output_path_template=output_path_template,
            input_output=input_output,
            concurrency=concurrency,
        )

async def download_lidar_hd_data_async(
    xmin: int,
    xmax: int,
    ymin: int,
    ymax: int,
    output_path_template: Path,
    input_output: InputOutput,
    concurrency: Optional[int],
) -> None:
    """Download LIDAR HD tiles covering a requested EPSG:2154 bounding box.

    Parameters
    ----------
    xmin : int
        Minimum X coordinate of the requested bounding box in EPSG:2154.
    xmax : int
        Maximum X coordinate of the requested bounding box in EPSG:2154.
    ymin : int
        Minimum Y coordinate of the requested bounding box in EPSG:2154.
    ymax : int
        Maximum Y coordinate of the requested bounding box in EPSG:2154.
    output_path_template : Path
        Output path template supporting the tile placeholders.
    input_output: InputOutput
        The handler for input and output file issues.
    concurrency : Optional[int]
        Maximum HTTP concurrency for tile discovery and downloading.
    """
    bbox = Box2154(Point2154(xmin, ymin), Point2154(xmax, ymax))
    tiles_boxes = bbox.get_tiles_boxes()

    logging.info(f"Generated {len(tiles_boxes)} tiles for the specified bounding box.")
    logging.debug("Tiles boxes:")
    for tile_box in tiles_boxes:
        logging.debug(
            f"  EPSG:2154({tile_box.p_min.x}, {tile_box.p_min.y}) -> EPSG:4326({tile_box.p_min.lon_4326}, {tile_box.p_min.lat_4326})"
        )

    # Discover WFS tiles intersecting the requested area.
    if concurrency is None:
        concurrency = DEFAULT_CONCURRENCY
    name_to_url_with_nones = await collect_existing_tiles(
        tiles=tiles_boxes, concurrency=concurrency
    )
    name_to_url = {
        name: url for name, url in name_to_url_with_nones.items() if url is not None
    }
    logging.info(
        f"Found {len(name_to_url)} existing tiles out of {len(name_to_url_with_nones)} candidates."
    )

    # Build the download path map, skipping files that already exist when requested.
    name_to_path: Dict[str, Path] = {}
    for tile_box in tiles_boxes:
        tile_name = _tile_box_to_name(tile_box)
        tile_url = name_to_url.get(tile_name)
        if tile_url is None:
            continue

        file_name = tile_url.split("/")[-1]

        output_path_template_str = str(output_path_template)
        output_path = output_path_template_str.format(
            xmin=tile_box.p_min.x,
            ymin=tile_box.p_min.y,
            xmax=tile_box.p_max.x,
            ymax=tile_box.p_max.y,
            xmin_km=tile_box.p_min.x // 1000,
            ymin_km=tile_box.p_min.y // 1000,
            xmax_km=tile_box.p_max.x // 1000,
            ymax_km=tile_box.p_max.y // 1000,
            file_name=file_name,
        )

        name_to_path[tile_name] = Path(output_path)

    # Handle input/output
    output_actions = input_output.get_output_actions(
        message_prefix="Downloaded LiDAR HD tiles",
        output_files=list(map(lambda path: [path], name_to_path.values())),
    )

    initial_keys = list(name_to_path.keys())
    for output_action, name in zip(output_actions, initial_keys):
        if output_action == OutputActionEnum.SKIP:
            name_to_url.pop(name)
            name_to_path.pop(name)
        elif output_action == OutputActionEnum.ERROR:
            raise Exception(f"Unexpected error when handling input/output for {name}")

    # Return if all tiles have already been downloaded
    if len(name_to_url) == 0:
        logging.info("No tiles to download after filtering. Exiting.")
        return

    # Download the tiles in parallel with retries and progress bars
    downloaded_count, failures, downloaded_files = await download_tiles(
        name_to_url,
        name_to_path,
        concurrency=concurrency,
    )
    logging.info(f"Downloaded {downloaded_count}/{len(name_to_url)} files.")

    # Log any download failures
    if failures:
        logging.error(f"{len(failures)} downloads failed:")
        for tile_name, error in failures:
            logging.error(f"  - {tile_name}: {error}")

    # Validate the downloaded files using PDAL and log results
    valid_count, invalid_files = validate_downloaded_files(downloaded_files)
    logging.info(f"Valid files: {valid_count}/{len(downloaded_files)}")
    if invalid_files:
        logging.error(f"{len(invalid_files)} invalid files:")
        for file_name, error in invalid_files:
            logging.error(f"  - {file_name}: {error}")

def download_lidar_hd_data_implementation(
    xmin: int,
    xmax: int,
    ymin: int,
    ymax: int,
    output_path_template: Path,
    input_output: InputOutput,
    concurrency: Optional[int],
):
    """Download LiDAR HD data for a specified bounding box.

    Parameters
    ----------
    xmin : int
        Minimum X coordinate of the requested bounding box in EPSG:2154.
    xmax : int
        Maximum X coordinate of the requested bounding box in EPSG:2154.
    ymin : int
        Minimum Y coordinate of the requested bounding box in EPSG:2154.
    ymax : int
        Maximum Y coordinate of the requested bounding box in EPSG:2154.
    output_path_template : Path
        Path to save the downloaded files. The path can contain the values {xmin}, {ymin}, {xmax}, {ymax}, {file_name} which will be replaced with the corresponding values. The values also have their kilometre equivalents {xmin_km}, {ymin_km}, {xmax_km}, {ymax_km}.
    input_output: InputOutput
        The handler for input and output file issues.
    concurrency: Optional[int]
        Maximum number of concurrent download workers.
    """
    asyncio.run(
        download_lidar_hd_data_async(
            xmin=xmin,
            xmax=xmax,
            ymin=ymin,
            ymax=ymax,
            output_path_template=output_path_template,
            input_output=input_output,
            concurrency=concurrency,
        )
    )

async def download_tiles(
    name_to_url: Dict[str, str],
    name_to_path: Dict[str, Path],
    concurrency: int = DEFAULT_CONCURRENCY,
) -> Tuple[int, List[Tuple[str, str]], List[Path]]:
    """Download the selected tiles in parallel.

    Parameters
    ----------
    name_to_url
        Mapping from tile name to download URL.
    name_to_path
        Mapping from tile name to output path.
    concurrency
        Maximum number of concurrent download workers.

    Returns
    -------
    tuple[int, list[tuple[str, str]], list[pathlib.Path]]
        Downloaded count, failures, and successfully written files.
    """
    timeout = httpx.Timeout(timeout=None, connect=30.0, write=60.0)
    limits = httpx.Limits(
        max_connections=concurrency, max_keepalive_connections=concurrency
    )
    failures: List[Tuple[str, str]] = []
    downloaded_count = 0
    downloaded_files: List[Path] = []
    progress_lock = asyncio.Lock()

    ordered_tiles = list(name_to_url.items())
    worker_count = min(concurrency, len(ordered_tiles)) if ordered_tiles else 0

    queue: asyncio.Queue[Optional[Tuple[str, str]]] = asyncio.Queue()
    for item in ordered_tiles:
        queue.put_nowait(item)
    for _ in range(worker_count):
        queue.put_nowait(None)

    async with httpx.AsyncClient(
        timeout=timeout, limits=limits, headers=STAC_HEADERS
    ) as client:
        results: List[Tuple[str, bool, Optional[str], Optional[Path]]] = []
        with tqdm(
            total=len(ordered_tiles),
            desc="Downloading tiles",
            unit="file",
            position=0,
            leave=True,
            dynamic_ncols=True,
        ) as overall_bar:
            workers = [
                asyncio.create_task(
                    _download_worker(
                        worker_id=worker_id,
                        queue=queue,
                        client=client,
                        name_to_path=name_to_path,
                        progress_lock=progress_lock,
                        overall_bar=overall_bar,
                        results=results,
                    )
                )
                for worker_id in range(worker_count)
            ]

            await queue.join()
            await asyncio.gather(*workers)

        for tile_name, success, error, file_path in results:
            if success:
                downloaded_count += 1
                if file_path is not None:
                    downloaded_files.append(file_path)
            else:
                failures.append((tile_name, error or "unknown error"))

    with suppress(Exception):
        tqdm.write("")

    return downloaded_count, failures, downloaded_files

def validate_downloaded_files(
    downloaded_files: List[Path],
) -> Tuple[int, List[Tuple[str, str]]]:
    """Validate downloaded LAZ files using PDAL.

    Parameters
    ----------
    downloaded_files
        Paths to files that were downloaded successfully.

    Returns
    -------
    tuple[int, list[tuple[str, str]]]
        Number of valid files and the invalid file/error pairs.
    """
    valid_count = 0
    invalid_files: List[Tuple[str, str]] = []

    if not downloaded_files:
        return valid_count, invalid_files

    with tqdm(
        total=len(downloaded_files), desc="Validating files", unit="file"
    ) as pbar:
        for file_path in downloaded_files:
            is_valid = False
            error_message = "unknown validation error"

            proc = subprocess.run(
                ["pdal", "info", "--summary", str(file_path)],
                capture_output=True,
                text=True,
            )
            if proc.returncode == 0:
                is_valid = True
            else:
                stderr_text = proc.stderr.strip()
                error_message = stderr_text if stderr_text else "pdal info failed"

            if is_valid:
                valid_count += 1
            else:
                invalid_files.append((file_path.name, error_message))

            pbar.update(1)

    return valid_count, invalid_files

def intersections_call(
    bd_topo_file: Path,
    output_edges_file: Path,
    output_intersections_file: Path,
    output_building_groups_file: Path,
    input_output: InputOutput,
    verbose: Verbose,
):
    """_summary_

    Parameters
    ----------
    bd_topo_file : Path
        Input BD TOPO file (Parquet) containing the building geometries.
    output_edges_file : Path
        Output Parquet file where the edges will be saved.
    output_intersections_file : Path
        Output Parquet file where the intersections will be saved.
    output_building_groups_file : Path
        Output Parquet file where the building groups will be saved.
    input_output: InputOutput
        The handler for input and output file issues.
    verbose: Verbose
        The verbosity level for logging.
    """
    with LoggingContext(verbose=verbose):
        intersections_implementation(
            bd_topo_file=bd_topo_file,
            output_edges_file=output_edges_file,
            output_intersections_file=output_intersections_file,
            output_building_groups_file=output_building_groups_file,
            input_output=input_output,
        )