#!/usr/bin/env python

# elements.py
# generating derived elements (street edge, block)
import os
import warnings

import geopandas as gpd
import libpysal
import numpy as np
import pandas as pd
import shapely
from packaging.version import Version
from scipy.spatial import Voronoi
from shapely.geometry.base import BaseGeometry
from shapely.ops import polygonize
from tqdm.auto import tqdm

from .utils import deprecated

__all__ = [
    "Tessellation",
    "Blocks",
    "get_network_id",
    "get_node_id",
    "enclosures",
    "get_network_ratio",
]

GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0")
GPD_GE_10 = Version(gpd.__version__) >= Version("1.0dev")


class Tessellation:
    """
    Generates tessellation. Three versions of tessellation can be created:

        1. Morphological tessellation around given buildings
            ``gdf`` within set ``limit``.
        2. Proximity bands around given street network ``gdf``
            within set ``limit``.
        3. Enclosed tessellation based on given buildings
            ``gdf`` within ``enclosures``.

    Pass either ``limit`` to create morphological tessellation or proximity bands or
    ``enclosures`` to create enclosed tessellation.

    See :cite:`fleischmann2020` for details of implementation of morphological
    tessellation and :cite:`araldi2019` for proximity bands.

    Tessellation requires data of relatively high level of precision
    and there are three particular patterns causing issues:

        1. Features will collapse into empty polygon - these
            do not have tessellation cell in the end.
        2. Features will split into MultiPolygons - in some cases,
            features with narrow links between parts split into two
            during 'shrinking'. In most cases that is not an issue
            and the resulting tessellation is correct anyway, but
            sometimes this results in a cell being a MultiPolygon,
            which is not correct.
        3. Overlapping features - features which overlap even
            after 'shrinking' cause invalid tessellation geometry.

    All three types can be tested prior :class:`momepy.Tessellation` using
    :class:`momepy.CheckTessellationInput`.

    Parameters
    ----------
    gdf : GeoDataFrame
       A GeoDataFrame containing building footprints or street network.
    unique_id : str
        The name of the column with the unique ID.
    limit : MultiPolygon or Polygon (default None)
        MultiPolygon or Polygon defining the study area limiting morphological
        tessellation or proximity bands (otherwise it could go to infinity).
    shrink : float (default 0.4)
        The distance for negative buffer to generate space between adjacent polygons
        (if geometry type of gdf is (Multi)Polygon).
    segment : float (default 0.5)
        The maximum distance between points after discretization.
    verbose : bool (default True)
        If ``True``, shows progress bars in loops and indication of steps.
    enclosures : GeoDataFrame (default None)
        The enclosures geometry, which can be generated
        using :func:`momepy.enclosures`.
    enclosure_id : str (default 'eID')
        The name of the ``enclosure_id`` containing unique identifer for each row in
        ``enclosures``. Applies only if ``enclosures`` are passed.
    threshold : float (default 0.05)
        The minimum threshold for a building to be considered within an enclosure.
        Threshold is a ratio of building area which needs to be within an enclosure to
        inlude it in the tessellation of that enclosure. Resolves sliver geometry
        issues. Applies only if ``enclosures`` are passed.
    use_dask : bool (default True)
        Use parallelised algorithm based on ``dask.dataframe``. Requires dask.
        Applies only if ``enclosures`` are passed.
    n_chunks : None
        The number of chunks to be used in parallelization. Ideal is one chunk per
        thread. Applies only if ``enclosures`` are passed. Default automatically
        uses ``n == dask.system.cpu_count``.

    Attributes
    ----------
    tessellation : GeoDataFrame
        A GeoDataFrame containing resulting tessellation.
        For enclosed tessellation, gdf contains three columns:

            - ``geometry``,
            - ``unique_id`` matching with parental building,
            - ``enclosure_id`` matching with enclosure integer index

    gdf : GeoDataFrame
        The original GeoDataFrame.
    id : Series
        A Series containing used unique ID.
    limit : MultiPolygon or Polygon
        MultiPolygon or Polygon defining the study area limiting morphological
        tessellation or proximity bands.
    shrink : float
        The distance for negative buffer to generate space between adjacent polygons.
    segment : float
        The maximum distance between points after discretization.
    collapsed : list
        A list of ``unique_id``s of collapsed features (if there are any).
        Applies only if ``limit`` is passed.
    multipolygons : list
        A list of ``unique_id``s of features causing MultiPolygons (if there are any).
        Applies only if ``limit`` is passed.

    Examples
    --------
    >>> tess = mm.Tessellation(
    ... buildings_df, 'uID', limit=mm.buffered_limit(buildings_df)
    ... )
    Inward offset...
    Generating input point array...
    Generating Voronoi diagram...
    Generating GeoDataFrame...
    Dissolving Voronoi polygons...
    >>> tess.tessellation.head()
        uID	geometry
    0	1	POLYGON ((1603586.677274485 6464344.667944215,...
    1	2	POLYGON ((1603048.399497852 6464176.180701573,...
    2	3	POLYGON ((1603071.342637536 6464158.863329805,...
    3	4	POLYGON ((1603055.834005827 6464093.614718676,...
    4	5	POLYGON ((1603106.417554705 6464130.215958447,...

    >>> enclosures = mm.enclosures(streets, admin_boundary, [railway, rivers])
    >>> encl_tess = mm.Tessellation(
    ... buildings_df, 'uID', enclosures=enclosures
    ... )
    >>> encl_tess.tessellation.head()
         uID                                           geometry  eID
    0  109.0  POLYGON ((1603369.789 6464340.661, 1603368.754...    0
    1  110.0  POLYGON ((1603368.754 6464340.097, 1603369.789...    0
    2  111.0  POLYGON ((1603458.666 6464332.614, 1603458.332...    0
    3  112.0  POLYGON ((1603462.235 6464285.609, 1603454.795...    0
    4  113.0  POLYGON ((1603524.561 6464388.609, 1603532.241...    0
    """

    def __init__(
        self,
        gdf,
        unique_id,
        limit=None,
        shrink=0.4,
        segment=0.5,
        verbose=True,
        enclosures=None,
        enclosure_id="eID",
        threshold=0.05,
        use_dask=True,
        n_chunks=None,
    ):
        if os.getenv("ALLOW_LEGACY_MOMEPY", "False").lower() not in (
            "true",
            "1",
            "yes",
        ):
            warnings.warn(
                "Class based API like `momepy.Tessellation` is deprecated. "
                "Replace it with `momepy.morphological_tessellation` or "
                "`momepy.enclosed_tessellation` to use functional API instead "
                "or pin momepy version <1.0. Class-based API will be removed in "
                "1.0. "
                # "See details at https://docs.momepy.org/en/stable/migration.html",
                "",
                FutureWarning,
                stacklevel=2,
            )
        self.gdf = gdf
        self.id = gdf[unique_id]
        self.limit = limit
        self.shrink = shrink
        self.segment = segment
        self.enclosure_id = enclosure_id

        if gdf.crs and gdf.crs.is_geographic:
            raise ValueError(
                "Geometry is in a geographic CRS. "
                "Use 'GeoDataFrame.to_crs()' to re-project geometries to a "
                "projected CRS before using Tessellation.",
            )

        if limit is not None and enclosures is not None:
            raise ValueError(
                "Both `limit` and `enclosures` cannot be passed together. "
                "Pass `limit` for morphological tessellation or `enclosures` "
                "for enclosed tessellation."
            )

        gdf = gdf.copy()

        if enclosures is not None:
            enclosures = enclosures.copy()

            bounds = enclosures.total_bounds
            centre_x = (bounds[0] + bounds[2]) / 2
            centre_y = (bounds[1] + bounds[3]) / 2

            gdf.geometry = gdf.geometry.translate(xoff=-centre_x, yoff=-centre_y)
            enclosures.geometry = enclosures.geometry.translate(
                xoff=-centre_x, yoff=-centre_y
            )

            self.tessellation = self._enclosed_tessellation(
                gdf,
                enclosures,
                unique_id,
                threshold,
                use_dask,
                n_chunks,
            )
        else:
            if isinstance(limit, gpd.GeoSeries | gpd.GeoDataFrame):
                limit = limit.union_all() if GPD_GE_10 else limit.unary_union

            bounds = shapely.bounds(limit)
            centre_x = (bounds[0] + bounds[2]) / 2
            centre_y = (bounds[1] + bounds[3]) / 2
            gdf.geometry = gdf.geometry.translate(xoff=-centre_x, yoff=-centre_y)

            # add convex hull buffered large distance to eliminate infinity issues
            limit = (
                gpd.GeoSeries(limit, crs=gdf.crs)
                .translate(xoff=-centre_x, yoff=-centre_y)
                .array[0]
            )

            self.tessellation = self._morphological_tessellation(
                gdf, unique_id, limit, shrink, segment, verbose
            )

        self.tessellation["geometry"] = self.tessellation["geometry"].translate(
            xoff=centre_x, yoff=centre_y
        )

    def _morphological_tessellation(
        self, gdf, unique_id, limit, shrink, segment, verbose, check=True
    ):
        objects = gdf

        if shrink != 0:
            print("Inward offset...") if verbose else None
            mask = objects.geom_type.isin(["Polygon", "MultiPolygon"])
            objects.loc[mask, objects.geometry.name] = objects[mask].buffer(
                -shrink, cap_style=2, join_style=2
            )
        objects = objects.reset_index(drop=True).explode(ignore_index=True)
        objects = objects.set_index(unique_id)

        print("Generating input point array...") if verbose else None
        points, ids = self._dense_point_array(
            objects.geometry.array, distance=segment, index=objects.index
        )

        hull = shapely.convex_hull(limit)
        bounds = shapely.bounds(hull)
        width = bounds[2] - bounds[0]
        leng = bounds[3] - bounds[1]
        hull = shapely.buffer(hull, 2 * width if width > leng else 2 * leng)

        hull_p, hull_ix = self._dense_point_array(
            [hull], distance=shapely.length(hull) / 100, index=[0]
        )
        points = np.append(points, hull_p, axis=0)
        ids = ids + ([-1] * len(hull_ix))

        print("Generating Voronoi diagram...") if verbose else None
        voronoi_diagram = Voronoi(np.array(points))

        print("Generating GeoDataFrame...") if verbose else None
        regions_gdf = self._regions(voronoi_diagram, unique_id, ids, crs=gdf.crs)

        print("Dissolving Voronoi polygons...") if verbose else None
        morphological_tessellation = regions_gdf[[unique_id, "geometry"]].dissolve(
            by=unique_id, as_index=False
        )

        morphological_tessellation = gpd.clip(
            morphological_tessellation, gpd.GeoSeries(limit, crs=gdf.crs)
        )

        if check:
            self._check_result(morphological_tessellation, gdf, unique_id=unique_id)

        return morphological_tessellation

    def _dense_point_array(self, geoms, distance, index):
        """
        geoms : array of shapely lines
        """
        # interpolate lines to represent them as points for Voronoi
        points = []
        ids = []

        if shapely.get_type_id(geoms[0]) not in [1, 2, 5]:
            lines = shapely.boundary(geoms)
        else:
            lines = geoms
        lengths = shapely.length(lines)
        for ix, line, length in zip(index, lines, lengths, strict=True):
            if length > distance:  # some polygons might have collapsed
                pts = shapely.line_interpolate_point(
                    line,
                    np.linspace(0.1, length - 0.1, num=int((length - 0.1) // distance)),
                )  # .1 offset to keep a gap between two segments
                points.append(shapely.get_coordinates(pts))
                ids += [ix] * len(pts)

        points = np.vstack(points)

        return points, ids

        # here we might also want to append original coordinates of each line
        # to get a higher precision on the corners

    def _regions(self, voronoi_diagram, unique_id, ids, crs):
        """Generate GeoDataFrame of Voronoi regions from scipy.spatial.Voronoi."""
        vertices = pd.Series(voronoi_diagram.regions).take(voronoi_diagram.point_region)
        polygons = []
        for region in vertices:
            if -1 not in region:
                polygons.append(shapely.polygons(voronoi_diagram.vertices[region]))
            else:
                polygons.append(None)

        regions_gdf = gpd.GeoDataFrame(
            {unique_id: ids}, geometry=polygons, crs=crs
        ).dropna()
        regions_gdf = regions_gdf.loc[
            regions_gdf[unique_id] != -1
        ]  # delete hull-based cells

        return regions_gdf

    def _check_result(self, tesselation, orig_gdf, unique_id):
        """Check whether result matches buildings and contains only Polygons."""
        # check against input layer
        ids_original = list(orig_gdf[unique_id])
        ids_generated = list(tesselation[unique_id])
        if len(ids_original) != len(ids_generated):
            self.collapsed = set(ids_original).difference(ids_generated)
            warnings.warn(
                message=(
                    "Tessellation does not fully match buildings. "
                    f"{len(self.collapsed)} element(s) collapsed "
                    f"during generation - unique_id: {self.collapsed}."
                ),
                category=UserWarning,
                stacklevel=4,
            )

        # check MultiPolygons - usually caused by error in input geometry
        self.multipolygons = tesselation[
            tesselation.geometry.geom_type == "MultiPolygon"
        ][unique_id]
        if len(self.multipolygons) > 0:
            warnings.warn(
                message=(
                    "Tessellation contains MultiPolygon elements. Initial "
                    "objects should  be edited. `unique_id` of affected "
                    f"elements: {list(self.multipolygons)}."
                ),
                category=UserWarning,
                stacklevel=4,
            )

    def _enclosed_tessellation(
        self,
        buildings,
        enclosures,
        unique_id,
        threshold=0.05,
        use_dask=True,
        n_chunks=None,
        **kwargs,
    ):
        """
        Generate enclosed tessellation based on barriers
        defining enclosures and building footprints.

        Parameters
        ----------
        buildings : GeoDataFrame
            A GeoDataFrame containing building footprints.
            Expects (Multi)Polygon geometry.
        enclosures : GeoDataFrame
            Enclosures geometry. Can  be generated using :func:`momepy.enclosures`.
        unique_id : str
            The name of the column with the unique ID of ``buildings`` gdf.
        threshold : float (default 0.05)
            The minimum threshold for a building to be considered within an enclosure.
            Threshold is a ratio of building area which needs to be within an enclosure
            to inlude it in the tessellation of that enclosure.
            Resolves sliver geometry issues.
        use_dask : bool (default True)
            Use parallelised algorithm based on ``dask.dataframe``. Requires dask.
        n_chunks : None
            The number of chunks to be used in parallelization. Ideal is one chunk per
            thread. Applies only if ``enclosures`` are passed. Default automatically
            uses ``n == dask.system.cpu_count``.
        **kwargs : dict
            Keyword arguments passed to Tessellation algorithm
            (such as ``'shrink'`` or ``'segment'``).

        Returns
        -------
        tessellation : GeoDataFrame
            A GeoDataFrame containing three columns:

                - ``geometry``,
                - ``unique_id`` matching with parent building,
                - ``enclosure_id`` matching with enclosure integer index

        Examples
        --------
        >>> enclosures = mm.enclosures(streets, admin_boundary, [railway, rivers])
        >>> enclosed_tess = mm.enclosed_tessellation(buildings, enclosures)
        """
        enclosures = enclosures.reset_index(drop=True)
        enclosures["position"] = range(len(enclosures))

        # determine which polygons should be split
        if GPD_GE_013:
            inp, res = buildings.sindex.query(
                enclosures.geometry, predicate="intersects"
            )
        else:
            inp, res = buildings.sindex.query_bulk(
                enclosures.geometry, predicate="intersects"
            )
        unique, counts = np.unique(inp, return_counts=True)
        splits = unique[counts > 1]
        single = unique[counts == 1]

        if use_dask:
            try:
                import dask.bag as db
                from dask.system import cpu_count
            except ImportError:
                use_dask = False

                warnings.warn(
                    message=(
                        "dask.dataframe could not be imported. "
                        f"Setting `use_dask={use_dask}`."
                    ),
                    category=UserWarning,
                    stacklevel=3,
                )

        if use_dask:
            if n_chunks is None:
                n_chunks = cpu_count() - 1 if cpu_count() > 1 else 1
            # initialize dask.bag
            bag = db.from_sequence(splits, npartitions=n_chunks)
            # generate enclosed tessellation using dask
            new = bag.map(
                self._tess,
                enclosures,
                buildings,
                inp,
                res,
                threshold,
                unique_id,
            ).compute()

        else:
            new = [
                self._tess(
                    i,
                    enclosures,
                    buildings,
                    inp,
                    res,
                    threshold=threshold,
                    unique_id=unique_id,
                    **kwargs,
                )
                for i in splits
            ]

        # finalise the result
        clean_blocks = enclosures.drop(splits)
        clean_blocks.loc[single, unique_id] = clean_blocks.loc[
            single, "position"
        ].apply(lambda ix: buildings.iloc[res[inp == ix][0]][unique_id])
        return pd.concat(new + [clean_blocks.drop(columns="position")]).reset_index(
            drop=True
        )

    def _tess(
        self,
        ix,
        enclosure,
        buildings,
        query_inp,
        query_res,
        threshold,
        unique_id,
    ):
        poly = enclosure.geometry.array[ix]
        blg = buildings.iloc[query_res[query_inp == ix]]
        within = blg[
            shapely.area(shapely.intersection(blg.geometry.array, poly))
            > (shapely.area(blg.geometry.array) * threshold)
        ].copy()
        if len(within) > 1:
            tess = self._morphological_tessellation(
                within,
                unique_id,
                poly,
                shrink=self.shrink,
                segment=self.segment,
                verbose=False,
                check=False,
            )
            tess[self.enclosure_id] = enclosure[self.enclosure_id].iloc[ix]
            return tess
        return gpd.GeoDataFrame(
            {self.enclosure_id: enclosure[self.enclosure_id].iloc[ix], unique_id: None},
            geometry=[poly],
            index=[0],
        )


@deprecated("generate_blocks")
class Blocks:
    """
    Generate blocks based on buildings, tessellation, and street network.
    Dissolves tessellation cells based on street-network based polygons.
    Links resulting ID to ``buildings`` and ``tessellation`` as attributes.

    Parameters
    ----------
    tessellation : GeoDataFrame
        A GeoDataFrame containing morphological tessellation.
    edges : GeoDataFrame
        A GeoDataFrame containing a street network.
    buildings : GeoDataFrame
        A GeoDataFrame containing buildings.
    id_name : str
        The name of the unique blocks ID column to be generated.
    unique_id : str
        The name of the column with the unique ID. If there is none, it can be
        generated with :func:`momepy.unique_id`. This should be the same for
        cells and buildings; ID's should match.

    Attributes
    ----------
    blocks : GeoDataFrame
        A GeoDataFrame containing generated blocks.
    buildings_id : Series
        A Series derived from buildings with block ID.
    tessellation_id : Series
        A Series derived from morphological tessellation with block ID.
    tessellation : GeoDataFrame
        A GeoDataFrame containing original tessellation.
    edges : GeoDataFrame
        A GeoDataFrame containing original edges.
    buildings : GeoDataFrame
        A GeoDataFrame containing original buildings.
    id_name : str
        The name of the unique blocks ID column.
    unique_id : str
        The name of the column with unique ID.

    Examples
    --------
    >>> blocks = mm.Blocks(tessellation_df, streets_df, buildings_df, 'bID', 'uID')
    >>> blocks.blocks.head()
        bID	geometry
    0	1.0	POLYGON ((1603560.078648818 6464202.366899694,...
    1	2.0	POLYGON ((1603457.225976106 6464299.454696888,...
    2	3.0	POLYGON ((1603056.595487018 6464093.903488506,...
    3	4.0	POLYGON ((1603260.943782872 6464141.327631323,...
    4	5.0	POLYGON ((1603183.399594798 6463966.109982309,...
    """

    def __init__(self, tessellation, edges, buildings, id_name, unique_id):
        self.tessellation = tessellation
        self.edges = edges
        self.buildings = buildings
        self.id_name = id_name
        self.unique_id = unique_id

        if id_name in buildings.columns:
            raise ValueError(
                f"'{id_name}' column cannot be in the buildings GeoDataFrame."
            )

        cut = gpd.overlay(
            tessellation,
            gpd.GeoDataFrame(geometry=edges.buffer(0.001)),
            how="difference",
        )
        cut = cut.explode(ignore_index=True)
        weights = libpysal.weights.Queen.from_dataframe(
            cut, silence_warnings=True, use_index=False
        )
        cut["component"] = weights.component_labels
        buildings_c = buildings.copy()
        buildings_c.geometry = buildings_c.representative_point()  # make points
        centroids_temp_id = gpd.sjoin(
            buildings_c,
            cut[[cut.geometry.name, "component"]],
            how="left",
            predicate="within",
        )

        cells_copy = tessellation[[unique_id, tessellation.geometry.name]].merge(
            centroids_temp_id[[unique_id, "component"]], on=unique_id, how="left"
        )
        blocks = cells_copy.dissolve(by="component").explode(ignore_index=True)
        blocks[id_name] = range(len(blocks))
        blocks = blocks[[id_name, blocks.geometry.name]]

        centroids_w_bl_id2 = gpd.sjoin(
            buildings_c, blocks, how="left", predicate="within"
        )

        self.buildings_id = centroids_w_bl_id2[id_name]

        cells_m = tessellation[[unique_id]].merge(
            centroids_w_bl_id2[[unique_id, id_name]], on=unique_id, how="left"
        )

        self.tessellation_id = cells_m[id_name]
        self.tessellation_id.index = self.tessellation.index

        self.blocks = blocks


@deprecated("get_nearest_street")
def get_network_id(left, right, network_id, min_size=100, verbose=True):
    """
    Snap each element (preferably building) to the closest
    street network segment and save its ID. Also, adds network ID to elements.

    Parameters
    ----------
    left : GeoDataFrame
        A GeoDataFrame containing objects to snap.
    right : GeoDataFrame
        A GeoDataFrame containing a street network with unique network IDs.
        If there is none, it can be generated with :func:`momepy.unique_id`.
    network_id : str, list, np.array, pd.Series (default None)
        The name of the streets dataframe column, ``np.array``, or ``pd.Series``
        with network unique IDs.
    min_size : int (default 100)
        A minimum size should be a valuee such that if you build a box centered in each
        building centroid with edges of size ``2*min_size``, you know a priori that at
        least one segment is intersected with the box.
    verbose : bool (default True)
        If ``True``, shows progress bars in loops and indication of steps.

    Returns
    -------
    elements_nID : Series
        A Series containing network ID for elements.

    Examples
    --------
    >>> buildings_df['nID'] = momepy.get_network_id(buildings_df, streets_df, 'nID')
    Generating centroids...
    Generating rtree...
    Snapping: 100%|██████████| 144/144 [00:00<00:00, 2718.98it/s]
    >>> buildings_df['nID'][0]
    1

    See also
    --------
    momepy.get_network_ratio
    momepy.get_node_id
    """
    infty = 1000000000000
    left = left.copy()
    right = right.copy()

    if not isinstance(network_id, str):
        right["mm_nid"] = network_id
        network_id = "mm_nid"

    buildings_c = left.copy()

    buildings_c[buildings_c.geometry.name] = buildings_c.centroid  # make centroids

    idx = right.sindex

    # TODO: use sjoin nearest once done
    result = []
    for p in tqdm(
        buildings_c.geometry,
        total=buildings_c.shape[0],
        desc="Snapping",
        disable=not verbose,
    ):
        pbox = (p.x - min_size, p.y - min_size, p.x + min_size, p.y + min_size)
        hits = list(idx.intersection(pbox))
        d = infty
        nid = None
        for h in hits:
            new_d = p.distance(right.geometry.iloc[h])
            if d >= new_d:
                d = new_d
                nid = right[network_id].iloc[h]
        if nid is None:
            result.append(np.nan)
        else:
            result.append(nid)

    series = pd.Series(result, index=left.index)

    if series.isnull().any():
        warnings.warn(
            message=(
                "Some objects were not attached to the network. Set larger "
                f"`min_size``. {sum(series.isnull())} affected elements."
            ),
            category=UserWarning,
            stacklevel=2,
        )
    return series


@deprecated("get_nearest_node")
def get_node_id(
    objects,
    nodes,
    edges,
    node_id,
    edge_id=None,
    edge_keys=None,
    edge_values=None,
    verbose=True,
):
    """
    Snap each building to the closest street network node on the closest network edge.
    Adds node ID to objects (preferably buildings). Gets ID of edge
    (:func:`momepy.get_network_id` or :func:`get_network_ratio`), and determines
    which of its end points is closer to the building centroid. Pass either ``edge_id``
    with a single value or ``edge_keys`` and ``edge_values`` with ratios.

    Parameters
    ----------
    objects : GeoDataFrame
        A GeoDataFrame containing objects to snap.
    nodes : GeoDataFrame
        A GeoDataFrame containing street nodes with unique node IDs.
        If there is none, it can be generated by :func:`momepy.unique_id`.
    edges : GeoDataFrame
        A GeoDataFrame containing street edges with unique edge IDs and IDs
        of start and end points of each segment. Start and endpoints are default
        outcome of :func:`momepy.nx_to_gdf`.
    node_id : str, list, np.array, pd.Series
        The name of the ``nodes`` dataframe column, ``np.array``,
        or ``pd.Series`` with a unique ID.
    edge_id : str (default None)
        The name of the objects dataframe column with unique edge IDs
        (like an outcome of :func:`momepy.get_network_id`).
    edge_keys : str (default None)
        The name of the objects dataframe column with ``edgeID_keys``
        (like an outcome of :func:`momepy.get_network_ratio`).
    edge_values : str (default None)
        The name of the objects dataframe column with ``edgeID_values``
        (like an outcome of :func:`momepy.get_network_ratio`).
    verbose : bool (default True)
        If ``True``, shows progress bars in loops and indication of steps.

    Returns
    -------
    node_ids : Series
        A Series containing node the ID for objects.
    """
    nodes = nodes.set_index(node_id)

    if not isinstance(node_id, str):
        nodes["mm_noid"] = node_id
        node_id = "mm_noid"

    results_list = []
    if edge_id is not None:
        edges = edges.set_index(edge_id)
        centroids = objects.centroid
        for eid, centroid in tqdm(
            zip(objects[edge_id], centroids, strict=True),
            total=objects.shape[0],
            disable=not verbose,
        ):
            if pd.isna(eid):
                results_list.append(pd.NA)
            else:
                edge = edges.loc[eid]
                start_id = edge.node_start
                start = nodes.loc[start_id].geometry
                sd = centroid.distance(start)
                end_id = edge.node_end
                end = nodes.loc[end_id].geometry
                ed = centroid.distance(end)
                if sd > ed:
                    results_list.append(end_id)
                else:
                    results_list.append(start_id)

    elif edge_keys is not None and edge_values is not None:
        for edge_i, edge_r, geom in tqdm(
            zip(
                objects[edge_keys], objects[edge_values], objects.geometry, strict=True
            ),
            total=objects.shape[0],
            disable=not verbose,
        ):
            edge = edges.iloc[edge_i[edge_r.index(max(edge_r))]]
            start_id = edge.node_start
            start = nodes.loc[start_id].geometry
            sd = geom.distance(start)
            end_id = edge.node_end
            end = nodes.loc[end_id].geometry
            ed = geom.distance(end)
            if sd > ed:
                results_list.append(end_id)
            else:
                results_list.append(start_id)

    series = pd.Series(results_list, index=objects.index)
    return series


def get_network_ratio(df, edges, initial_buffer=500):
    """
    Link polygons to network edges based on the proportion of overlap (if a cell
    intersects more than one edge). Useful if you need to link enclosed tessellation to
    street network. Ratios can be used as weights when linking network-based values
    to cells. For a purely distance-based link use :func:`momepy.get_network_id`.
    Links are based on the integer position of edge (``iloc``).

    Parameters
    ----------
    df : GeoDataFrame
        A GeoDataFrame containing objects to snap (typically enclosed tessellation).
    edges : GeoDataFrame
        A GeoDataFrame containing a street network.
    initial_buffer : float
        The initial buffer used to link non-intersecting cells.

    Returns
    -------
    result : DataFrame
        The resultant DataFrame.

    See also
    --------
    momepy.get_network_id
    momepy.get_node_id

    Examples
    --------
    >>> links = mm.get_network_ratio(enclosed_tessellation, streets)
    >>> links.head()
      edgeID_keys                              edgeID_values
    0        [34]                                      [1.0]
    1     [0, 34]  [0.38508998545027145, 0.6149100145497285]
    2        [32]                                        [1]
    3         [0]                                      [1.0]
    4        [26]                                        [1]
    """

    (df_ix, edg_ix), dist = edges.sindex.nearest(
        df.geometry, max_distance=initial_buffer, return_distance=True
    )

    touching = dist < 0.1

    intersections = (
        df.iloc[df_ix[touching]]
        .intersection(edges.buffer(0.0001).iloc[edg_ix[touching]], align=False)
        .reset_index()
    )

    mask = intersections.area > 0.0001

    df_ix_touching = df_ix[touching][mask]
    lengths = intersections[mask].area
    grouped = lengths.groupby(df_ix_touching)
    totals = grouped.sum()
    ints_vect = []
    for name, group in grouped:
        ratios = group / totals.loc[name]
        ints_vect.append(
            {edg_ix[touching][item[0]]: item[1] for item in ratios.items()}
        )

    ratios = pd.Series(ints_vect, index=df.index[list(grouped.groups.keys())])

    near = []
    df_ix_non = df_ix[~touching]
    grouped = pd.Series(dist[~touching]).groupby(df_ix_non)
    for _, group in grouped:
        near.append({edg_ix[~touching][group.idxmin()]: 1.0})

    near = pd.Series(near, index=df.index[list(grouped.groups.keys())])

    ratios = pd.concat([ratios, near])

    nans = df[~df.index.isin(ratios.index)]
    if not nans.empty:
        df_ix, edg_ix = edges.sindex.nearest(
            nans.geometry, return_all=False, max_distance=None
        )
        additional = pd.Series([{i: 1.0} for i in edg_ix], index=nans.index)

        ratios = pd.concat([ratios, additional])

    result = pd.DataFrame()
    result["edgeID_keys"] = ratios.apply(lambda d: list(d.keys()))
    result["edgeID_values"] = ratios.apply(lambda d: list(d.values()))

    return result


def enclosures(
    primary_barriers,
    limit=None,
    additional_barriers=None,
    enclosure_id="eID",
    clip=False,
):
    """
    Generate enclosures based on passed barriers. Enclosures are areas enclosed from
    all sides by at least one type of a barrier. Barriers are typically roads,
    railways, natural features like rivers and other water bodies or coastline.
    Enclosures are a result of polygonization of the  ``primary_barrier`` and ``limit``
    and its subdivision based on additional_barriers.

    Parameters
    ----------
    primary_barriers : GeoDataFrame, GeoSeries
        A GeoDataFrame or GeoSeries containing primary barriers.
        (Multi)LineString geometry is expected.
    limit : GeoDataFrame, GeoSeries, shapely geometry (default None)
        A GeoDataFrame, GeoSeries or shapely geometry containing external limit
        of enclosures, i.e. the area which gets partitioned. If ``None`` is passed,
        the internal area of ``primary_barriers`` will be used.
    additional_barriers : GeoDataFrame
        A GeoDataFrame or GeoSeries containing additional barriers.
        (Multi)LineString geometry is expected.
    enclosure_id : str (default 'eID')
        The name of the ``enclosure_id`` (to be created).
    clip : bool (default False)
        If ``True``, returns enclosures with representative point within the limit
        (if given). Requires ``limit`` composed of Polygon or MultiPolygon geometries.

    Returns
    -------
    enclosures : GeoDataFrame
       A GeoDataFrame containing enclosure geometries and ``enclosure_id``.

    Examples
    --------
    >>> enclosures = mm.enclosures(streets, admin_boundary, [railway, rivers])
    """
    if limit is not None:
        if isinstance(limit, BaseGeometry):
            limit = gpd.GeoSeries([limit], crs=primary_barriers.crs)
        if limit.geom_type.isin(["Polygon", "MultiPolygon"]).any():
            limit_b = limit.boundary
        else:
            limit_b = limit
        barriers = pd.concat([primary_barriers.geometry, limit_b.geometry])
    else:
        barriers = primary_barriers
    unioned = barriers.union_all() if GPD_GE_10 else barriers.unary_union
    polygons = polygonize(unioned)
    enclosures = gpd.GeoSeries(list(polygons), crs=primary_barriers.crs)

    if additional_barriers is not None:
        if not isinstance(additional_barriers, list):
            raise TypeError(
                "`additional_barriers` expects a list of GeoDataFrames "
                f"or GeoSeries. Got {type(additional_barriers)}."
            )
        additional = pd.concat([gdf.geometry for gdf in additional_barriers])

        if GPD_GE_013:
            inp, res = enclosures.sindex.query(
                additional.geometry, predicate="intersects"
            )
        else:
            inp, res = enclosures.sindex.query_bulk(
                additional.geometry, predicate="intersects"
            )
        unique = np.unique(res)

        new = []

        for i in unique:
            poly = enclosures.array[i]  # get enclosure polygon
            crossing = inp[res == i]  # get relevant additional barriers
            buf = shapely.buffer(poly, 0.01)  # to avoid floating point errors
            crossing_ins = shapely.intersection(
                buf, additional.array[crossing]
            )  # keeping only parts of additional barriers within polygon
            union = shapely.union_all(
                np.append(crossing_ins, shapely.boundary(poly))
            )  # union
            polygons = shapely.get_parts(shapely.polygonize([union]))  # polygonize
            within = shapely.covered_by(
                polygons, buf
            )  # keep only those within original polygon
            new += list(polygons[within])

        final_enclosures = pd.concat(
            [
                gpd.GeoSeries(enclosures).drop(unique),
                gpd.GeoSeries(new, crs=primary_barriers.crs),
            ]
        ).reset_index(drop=True)

        final_enclosures = gpd.GeoDataFrame(
            {enclosure_id: range(len(final_enclosures))}, geometry=final_enclosures
        )

    else:
        final_enclosures = gpd.GeoDataFrame(
            {enclosure_id: range(len(enclosures))}, geometry=enclosures
        )

    if clip and limit is not None:
        if not limit.geom_type.isin(["Polygon", "MultiPolygon"]).all():
            raise TypeError(
                "`limit` requires a GeoDataFrame or GeoSeries with Polygon or "
                "MultiPolygon geometry to be used with `clip=True`."
            )
        if GPD_GE_013:
            _, encl_index = final_enclosures.representative_point().sindex.query(
                limit.geometry, predicate="contains"
            )
        else:
            _, encl_index = final_enclosures.representative_point().sindex.query_bulk(
                limit.geometry, predicate="contains"
            )
        keep = np.unique(encl_index)
        return final_enclosures.iloc[keep]

    return final_enclosures