Minimal Tile Coverage

minimal_tile_coverage.py

import numpy as np
from scipy.optimize import linprog
import pandas as pd
import geopandas as gpd
from pyproj.crs.crs import CRS
from shapely.ops import(
    unary_union,
    polygonize,
)

EPSILON = 1e-7

def minimal_tile_coverage(
    tiles: list[gpd.GeoDataFrame],
    roi: gpd.GeoDataFrame,
) -> list[gpd.GeoDataFrame]:
    """Minimal Tile Coverage

    Given a set of tiles (GeoDataFrames with a single feature each)
    and a ROI (GeoDataFrame with a single feature), finds out
    the minimal set of tiles completely covering the ROI by solving the
    underlying linear programming problem.

    Args:
        tiles (list[gpd.GeoDataFrame]): Tiles used for covering the ROI.
        roi (gpd.GeoDataFrame): The ROI to cover.
    
    Returns:
        gpd.GeoDataFrame: The list with the fewest number of tiles
            completely covering the ROI.
    """
    # Reproject to common epsg
    tiles = [tile.to_crs(roi.crs) for tile in tiles]
    
    # Create the set of overlapping parts
    boundaries = unary_union([tile.geometry[0].exterior for tile in tiles])
    polygons = [
        p for p in polygonize(boundaries)
        if roi.intersects(p).values[0]
    ]
    tiles_gdf = pd.concat(tiles)

    # Create constraints to solve
    A_ub = -np.array(
        [
            tiles_gdf.covers(p.buffer(-EPSILON)).values
            for p in polygons
        ],
        dtype=int,
    )
    b_ub = -np.ones(A_ub.shape[0:1])
    bounds = (0, 1)
    c = np.ones((len(tiles),))

    # Solve
    res = linprog(
        c = c,
        A_ub = A_ub,
        b_ub = b_ub,
        bounds = bounds,
        integrality = 1,
    )

    # Return selected tiles
    return [
        t for t,p in zip(tiles, res.x) if p
    ]