initial commit

billinglayer/python/shapely/__init__.py

@@ -0,0 +1,33 @@
+from .lib import GEOSException  # NOQA
+from .lib import Geometry  # NOQA
+from .lib import geos_version, geos_version_string  # NOQA
+from .lib import geos_capi_version, geos_capi_version_string  # NOQA
+from .errors import setup_signal_checks  # NOQA
+from ._geometry import *  # NOQA
+from .creation import *  # NOQA
+from .constructive import *  # NOQA
+from .predicates import *  # NOQA
+from .measurement import *  # NOQA
+from .set_operations import *  # NOQA
+from .linear import *  # NOQA
+from .coordinates import *  # NOQA
+from .strtree import *  # NOQA
+from .io import *  # NOQA
+
+# Submodule always needs to be imported to ensure Geometry subclasses are registered
+from shapely.geometry import (  # NOQA
+    Point,
+    LineString,
+    Polygon,
+    MultiPoint,
+    MultiLineString,
+    MultiPolygon,
+    GeometryCollection,
+    LinearRing,
+)
+
+from . import _version
+
+__version__ = _version.get_versions()["version"]
+
+setup_signal_checks()

billinglayer/python/shapely/_enum.py

@@ -0,0 +1,23 @@
+from enum import IntEnum
+
+
+class ParamEnum(IntEnum):
+    """Wraps IntEnum to provide validation of a requested item.
+
+    Intended for enums used for function parameters.
+
+    Use enum.get_value(item) for this behavior instead of builtin enum[item].
+    """
+
+    @classmethod
+    def get_value(cls, item):
+        """Validate incoming item and raise a ValueError with valid options if not present."""
+        try:
+            return cls[item].value
+        except KeyError:
+            valid_options = {e.name for e in cls}
+            raise ValueError(
+                "'{}' is not a valid option, must be one of '{}'".format(
+                    item, "', '".join(valid_options)
+                )
+            )

billinglayer/python/shapely/_geometry.py

@@ -0,0 +1,886 @@
+import warnings
+from enum import IntEnum
+
+import numpy as np
+
+from . import _geometry_helpers, geos_version, lib
+from ._enum import ParamEnum
+from .decorators import multithreading_enabled, requires_geos
+
+__all__ = [
+    "GeometryType",
+    "get_type_id",
+    "get_dimensions",
+    "get_coordinate_dimension",
+    "get_num_coordinates",
+    "get_srid",
+    "set_srid",
+    "get_x",
+    "get_y",
+    "get_z",
+    "get_exterior_ring",
+    "get_num_points",
+    "get_num_interior_rings",
+    "get_num_geometries",
+    "get_point",
+    "get_interior_ring",
+    "get_geometry",
+    "get_parts",
+    "get_rings",
+    "get_precision",
+    "set_precision",
+    "force_2d",
+    "force_3d",
+]
+
+
+class GeometryType(IntEnum):
+    """The enumeration of GEOS geometry types"""
+
+    MISSING = -1
+    POINT = 0
+    LINESTRING = 1
+    LINEARRING = 2
+    POLYGON = 3
+    MULTIPOINT = 4
+    MULTILINESTRING = 5
+    MULTIPOLYGON = 6
+    GEOMETRYCOLLECTION = 7
+
+
+# generic
+
+
+@multithreading_enabled
+def get_type_id(geometry, **kwargs):
+    """Returns the type ID of a geometry.
+
+    - None (missing) is -1
+    - POINT is 0
+    - LINESTRING is 1
+    - LINEARRING is 2
+    - POLYGON is 3
+    - MULTIPOINT is 4
+    - MULTILINESTRING is 5
+    - MULTIPOLYGON is 6
+    - GEOMETRYCOLLECTION is 7
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    GeometryType
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point
+    >>> get_type_id(LineString([(0, 0), (1, 1), (2, 2), (3, 3)]))
+    1
+    >>> get_type_id([Point(1, 2), Point(2, 3)]).tolist()
+    [0, 0]
+    """
+    return lib.get_type_id(geometry, **kwargs)
+
+
+@multithreading_enabled
+def get_dimensions(geometry, **kwargs):
+    """Returns the inherent dimensionality of a geometry.
+
+    The inherent dimension is 0 for points, 1 for linestrings and linearrings,
+    and 2 for polygons. For geometrycollections it is the max of the containing
+    elements. Empty collections and None values return -1.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import GeometryCollection, Point, Polygon
+    >>> point = Point(0, 0)
+    >>> get_dimensions(point)
+    0
+    >>> polygon = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])
+    >>> get_dimensions(polygon)
+    2
+    >>> get_dimensions(GeometryCollection([point, polygon]))
+    2
+    >>> get_dimensions(GeometryCollection([]))
+    -1
+    >>> get_dimensions(None)
+    -1
+    """
+    return lib.get_dimensions(geometry, **kwargs)
+
+
+@multithreading_enabled
+def get_coordinate_dimension(geometry, **kwargs):
+    """Returns the dimensionality of the coordinates in a geometry (2 or 3).
+
+    Returns -1 for missing geometries (``None`` values). Note that if the first Z
+    coordinate equals ``nan``, this function will return ``2``.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import Point
+    >>> get_coordinate_dimension(Point(0, 0))
+    2
+    >>> get_coordinate_dimension(Point(0, 0, 1))
+    3
+    >>> get_coordinate_dimension(None)
+    -1
+    >>> get_coordinate_dimension(Point(0, 0, float("nan")))
+    2
+    """
+    return lib.get_coordinate_dimension(geometry, **kwargs)
+
+
+@multithreading_enabled
+def get_num_coordinates(geometry, **kwargs):
+    """Returns the total number of coordinates in a geometry.
+
+    Returns 0 for not-a-geometry values.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import GeometryCollection, LineString, Point
+    >>> point = Point(0, 0)
+    >>> get_num_coordinates(point)
+    1
+    >>> get_num_coordinates(Point(0, 0, 0))
+    1
+    >>> line = LineString([(0, 0), (1, 1)])
+    >>> get_num_coordinates(line)
+    2
+    >>> get_num_coordinates(GeometryCollection([point, line]))
+    3
+    >>> get_num_coordinates(None)
+    0
+    """
+    return lib.get_num_coordinates(geometry, **kwargs)
+
+
+@multithreading_enabled
+def get_srid(geometry, **kwargs):
+    """Returns the SRID of a geometry.
+
+    Returns -1 for not-a-geometry values.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    set_srid
+
+    Examples
+    --------
+    >>> from shapely import Point
+    >>> point = Point(0, 0)
+    >>> get_srid(point)
+    0
+    >>> with_srid = set_srid(point, 4326)
+    >>> get_srid(with_srid)
+    4326
+    """
+    return lib.get_srid(geometry, **kwargs)
+
+
+@multithreading_enabled
+def set_srid(geometry, srid, **kwargs):
+    """Returns a geometry with its SRID set.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    srid : int
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_srid
+
+    Examples
+    --------
+    >>> from shapely import Point
+    >>> point = Point(0, 0)
+    >>> get_srid(point)
+    0
+    >>> with_srid = set_srid(point, 4326)
+    >>> get_srid(with_srid)
+    4326
+    """
+    return lib.set_srid(geometry, np.intc(srid), **kwargs)
+
+
+# points
+
+
+@multithreading_enabled
+def get_x(point, **kwargs):
+    """Returns the x-coordinate of a point
+
+    Parameters
+    ----------
+    point : Geometry or array_like
+        Non-point geometries will result in NaN being returned.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_y, get_z
+
+    Examples
+    --------
+    >>> from shapely import MultiPoint, Point
+    >>> get_x(Point(1, 2))
+    1.0
+    >>> get_x(MultiPoint([(1, 1), (1, 2)]))
+    nan
+    """
+    return lib.get_x(point, **kwargs)
+
+
+@multithreading_enabled
+def get_y(point, **kwargs):
+    """Returns the y-coordinate of a point
+
+    Parameters
+    ----------
+    point : Geometry or array_like
+        Non-point geometries will result in NaN being returned.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_x, get_z
+
+    Examples
+    --------
+    >>> from shapely import MultiPoint, Point
+    >>> get_y(Point(1, 2))
+    2.0
+    >>> get_y(MultiPoint([(1, 1), (1, 2)]))
+    nan
+    """
+    return lib.get_y(point, **kwargs)
+
+
+@requires_geos("3.7.0")
+@multithreading_enabled
+def get_z(point, **kwargs):
+    """Returns the z-coordinate of a point.
+
+    Parameters
+    ----------
+    point : Geometry or array_like
+        Non-point geometries or geometries without 3rd dimension will result
+        in NaN being returned.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_x, get_y
+
+    Examples
+    --------
+    >>> from shapely import MultiPoint, Point
+    >>> get_z(Point(1, 2, 3))
+    3.0
+    >>> get_z(Point(1, 2))
+    nan
+    >>> get_z(MultiPoint([(1, 1, 1), (2, 2, 2)]))
+    nan
+    """
+    return lib.get_z(point, **kwargs)
+
+
+# linestrings
+
+
+@multithreading_enabled
+def get_point(geometry, index, **kwargs):
+    """Returns the nth point of a linestring or linearring.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    index : int or array_like
+        Negative values count from the end of the linestring backwards.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_num_points
+
+    Examples
+    --------
+    >>> from shapely import LinearRing, LineString, MultiPoint, Point
+    >>> line = LineString([(0, 0), (1, 1), (2, 2), (3, 3)])
+    >>> get_point(line, 1)
+    <POINT (1 1)>
+    >>> get_point(line, -2)
+    <POINT (2 2)>
+    >>> get_point(line, [0, 3]).tolist()
+    [<POINT (0 0)>, <POINT (3 3)>]
+
+    The functcion works the same for LinearRing input:
+
+    >>> get_point(LinearRing([(0, 0), (1, 1), (2, 2), (0, 0)]), 1)
+    <POINT (1 1)>
+
+    For non-linear geometries it returns None:
+
+    >>> get_point(MultiPoint([(0, 0), (1, 1), (2, 2), (3, 3)]), 1) is None
+    True
+    >>> get_point(Point(1, 1), 0) is None
+    True
+    """
+    return lib.get_point(geometry, np.intc(index), **kwargs)
+
+
+@multithreading_enabled
+def get_num_points(geometry, **kwargs):
+    """Returns number of points in a linestring or linearring.
+
+    Returns 0 for not-a-geometry values.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+        The number of points in geometries other than linestring or linearring
+        equals zero.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_point
+    get_num_geometries
+
+    Examples
+    --------
+    >>> from shapely import LineString, MultiPoint
+    >>> get_num_points(LineString([(0, 0), (1, 1), (2, 2), (3, 3)]))
+    4
+    >>> get_num_points(MultiPoint([(0, 0), (1, 1), (2, 2), (3, 3)]))
+    0
+    >>> get_num_points(None)
+    0
+    """
+    return lib.get_num_points(geometry, **kwargs)
+
+
+# polygons
+
+
+@multithreading_enabled
+def get_exterior_ring(geometry, **kwargs):
+    """Returns the exterior ring of a polygon.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_interior_ring
+
+    Examples
+    --------
+    >>> from shapely import Point, Polygon
+    >>> get_exterior_ring(Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)]))
+    <LINEARRING (0 0, 0 10, 10 10, 10 0, 0 0)>
+    >>> get_exterior_ring(Point(1, 1)) is None
+    True
+    """
+    return lib.get_exterior_ring(geometry, **kwargs)
+
+
+@multithreading_enabled
+def get_interior_ring(geometry, index, **kwargs):
+    """Returns the nth interior ring of a polygon.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    index : int or array_like
+        Negative values count from the end of the interior rings backwards.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_exterior_ring
+    get_num_interior_rings
+
+    Examples
+    --------
+    >>> from shapely import Point, Polygon
+    >>> polygon_with_hole = Polygon(
+    ...     [(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)],
+    ...     holes=[[(2, 2), (2, 4), (4, 4), (4, 2), (2, 2)]]
+    ... )
+    >>> get_interior_ring(polygon_with_hole, 0)
+    <LINEARRING (2 2, 2 4, 4 4, 4 2, 2 2)>
+    >>> get_interior_ring(polygon_with_hole, 1) is None
+    True
+    >>> polygon = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])
+    >>> get_interior_ring(polygon, 0) is None
+    True
+    >>> get_interior_ring(Point(0, 0), 0) is None
+    True
+    """
+    return lib.get_interior_ring(geometry, np.intc(index), **kwargs)
+
+
+@multithreading_enabled
+def get_num_interior_rings(geometry, **kwargs):
+    """Returns number of internal rings in a polygon
+
+    Returns 0 for not-a-geometry values.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+        The number of interior rings in non-polygons equals zero.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_exterior_ring
+    get_interior_ring
+
+    Examples
+    --------
+    >>> from shapely import Point, Polygon
+    >>> polygon = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])
+    >>> get_num_interior_rings(polygon)
+    0
+    >>> polygon_with_hole = Polygon(
+    ...     [(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)],
+    ...     holes=[[(2, 2), (2, 4), (4, 4), (4, 2), (2, 2)]]
+    ... )
+    >>> get_num_interior_rings(polygon_with_hole)
+    1
+    >>> get_num_interior_rings(Point(0, 0))
+    0
+    >>> get_num_interior_rings(None)
+    0
+    """
+    return lib.get_num_interior_rings(geometry, **kwargs)
+
+
+# collections
+
+
+@multithreading_enabled
+def get_geometry(geometry, index, **kwargs):
+    """Returns the nth geometry from a collection of geometries.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    index : int or array_like
+        Negative values count from the end of the collection backwards.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Notes
+    -----
+    - simple geometries act as length-1 collections
+    - out-of-range values return None
+
+    See also
+    --------
+    get_num_geometries, get_parts
+
+    Examples
+    --------
+    >>> from shapely import Point, MultiPoint
+    >>> multipoint = MultiPoint([(0, 0), (1, 1), (2, 2), (3, 3)])
+    >>> get_geometry(multipoint, 1)
+    <POINT (1 1)>
+    >>> get_geometry(multipoint, -1)
+    <POINT (3 3)>
+    >>> get_geometry(multipoint, 5) is None
+    True
+    >>> get_geometry(Point(1, 1), 0)
+    <POINT (1 1)>
+    >>> get_geometry(Point(1, 1), 1) is None
+    True
+    """
+    return lib.get_geometry(geometry, np.intc(index), **kwargs)
+
+
+def get_parts(geometry, return_index=False):
+    """Gets parts of each GeometryCollection or Multi* geometry object; returns
+    a copy of each geometry in the GeometryCollection or Multi* geometry object.
+
+    Note: This does not return the individual parts of Multi* geometry objects in
+    a GeometryCollection.  You may need to call this function multiple times to
+    return individual parts of Multi* geometry objects in a GeometryCollection.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    return_index : bool, default False
+        If True, will return a tuple of ndarrays of (parts, indexes), where indexes
+        are the indexes of the original geometries in the source array.
+
+    Returns
+    -------
+    ndarray of parts or tuple of (parts, indexes)
+
+    See also
+    --------
+    get_geometry, get_rings
+
+    Examples
+    --------
+    >>> from shapely import MultiPoint
+    >>> get_parts(MultiPoint([(0, 1), (2, 3)])).tolist()
+    [<POINT (0 1)>, <POINT (2 3)>]
+    >>> parts, index = get_parts([MultiPoint([(0, 1)]), MultiPoint([(4, 5), (6, 7)])], \
+return_index=True)
+    >>> parts.tolist()
+    [<POINT (0 1)>, <POINT (4 5)>, <POINT (6 7)>]
+    >>> index.tolist()
+    [0, 1, 1]
+    """
+    geometry = np.asarray(geometry, dtype=np.object_)
+    geometry = np.atleast_1d(geometry)
+
+    if geometry.ndim != 1:
+        raise ValueError("Array should be one dimensional")
+
+    if return_index:
+        return _geometry_helpers.get_parts(geometry)
+
+    return _geometry_helpers.get_parts(geometry)[0]
+
+
+def get_rings(geometry, return_index=False):
+    """Gets rings of Polygon geometry object.
+
+    For each Polygon, the first returned ring is always the exterior ring
+    and potential subsequent rings are interior rings.
+
+    If the geometry is not a Polygon, nothing is returned (empty array for
+    scalar geometry input or no element in output array for array input).
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    return_index : bool, default False
+        If True, will return a tuple of ndarrays of (rings, indexes), where
+        indexes are the indexes of the original geometries in the source array.
+
+    Returns
+    -------
+    ndarray of rings or tuple of (rings, indexes)
+
+    See also
+    --------
+    get_exterior_ring, get_interior_ring, get_parts
+
+    Examples
+    --------
+    >>> from shapely import Polygon
+    >>> polygon_with_hole = Polygon(
+    ...     [(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)],
+    ...     holes=[[(2, 2), (2, 4), (4, 4), (4, 2), (2, 2)]]
+    ... )
+    >>> get_rings(polygon_with_hole).tolist()
+    [<LINEARRING (0 0, 0 10, 10 10, 10 0, 0 0)>,
+     <LINEARRING (2 2, 2 4, 4 4, 4 2, 2 2)>]
+
+    With ``return_index=True``:
+
+    >>> polygon = Polygon([(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)])
+    >>> rings, index = get_rings([polygon, polygon_with_hole], return_index=True)
+    >>> rings.tolist()
+    [<LINEARRING (0 0, 2 0, 2 2, 0 2, 0 0)>,
+     <LINEARRING (0 0, 0 10, 10 10, 10 0, 0 0)>,
+     <LINEARRING (2 2, 2 4, 4 4, 4 2, 2 2)>]
+    >>> index.tolist()
+    [0, 1, 1]
+    """
+    geometry = np.asarray(geometry, dtype=np.object_)
+    geometry = np.atleast_1d(geometry)
+
+    if geometry.ndim != 1:
+        raise ValueError("Array should be one dimensional")
+
+    if return_index:
+        return _geometry_helpers.get_parts(geometry, extract_rings=True)
+
+    return _geometry_helpers.get_parts(geometry, extract_rings=True)[0]
+
+
+@multithreading_enabled
+def get_num_geometries(geometry, **kwargs):
+    """Returns number of geometries in a collection.
+
+    Returns 0 for not-a-geometry values.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+        The number of geometries in points, linestrings, linearrings and
+        polygons equals one.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_num_points
+    get_geometry
+
+    Examples
+    --------
+    >>> from shapely import MultiPoint, Point
+    >>> get_num_geometries(MultiPoint([(0, 0), (1, 1), (2, 2), (3, 3)]))
+    4
+    >>> get_num_geometries(Point(1, 1))
+    1
+    >>> get_num_geometries(None)
+    0
+    """
+    return lib.get_num_geometries(geometry, **kwargs)
+
+
+@requires_geos("3.6.0")
+@multithreading_enabled
+def get_precision(geometry, **kwargs):
+    """Get the precision of a geometry.
+
+    If a precision has not been previously set, it will be 0 (double
+    precision). Otherwise, it will return the precision grid size that was
+    set on a geometry.
+
+    Returns NaN for not-a-geometry values.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    set_precision
+
+    Examples
+    --------
+    >>> from shapely import Point
+    >>> point = Point(1, 1)
+    >>> get_precision(point)
+    0.0
+    >>> geometry = set_precision(point, 1.0)
+    >>> get_precision(geometry)
+    1.0
+    >>> get_precision(None)
+    nan
+    """
+    return lib.get_precision(geometry, **kwargs)
+
+
+class SetPrecisionMode(ParamEnum):
+    valid_output = 0
+    pointwise = 1
+    keep_collapsed = 2
+
+
+@requires_geos("3.6.0")
+@multithreading_enabled
+def set_precision(geometry, grid_size, mode="valid_output", **kwargs):
+    """Returns geometry with the precision set to a precision grid size.
+
+    By default, geometries use double precision coordinates (grid_size = 0).
+
+    Coordinates will be rounded if a precision grid is less precise than the
+    input geometry. Duplicated vertices will be dropped from lines and
+    polygons for grid sizes greater than 0. Line and polygon geometries may
+    collapse to empty geometries if all vertices are closer together than
+    grid_size. Z values, if present, will not be modified.
+
+    Note: subsequent operations will always be performed in the precision of
+    the geometry with higher precision (smaller "grid_size"). That same
+    precision will be attached to the operation outputs.
+
+    Also note: input geometries should be geometrically valid; unexpected
+    results may occur if input geometries are not.
+
+    Returns None if geometry is None.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    grid_size : float
+        Precision grid size. If 0, will use double precision (will not modify
+        geometry if precision grid size was not previously set). If this
+        value is more precise than input geometry, the input geometry will
+        not be modified.
+    mode :  {'valid_output', 'pointwise', 'keep_collapsed'}, default 'valid_output'
+        This parameter determines how to handle invalid output geometries. There are three modes:
+
+        1. `'valid_output'` (default):  The output is always valid. Collapsed geometry elements
+           (including both polygons and lines) are removed. Duplicate vertices are removed.
+        2. `'pointwise'`: Precision reduction is performed pointwise. Output geometry
+           may be invalid due to collapse or self-intersection. Duplicate vertices are not
+           removed. In GEOS this option is called NO_TOPO.
+
+           .. note::
+
+             'pointwise' mode requires at least GEOS 3.10. It is accepted in earlier versions,
+             but the results may be unexpected.
+        3. `'keep_collapsed'`: Like the default mode, except that collapsed linear geometry
+           elements are preserved. Collapsed polygonal input elements are removed. Duplicate
+           vertices are removed.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_precision
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point
+    >>> set_precision(Point(0.9, 0.9), 1.0)
+    <POINT (1 1)>
+    >>> set_precision(Point(0.9, 0.9, 0.9), 1.0)
+    <POINT Z (1 1 0.9)>
+    >>> set_precision(LineString([(0, 0), (0, 0.1), (0, 1), (1, 1)]), 1.0)
+    <LINESTRING (0 0, 0 1, 1 1)>
+    >>> set_precision(LineString([(0, 0), (0, 0.1), (0.1, 0.1)]), 1.0, mode="valid_output")
+    <LINESTRING Z EMPTY>
+    >>> set_precision(LineString([(0, 0), (0, 0.1), (0.1, 0.1)]), 1.0, mode="pointwise")
+    <LINESTRING (0 0, 0 0, 0 0)>
+    >>> set_precision(LineString([(0, 0), (0, 0.1), (0.1, 0.1)]), 1.0, mode="keep_collapsed")
+    <LINESTRING (0 0, 0 0)>
+    >>> set_precision(None, 1.0) is None
+    True
+    """
+    if isinstance(mode, str):
+        mode = SetPrecisionMode.get_value(mode)
+    elif not np.isscalar(mode):
+        raise TypeError("mode only accepts scalar values")
+    if mode == SetPrecisionMode.pointwise and geos_version < (3, 10, 0):
+        warnings.warn(
+            "'pointwise' is only supported for GEOS 3.10",
+            UserWarning,
+            stacklevel=2,
+        )
+    return lib.set_precision(geometry, grid_size, np.intc(mode), **kwargs)
+
+
+@multithreading_enabled
+def force_2d(geometry, **kwargs):
+    """Forces the dimensionality of a geometry to 2D.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point, Polygon, from_wkt
+    >>> force_2d(Point(0, 0, 1))
+    <POINT (0 0)>
+    >>> force_2d(Point(0, 0))
+    <POINT (0 0)>
+    >>> force_2d(LineString([(0, 0, 0), (0, 1, 1), (1, 1, 2)]))
+    <LINESTRING (0 0, 0 1, 1 1)>
+    >>> force_2d(from_wkt("POLYGON Z EMPTY"))
+    <POLYGON EMPTY>
+    >>> force_2d(None) is None
+    True
+    """
+    return lib.force_2d(geometry, **kwargs)
+
+
+@multithreading_enabled
+def force_3d(geometry, z=0.0, **kwargs):
+    """Forces the dimensionality of a geometry to 3D.
+
+    2D geometries will get the provided Z coordinate; Z coordinates of 3D geometries
+    are unchanged (unless they are nan).
+
+    Note that for empty geometries, 3D is only supported since GEOS 3.9 and then
+    still only for simple geometries (non-collections).
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    z : float or array_like, default 0.0
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point
+    >>> force_3d(Point(0, 0), z=3)
+    <POINT Z (0 0 3)>
+    >>> force_3d(Point(0, 0, 0), z=3)
+    <POINT Z (0 0 0)>
+    >>> force_3d(LineString([(0, 0), (0, 1), (1, 1)]))
+    <LINESTRING Z (0 0 0, 0 1 0, 1 1 0)>
+    >>> force_3d(None) is None
+    True
+    """
+    if np.isnan(z).any():
+        raise ValueError("It is not allowed to set the Z coordinate to NaN.")
+    return lib.force_3d(geometry, z, **kwargs)

billinglayer/python/shapely/_geos.pxd

@@ -0,0 +1,51 @@
+"""
+Provides a wrapper for GEOS types and functions.
+
+Note: GEOS functions in Cython must be called using the get_geos_handle context manager.
+Example:
+    with get_geos_handle() as geos_handle:
+        SomeGEOSFunc(geos_handle, ...<other params>)
+"""
+
+cdef extern from "geos_c.h":
+    # Types
+    ctypedef void *GEOSContextHandle_t
+    ctypedef struct GEOSGeometry
+    ctypedef struct GEOSCoordSequence
+    ctypedef void (*GEOSMessageHandler_r)(const char *message, void *userdata)
+
+    # GEOS Context & Messaging
+    GEOSContextHandle_t GEOS_init_r() nogil
+    void GEOS_finish_r(GEOSContextHandle_t handle) nogil
+    void GEOSContext_setErrorMessageHandler_r(GEOSContextHandle_t handle, GEOSMessageHandler_r ef, void* userData) nogil
+    void GEOSContext_setNoticeMessageHandler_r(GEOSContextHandle_t handle, GEOSMessageHandler_r nf, void* userData) nogil
+
+    # Geometry functions
+    const GEOSGeometry* GEOSGetGeometryN_r(GEOSContextHandle_t handle, const GEOSGeometry* g, int n) nogil
+    const GEOSGeometry* GEOSGetExteriorRing_r(GEOSContextHandle_t handle, const GEOSGeometry* g) nogil
+    const GEOSGeometry* GEOSGetInteriorRingN_r(GEOSContextHandle_t handle, const GEOSGeometry* g, int n) nogil
+    int GEOSGeomTypeId_r(GEOSContextHandle_t handle, GEOSGeometry* g) nogil
+
+    # Geometry creation / destruction
+    GEOSGeometry* GEOSGeom_clone_r(GEOSContextHandle_t handle, const GEOSGeometry* g) nogil
+    GEOSGeometry* GEOSGeom_createPoint_r(GEOSContextHandle_t handle, GEOSCoordSequence* s) nogil
+    GEOSGeometry* GEOSGeom_createLineString_r(GEOSContextHandle_t handle, GEOSCoordSequence* s) nogil
+    GEOSGeometry* GEOSGeom_createLinearRing_r(GEOSContextHandle_t handle, GEOSCoordSequence* s) nogil
+    GEOSGeometry* GEOSGeom_createEmptyPolygon_r(GEOSContextHandle_t handle) nogil
+    GEOSGeometry* GEOSGeom_createPolygon_r(GEOSContextHandle_t handle, GEOSGeometry* shell, GEOSGeometry** holes, unsigned int nholes) nogil
+    GEOSGeometry* GEOSGeom_createCollection_r(GEOSContextHandle_t handle, int type, GEOSGeometry** geoms, unsigned int ngeoms) nogil
+    void GEOSGeom_destroy_r(GEOSContextHandle_t handle, GEOSGeometry* g) nogil
+
+    # Coordinate sequences
+    GEOSCoordSequence* GEOSCoordSeq_create_r(GEOSContextHandle_t handle, unsigned int size, unsigned int dims) nogil
+    void GEOSCoordSeq_destroy_r(GEOSContextHandle_t handle, GEOSCoordSequence* s) nogil
+    int GEOSCoordSeq_setX_r(GEOSContextHandle_t handle, GEOSCoordSequence* s, unsigned int idx, double val) nogil
+    int GEOSCoordSeq_setY_r(GEOSContextHandle_t handle, GEOSCoordSequence* s, unsigned int idx, double val) nogil
+    int GEOSCoordSeq_setZ_r(GEOSContextHandle_t handle, GEOSCoordSequence* s, unsigned int idx, double val) nogil
+
+
+cdef class get_geos_handle:
+    cdef GEOSContextHandle_t handle
+    cdef char* last_error
+    cdef char* last_warning
+    cdef GEOSContextHandle_t __enter__(self)

billinglayer/python/shapely/_pygeos_api.pxd

@@ -0,0 +1,35 @@
+"""
+Provides a wrapper for the shapely.lib C API for use in Cython.
+Internally, the shapely C extension uses a PyCapsule to provide run-time access
+to function pointers within the C API.
+
+To use these functions, you must first call the following function in each Cython module:
+`import_shapely_c_api()`
+
+This uses a macro to dynamically load the functions from pointers in the PyCapsule.
+Each C function in shapely.lib exposed in the C API must be specially-wrapped to enable
+this capability.
+
+Segfaults will occur if the C API is not imported properly.
+"""
+
+cimport numpy as np
+from cpython.ref cimport PyObject
+
+from shapely._geos cimport GEOSContextHandle_t, GEOSCoordSequence, GEOSGeometry
+
+
+cdef extern from "c_api.h":
+    # shapely.lib C API loader; returns -1 on error
+    # MUST be called before calling other C API functions
+    int import_shapely_c_api() except -1
+
+    # C functions provided by the shapely.lib C API
+    # Note: GeometryObjects are always managed as Python objects
+    # in Cython to avoid memory leaks, not PyObject* (even though
+    # they are declared that way in the header file).
+    object PyGEOS_CreateGeometry(GEOSGeometry *ptr, GEOSContextHandle_t ctx)
+    char PyGEOS_GetGEOSGeometry(PyObject *obj, GEOSGeometry **out) nogil
+    GEOSCoordSequence* PyGEOS_CoordSeq_FromBuffer(GEOSContextHandle_t ctx, const double* buf,
+                                                 unsigned int size, unsigned int dims,
+                                                 char ring_closure) nogil

billinglayer/python/shapely/_ragged_array.py

@@ -0,0 +1,452 @@
+"""
+This modules provides a conversion to / from a ragged (or "jagged") array
+representation of the geometries.
+
+A ragged array is an irregular array of arrays of which each element can have
+a different length. As a result, such an array cannot be represented as a
+standard, rectangular nD array.
+The coordinates of geometries can be represented as arrays of arrays of
+coordinate pairs (possibly multiple levels of nesting, depending on the
+geometry type).
+
+
+Geometries, as a ragged array of coordinates, can be efficiently represented
+as contiguous arrays of coordinates provided that there is another data
+structure that keeps track of which range of coordinate values corresponds
+to a given geometry. This can be done using offsets, counts, or indices.
+
+This module currently implements offsets into the coordinates array. This
+is the ragged array representation defined by the the Apache Arrow project
+as "variable size list array" (https://arrow.apache.org/docs/format/Columnar.html#variable-size-list-layout).
+See for example https://cfconventions.org/Data/cf-conventions/cf-conventions-1.9/cf-conventions.html#representations-features
+for different options.
+
+The exact usage of the Arrow list array with varying degrees of nesting for the
+different geometry types is defined by the GeoArrow project:
+https://github.com/geoarrow/geoarrow
+
+"""
+import numpy as np
+
+from . import creation
+from ._geometry import (
+    GeometryType,
+    get_coordinate_dimension,
+    get_parts,
+    get_rings,
+    get_type_id,
+)
+from .coordinates import get_coordinates
+from .predicates import is_empty
+
+__all__ = ["to_ragged_array", "from_ragged_array"]
+
+
+# # GEOS -> coords/offset arrays (to_ragged_array)
+
+
+def _get_arrays_point(arr, include_z):
+    # only one array of coordinates
+    coords = get_coordinates(arr, include_z=include_z)
+
+    # empty points are represented by NaNs
+    empties = is_empty(arr)
+    if empties.any():
+        indices = np.nonzero(empties)[0]
+        indices = indices - np.arange(len(indices))
+        coords = np.insert(coords, indices, np.nan, axis=0)
+
+    return coords, ()
+
+
+def _indices_to_offsets(indices, n):
+    offsets = np.insert(np.bincount(indices).cumsum(), 0, 0)
+    if len(offsets) != n + 1:
+        # last geometries might be empty or missing
+        offsets = np.pad(
+            offsets,
+            (0, n + 1 - len(offsets)),
+            "constant",
+            constant_values=offsets[-1],
+        )
+    return offsets
+
+
+def _get_arrays_multipoint(arr, include_z):
+    # explode/flatten the MultiPoints
+    _, part_indices = get_parts(arr, return_index=True)
+    # the offsets into the multipoint parts
+    offsets = _indices_to_offsets(part_indices, len(arr))
+
+    # only one array of coordinates
+    coords = get_coordinates(arr, include_z=include_z)
+
+    return coords, (offsets,)
+
+
+def _get_arrays_linestring(arr, include_z):
+    # the coords and offsets into the coordinates of the linestrings
+    coords, indices = get_coordinates(arr, return_index=True, include_z=include_z)
+    offsets = _indices_to_offsets(indices, len(arr))
+
+    return coords, (offsets,)
+
+
+def _get_arrays_multilinestring(arr, include_z):
+    # explode/flatten the MultiLineStrings
+    arr_flat, part_indices = get_parts(arr, return_index=True)
+    # the offsets into the multilinestring parts
+    offsets2 = _indices_to_offsets(part_indices, len(arr))
+
+    # the coords and offsets into the coordinates of the linestrings
+    coords, indices = get_coordinates(arr_flat, return_index=True, include_z=include_z)
+    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)
+
+    return coords, (offsets1, offsets2)
+
+
+def _get_arrays_polygon(arr, include_z):
+    # explode/flatten the Polygons into Rings
+    arr_flat, ring_indices = get_rings(arr, return_index=True)
+    # the offsets into the exterior/interior rings of the multipolygon parts
+    offsets2 = _indices_to_offsets(ring_indices, len(arr))
+
+    # the coords and offsets into the coordinates of the rings
+    coords, indices = get_coordinates(arr_flat, return_index=True, include_z=include_z)
+    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)
+
+    return coords, (offsets1, offsets2)
+
+
+def _get_arrays_multipolygon(arr, include_z):
+    # explode/flatten the MultiPolygons
+    arr_flat, part_indices = get_parts(arr, return_index=True)
+    # the offsets into the multipolygon parts
+    offsets3 = _indices_to_offsets(part_indices, len(arr))
+
+    # explode/flatten the Polygons into Rings
+    arr_flat2, ring_indices = get_rings(arr_flat, return_index=True)
+    # the offsets into the exterior/interior rings of the multipolygon parts
+    offsets2 = np.insert(np.bincount(ring_indices).cumsum(), 0, 0)
+
+    # the coords and offsets into the coordinates of the rings
+    coords, indices = get_coordinates(arr_flat2, return_index=True, include_z=include_z)
+    offsets1 = np.insert(np.bincount(indices).cumsum(), 0, 0)
+
+    return coords, (offsets1, offsets2, offsets3)
+
+
+def to_ragged_array(geometries, include_z=None):
+    """
+    Converts geometries to a ragged array representation using a contiguous
+    array of coordinates and offset arrays.
+
+    This function converts an array of geometries to a ragged array
+    (i.e. irregular array of arrays) of coordinates, represented in memory
+    using a single contiguous array of the coordinates, and
+    up to 3 offset arrays that keep track where each sub-array
+    starts and ends.
+
+    This follows the in-memory layout of the variable size list arrays defined
+    by Apache Arrow, as specified for geometries by the GeoArrow project:
+    https://github.com/geoarrow/geoarrow.
+
+    Parameters
+    ----------
+    geometries : array_like
+        Array of geometries (1-dimensional).
+    include_z : bool, default None
+        If False, return 2D geometries. If True, include the third dimension
+        in the output (if a geometry has no third dimension, the z-coordinates
+        will be NaN). By default, will infer the dimensionality from the
+        input geometries. Note that this inference can be unreliable with
+        empty geometries (for a guaranteed result, it is recommended to
+        specify the keyword).
+
+    Returns
+    -------
+    tuple of (geometry_type, coords, offsets)
+        geometry_type : GeometryType
+            The type of the input geometries (required information for
+            roundtrip).
+        coords : np.ndarray
+            Contiguous array of shape (n, 2) or (n, 3) of all coordinates
+            of all input geometries.
+        offsets: tuple of np.ndarray
+            Offset arrays that make it possible to reconstruct the
+            geometries from the flat coordinates array. The number of
+            offset arrays depends on the geometry type. See
+            https://github.com/geoarrow/geoarrow/blob/main/format.md
+            for details.
+
+    Notes
+    -----
+    Mixed singular and multi geometry types of the same basic type are
+    allowed (e.g., Point and MultiPoint) and all singular types will be
+    treated as multi types.
+    GeometryCollections and other mixed geometry types are not supported.
+
+    See also
+    --------
+    from_ragged_array
+
+    Examples
+    --------
+    Consider a Polygon with one hole (interior ring):
+
+    >>> import shapely
+    >>> polygon = shapely.Polygon(
+    ...     [(0, 0), (10, 0), (10, 10), (0, 10)],
+    ...     holes=[[(2, 2), (3, 2), (2, 3)]]
+    ... )
+    >>> polygon
+    <POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (2 2, 3 2, 2 3, 2 2))>
+
+    This polygon can be thought of as a list of rings (first ring is the
+    exterior ring, subsequent rings are the interior rings), and each ring
+    as a list of coordinate pairs. This is very similar to how GeoJSON
+    represents the coordinates:
+
+    >>> import json
+    >>> json.loads(shapely.to_geojson(polygon))["coordinates"]
+    [[[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]],
+     [[2.0, 2.0], [3.0, 2.0], [2.0, 3.0], [2.0, 2.0]]]
+
+    This function will return a similar list of lists of lists, but
+    using a single contiguous array of coordinates, and multiple arrays of
+    offsets:
+
+    >>> geometry_type, coords, offsets = shapely.to_ragged_array([polygon])
+    >>> geometry_type
+    <GeometryType.POLYGON: 3>
+    >>> coords
+    array([[ 0.,  0.],
+           [10.,  0.],
+           [10., 10.],
+           [ 0., 10.],
+           [ 0.,  0.],
+           [ 2.,  2.],
+           [ 3.,  2.],
+           [ 2.,  3.],
+           [ 2.,  2.]])
+
+    >>> offsets
+    (array([0, 5, 9]), array([0, 2]))
+
+    As an example how to interpret the offsets: the i-th ring in the
+    coordinates is represented by ``offsets[0][i]`` to ``offsets[0][i+1]``:
+
+    >>> exterior_ring_start, exterior_ring_end = offsets[0][0], offsets[0][1]
+    >>> coords[exterior_ring_start:exterior_ring_end]
+    array([[ 0.,  0.],
+           [10.,  0.],
+           [10., 10.],
+           [ 0., 10.],
+           [ 0.,  0.]])
+
+    """
+    geometries = np.asarray(geometries)
+    if include_z is None:
+        include_z = np.any(
+            get_coordinate_dimension(geometries[~is_empty(geometries)]) == 3
+        )
+
+    geom_types = np.unique(get_type_id(geometries))
+    # ignore missing values (type of -1)
+    geom_types = geom_types[geom_types >= 0]
+
+    if len(geom_types) == 1:
+        typ = GeometryType(geom_types[0])
+        if typ == GeometryType.POINT:
+            coords, offsets = _get_arrays_point(geometries, include_z)
+        elif typ == GeometryType.LINESTRING:
+            coords, offsets = _get_arrays_linestring(geometries, include_z)
+        elif typ == GeometryType.POLYGON:
+            coords, offsets = _get_arrays_polygon(geometries, include_z)
+        elif typ == GeometryType.MULTIPOINT:
+            coords, offsets = _get_arrays_multipoint(geometries, include_z)
+        elif typ == GeometryType.MULTILINESTRING:
+            coords, offsets = _get_arrays_multilinestring(geometries, include_z)
+        elif typ == GeometryType.MULTIPOLYGON:
+            coords, offsets = _get_arrays_multipolygon(geometries, include_z)
+        else:
+            raise ValueError(f"Geometry type {typ.name} is not supported")
+
+    elif len(geom_types) == 2:
+        if set(geom_types) == {GeometryType.POINT, GeometryType.MULTIPOINT}:
+            typ = GeometryType.MULTIPOINT
+            coords, offsets = _get_arrays_multipoint(geometries, include_z)
+        elif set(geom_types) == {GeometryType.LINESTRING, GeometryType.MULTILINESTRING}:
+            typ = GeometryType.MULTILINESTRING
+            coords, offsets = _get_arrays_multilinestring(geometries, include_z)
+        elif set(geom_types) == {GeometryType.POLYGON, GeometryType.MULTIPOLYGON}:
+            typ = GeometryType.MULTIPOLYGON
+            coords, offsets = _get_arrays_multipolygon(geometries, include_z)
+        else:
+            raise ValueError(
+                "Geometry type combination is not supported "
+                f"({[GeometryType(t).name for t in geom_types]})"
+            )
+    else:
+        raise ValueError(
+            "Geometry type combination is not supported "
+            f"({[GeometryType(t).name for t in geom_types]})"
+        )
+
+    return typ, coords, offsets
+
+
+# # coords/offset arrays -> GEOS (from_ragged_array)
+
+
+def _point_from_flatcoords(coords):
+    result = creation.points(coords)
+
+    # Older versions of GEOS (<= 3.9) don't automatically convert NaNs
+    # to empty points -> do manually
+    empties = np.isnan(coords).all(axis=1)
+    if empties.any():
+        result[empties] = creation.empty(1, geom_type=GeometryType.POINT).item()
+
+    return result
+
+
+def _multipoint_from_flatcoords(coords, offsets):
+    # recreate points
+    points = creation.points(coords)
+
+    # recreate multipoints
+    multipoint_parts = np.diff(offsets)
+    multipoint_indices = np.repeat(np.arange(len(multipoint_parts)), multipoint_parts)
+
+    result = np.empty(len(offsets) - 1, dtype=object)
+    result = creation.multipoints(points, indices=multipoint_indices, out=result)
+    result[multipoint_parts == 0] = creation.empty(
+        1, geom_type=GeometryType.MULTIPOINT
+    ).item()
+
+    return result
+
+
+def _linestring_from_flatcoords(coords, offsets):
+    # recreate linestrings
+    linestring_n = np.diff(offsets)
+    linestring_indices = np.repeat(np.arange(len(linestring_n)), linestring_n)
+
+    result = np.empty(len(offsets) - 1, dtype=object)
+    result = creation.linestrings(coords, indices=linestring_indices, out=result)
+    result[linestring_n == 0] = creation.empty(
+        1, geom_type=GeometryType.LINESTRING
+    ).item()
+    return result
+
+
+def _multilinestrings_from_flatcoords(coords, offsets1, offsets2):
+    # recreate linestrings
+    linestrings = _linestring_from_flatcoords(coords, offsets1)
+
+    # recreate multilinestrings
+    multilinestring_parts = np.diff(offsets2)
+    multilinestring_indices = np.repeat(
+        np.arange(len(multilinestring_parts)), multilinestring_parts
+    )
+
+    result = np.empty(len(offsets2) - 1, dtype=object)
+    result = creation.multilinestrings(
+        linestrings, indices=multilinestring_indices, out=result
+    )
+    result[multilinestring_parts == 0] = creation.empty(
+        1, geom_type=GeometryType.MULTILINESTRING
+    ).item()
+
+    return result
+
+
+def _polygon_from_flatcoords(coords, offsets1, offsets2):
+    # recreate rings
+    ring_lengths = np.diff(offsets1)
+    ring_indices = np.repeat(np.arange(len(ring_lengths)), ring_lengths)
+    rings = creation.linearrings(coords, indices=ring_indices)
+
+    # recreate polygons
+    polygon_rings_n = np.diff(offsets2)
+    polygon_indices = np.repeat(np.arange(len(polygon_rings_n)), polygon_rings_n)
+    result = np.empty(len(offsets2) - 1, dtype=object)
+    result = creation.polygons(rings, indices=polygon_indices, out=result)
+    result[polygon_rings_n == 0] = creation.empty(
+        1, geom_type=GeometryType.POLYGON
+    ).item()
+
+    return result
+
+
+def _multipolygons_from_flatcoords(coords, offsets1, offsets2, offsets3):
+    # recreate polygons
+    polygons = _polygon_from_flatcoords(coords, offsets1, offsets2)
+
+    # recreate multipolygons
+    multipolygon_parts = np.diff(offsets3)
+    multipolygon_indices = np.repeat(
+        np.arange(len(multipolygon_parts)), multipolygon_parts
+    )
+    result = np.empty(len(offsets3) - 1, dtype=object)
+    result = creation.multipolygons(polygons, indices=multipolygon_indices, out=result)
+    result[multipolygon_parts == 0] = creation.empty(
+        1, geom_type=GeometryType.MULTIPOLYGON
+    ).item()
+
+    return result
+
+
+def from_ragged_array(geometry_type, coords, offsets=None):
+    """
+    Creates geometries from a contiguous array of coordinates
+    and offset arrays.
+
+    This function creates geometries from the ragged array representation
+    as returned by ``to_ragged_array``.
+
+    This follows the in-memory layout of the variable size list arrays defined
+    by Apache Arrow, as specified for geometries by the GeoArrow project:
+    https://github.com/geoarrow/geoarrow.
+
+    See :func:`to_ragged_array` for more details.
+
+    Parameters
+    ----------
+    geometry_type : GeometryType
+        The type of geometry to create.
+    coords : np.ndarray
+        Contiguous array of shape (n, 2) or (n, 3) of all coordinates
+        for the geometries.
+    offsets: tuple of np.ndarray
+        Offset arrays that allow to reconstruct the geometries based on the
+        flat coordinates array. The number of offset arrays depends on the
+        geometry type. See
+        https://github.com/geoarrow/geoarrow/blob/main/format.md for details.
+
+    Returns
+    -------
+    np.ndarray
+        Array of geometries (1-dimensional).
+
+    See Also
+    --------
+    to_ragged_array
+
+    """
+    if geometry_type == GeometryType.POINT:
+        assert offsets is None or len(offsets) == 0
+        return _point_from_flatcoords(coords)
+    if geometry_type == GeometryType.LINESTRING:
+        return _linestring_from_flatcoords(coords, *offsets)
+    if geometry_type == GeometryType.POLYGON:
+        return _polygon_from_flatcoords(coords, *offsets)
+    elif geometry_type == GeometryType.MULTIPOINT:
+        return _multipoint_from_flatcoords(coords, *offsets)
+    elif geometry_type == GeometryType.MULTILINESTRING:
+        return _multilinestrings_from_flatcoords(coords, *offsets)
+    elif geometry_type == GeometryType.MULTIPOLYGON:
+        return _multipolygons_from_flatcoords(coords, *offsets)
+    else:
+        raise ValueError(f"Geometry type {geometry_type.name} is not supported")

billinglayer/python/shapely/_version.py

@@ -0,0 +1,21 @@
+
+# This file was generated by 'versioneer.py' (0.28) from
+# revision-control system data, or from the parent directory name of an
+# unpacked source archive. Distribution tarballs contain a pre-generated copy
+# of this file.
+
+import json
+
+version_json = '''
+{
+ "date": "2023-01-30T08:49:38+0100",
+ "dirty": false,
+ "error": null,
+ "full-revisionid": "5f0352897968e82283a61fe6b201df87993628fb",
+ "version": "2.0.1"
+}
+'''  # END VERSION_JSON
+
+
+def get_versions():
+    return json.loads(version_json)

billinglayer/python/shapely/affinity.py

@@ -0,0 +1,250 @@
+"""Affine transforms, both in general and specific, named transforms."""
+
+from math import cos, pi, sin, tan
+
+import numpy as np
+
+import shapely
+
+__all__ = ["affine_transform", "rotate", "scale", "skew", "translate"]
+
+
+def affine_transform(geom, matrix):
+    r"""Return a transformed geometry using an affine transformation matrix.
+
+    The coefficient matrix is provided as a list or tuple with 6 or 12 items
+    for 2D or 3D transformations, respectively.
+
+    For 2D affine transformations, the 6 parameter matrix is::
+
+        [a, b, d, e, xoff, yoff]
+
+    which represents the augmented matrix::
+
+        [x']   / a  b xoff \ [x]
+        [y'] = | d  e yoff | [y]
+        [1 ]   \ 0  0   1  / [1]
+
+    or the equations for the transformed coordinates::
+
+        x' = a * x + b * y + xoff
+        y' = d * x + e * y + yoff
+
+    For 3D affine transformations, the 12 parameter matrix is::
+
+        [a, b, c, d, e, f, g, h, i, xoff, yoff, zoff]
+
+    which represents the augmented matrix::
+
+        [x']   / a  b  c xoff \ [x]
+        [y'] = | d  e  f yoff | [y]
+        [z']   | g  h  i zoff | [z]
+        [1 ]   \ 0  0  0   1  / [1]
+
+    or the equations for the transformed coordinates::
+
+        x' = a * x + b * y + c * z + xoff
+        y' = d * x + e * y + f * z + yoff
+        z' = g * x + h * y + i * z + zoff
+    """
+    if len(matrix) == 6:
+        ndim = 2
+        a, b, d, e, xoff, yoff = matrix
+        if geom.has_z:
+            ndim = 3
+            i = 1.0
+            c = f = g = h = zoff = 0.0
+    elif len(matrix) == 12:
+        ndim = 3
+        a, b, c, d, e, f, g, h, i, xoff, yoff, zoff = matrix
+        if not geom.has_z:
+            ndim = 2
+    else:
+        raise ValueError("'matrix' expects either 6 or 12 coefficients")
+    if ndim == 2:
+        A = np.array([[a, b], [d, e]], dtype=float)
+        off = np.array([xoff, yoff], dtype=float)
+    else:
+        A = np.array([[a, b, c], [d, e, f], [g, h, i]], dtype=float)
+        off = np.array([xoff, yoff, zoff], dtype=float)
+
+    def _affine_coords(coords):
+        return np.matmul(A, coords.T).T + off
+
+    return shapely.transform(geom, _affine_coords, include_z=ndim == 3)
+
+
+def interpret_origin(geom, origin, ndim):
+    """Returns interpreted coordinate tuple for origin parameter.
+
+    This is a helper function for other transform functions.
+
+    The point of origin can be a keyword 'center' for the 2D bounding box
+    center, 'centroid' for the geometry's 2D centroid, a Point object or a
+    coordinate tuple (x0, y0, z0).
+    """
+    # get coordinate tuple from 'origin' from keyword or Point type
+    if origin == "center":
+        # bounding box center
+        minx, miny, maxx, maxy = geom.bounds
+        origin = ((maxx + minx) / 2.0, (maxy + miny) / 2.0)
+    elif origin == "centroid":
+        origin = geom.centroid.coords[0]
+    elif isinstance(origin, str):
+        raise ValueError(f"'origin' keyword {origin!r} is not recognized")
+    elif getattr(origin, "geom_type", None) == "Point":
+        origin = origin.coords[0]
+
+    # origin should now be tuple-like
+    if len(origin) not in (2, 3):
+        raise ValueError("Expected number of items in 'origin' to be " "either 2 or 3")
+    if ndim == 2:
+        return origin[0:2]
+    else:  # 3D coordinate
+        if len(origin) == 2:
+            return origin + (0.0,)
+        else:
+            return origin
+
+
+def rotate(geom, angle, origin="center", use_radians=False):
+    r"""Returns a rotated geometry on a 2D plane.
+
+    The angle of rotation can be specified in either degrees (default) or
+    radians by setting ``use_radians=True``. Positive angles are
+    counter-clockwise and negative are clockwise rotations.
+
+    The point of origin can be a keyword 'center' for the bounding box
+    center (default), 'centroid' for the geometry's centroid, a Point object
+    or a coordinate tuple (x0, y0).
+
+    The affine transformation matrix for 2D rotation is:
+
+      / cos(r) -sin(r) xoff \
+      | sin(r)  cos(r) yoff |
+      \   0       0      1  /
+
+    where the offsets are calculated from the origin Point(x0, y0):
+
+        xoff = x0 - x0 * cos(r) + y0 * sin(r)
+        yoff = y0 - x0 * sin(r) - y0 * cos(r)
+    """
+    if geom.is_empty:
+        return geom
+    if not use_radians:  # convert from degrees
+        angle = angle * pi / 180.0
+    cosp = cos(angle)
+    sinp = sin(angle)
+    if abs(cosp) < 2.5e-16:
+        cosp = 0.0
+    if abs(sinp) < 2.5e-16:
+        sinp = 0.0
+    x0, y0 = interpret_origin(geom, origin, 2)
+
+    # fmt: off
+    matrix = (cosp, -sinp, 0.0,
+              sinp, cosp, 0.0,
+              0.0, 0.0, 1.0,
+              x0 - x0 * cosp + y0 * sinp, y0 - x0 * sinp - y0 * cosp, 0.0)
+    # fmt: on
+    return affine_transform(geom, matrix)
+
+
+def scale(geom, xfact=1.0, yfact=1.0, zfact=1.0, origin="center"):
+    r"""Returns a scaled geometry, scaled by factors along each dimension.
+
+    The point of origin can be a keyword 'center' for the 2D bounding box
+    center (default), 'centroid' for the geometry's 2D centroid, a Point
+    object or a coordinate tuple (x0, y0, z0).
+
+    Negative scale factors will mirror or reflect coordinates.
+
+    The general 3D affine transformation matrix for scaling is:
+
+        / xfact  0    0   xoff \
+        |   0  yfact  0   yoff |
+        |   0    0  zfact zoff |
+        \   0    0    0     1  /
+
+    where the offsets are calculated from the origin Point(x0, y0, z0):
+
+        xoff = x0 - x0 * xfact
+        yoff = y0 - y0 * yfact
+        zoff = z0 - z0 * zfact
+    """
+    if geom.is_empty:
+        return geom
+    x0, y0, z0 = interpret_origin(geom, origin, 3)
+
+    # fmt: off
+    matrix = (xfact, 0.0, 0.0,
+              0.0, yfact, 0.0,
+              0.0, 0.0, zfact,
+              x0 - x0 * xfact, y0 - y0 * yfact, z0 - z0 * zfact)
+    # fmt: on
+    return affine_transform(geom, matrix)
+
+
+def skew(geom, xs=0.0, ys=0.0, origin="center", use_radians=False):
+    r"""Returns a skewed geometry, sheared by angles along x and y dimensions.
+
+    The shear angle can be specified in either degrees (default) or radians
+    by setting ``use_radians=True``.
+
+    The point of origin can be a keyword 'center' for the bounding box
+    center (default), 'centroid' for the geometry's centroid, a Point object
+    or a coordinate tuple (x0, y0).
+
+    The general 2D affine transformation matrix for skewing is:
+
+        /   1    tan(xs) xoff \
+        | tan(ys)  1     yoff |
+        \   0      0       1  /
+
+    where the offsets are calculated from the origin Point(x0, y0):
+
+        xoff = -y0 * tan(xs)
+        yoff = -x0 * tan(ys)
+    """
+    if geom.is_empty:
+        return geom
+    if not use_radians:  # convert from degrees
+        xs = xs * pi / 180.0
+        ys = ys * pi / 180.0
+    tanx = tan(xs)
+    tany = tan(ys)
+    if abs(tanx) < 2.5e-16:
+        tanx = 0.0
+    if abs(tany) < 2.5e-16:
+        tany = 0.0
+    x0, y0 = interpret_origin(geom, origin, 2)
+
+    # fmt: off
+    matrix = (1.0, tanx, 0.0,
+              tany, 1.0, 0.0,
+              0.0, 0.0, 1.0,
+              -y0 * tanx, -x0 * tany, 0.0)
+    # fmt: on
+    return affine_transform(geom, matrix)
+
+
+def translate(geom, xoff=0.0, yoff=0.0, zoff=0.0):
+    r"""Returns a translated geometry shifted by offsets along each dimension.
+
+    The general 3D affine transformation matrix for translation is:
+
+        / 1  0  0 xoff \
+        | 0  1  0 yoff |
+        | 0  0  1 zoff |
+        \ 0  0  0   1  /
+    """
+    if geom.is_empty:
+        return geom
+
+    # fmt: off
+    matrix = (1.0, 0.0, 0.0,
+              0.0, 1.0, 0.0,
+              0.0, 0.0, 1.0,
+              xoff, yoff, zoff)
+    # fmt: on
+    return affine_transform(geom, matrix)

billinglayer/python/shapely/algorithms/cga.py

@@ -0,0 +1,24 @@
+import numpy as np
+
+import shapely
+
+
+def signed_area(ring):
+    """Return the signed area enclosed by a ring in linear time using the
+    algorithm at: https://web.archive.org/web/20080209143651/http://cgafaq.info:80/wiki/Polygon_Area
+    """
+    coords = np.array(ring.coords)[:, :2]
+    xs, ys = np.vstack([coords, coords[1]]).T
+    return np.sum(xs[1:-1] * (ys[2:] - ys[:-2])) / 2.0
+
+
+def is_ccw_impl(name=None):
+    """Predicate implementation"""
+
+    def is_ccw_op(ring):
+        return signed_area(ring) >= 0.0
+
+    if shapely.geos_version >= (3, 7, 0):
+        return shapely.is_ccw
+    else:
+        return is_ccw_op

billinglayer/python/shapely/algorithms/polylabel.py

@@ -0,0 +1,139 @@
+from heapq import heappop, heappush
+
+from ..errors import TopologicalError
+from ..geometry import Point
+
+
+class Cell:
+    """A `Cell`'s centroid property is a potential solution to finding the pole
+    of inaccessibility for a given polygon. Rich comparison operators are used
+    for sorting `Cell` objects in a priority queue based on the potential
+    maximum distance of any theoretical point within a cell to a given
+    polygon's exterior boundary.
+    """
+
+    def __init__(self, x, y, h, polygon):
+        self.x = x
+        self.y = y
+        self.h = h  # half of cell size
+        self.centroid = Point(x, y)  # cell centroid, potential solution
+
+        # distance from cell centroid to polygon exterior
+        self.distance = self._dist(polygon)
+
+        # max distance to polygon exterior within a cell
+        self.max_distance = self.distance + h * 1.4142135623730951  # sqrt(2)
+
+    # rich comparison operators for sorting in minimum priority queue
+    def __lt__(self, other):
+        return self.max_distance > other.max_distance
+
+    def __le__(self, other):
+        return self.max_distance >= other.max_distance
+
+    def __eq__(self, other):
+        return self.max_distance == other.max_distance
+
+    def __ne__(self, other):
+        return self.max_distance != other.max_distance
+
+    def __gt__(self, other):
+        return self.max_distance < other.max_distance
+
+    def __ge__(self, other):
+        return self.max_distance <= other.max_distance
+
+    def _dist(self, polygon):
+        """Signed distance from Cell centroid to polygon outline. The returned
+        value is negative if the point is outside of the polygon exterior
+        boundary.
+        """
+        inside = polygon.contains(self.centroid)
+        distance = self.centroid.distance(polygon.exterior)
+        for interior in polygon.interiors:
+            distance = min(distance, self.centroid.distance(interior))
+        if inside:
+            return distance
+        return -distance
+
+
+def polylabel(polygon, tolerance=1.0):
+    """Finds pole of inaccessibility for a given polygon. Based on
+    Vladimir Agafonkin's https://github.com/mapbox/polylabel
+
+    Parameters
+    ----------
+    polygon : shapely.geometry.Polygon
+    tolerance : int or float, optional
+                `tolerance` represents the highest resolution in units of the
+                input geometry that will be considered for a solution. (default
+                value is 1.0).
+
+    Returns
+    -------
+    shapely.geometry.Point
+        A point representing the pole of inaccessibility for the given input
+        polygon.
+
+    Raises
+    ------
+    shapely.errors.TopologicalError
+        If the input polygon is not a valid geometry.
+
+    Example
+    -------
+    >>> from shapely import LineString
+    >>> polygon = LineString([(0, 0), (50, 200), (100, 100), (20, 50),
+    ... (-100, -20), (-150, -200)]).buffer(100)
+    >>> polylabel(polygon, tolerance=10).wkt
+    'POINT (59.35615556364569 121.83919629746435)'
+    """
+    if not polygon.is_valid:
+        raise TopologicalError("Invalid polygon")
+    minx, miny, maxx, maxy = polygon.bounds
+    width = maxx - minx
+    height = maxy - miny
+    cell_size = min(width, height)
+    h = cell_size / 2.0
+    cell_queue = []
+
+    # First best cell approximation is one constructed from the centroid
+    # of the polygon
+    x, y = polygon.centroid.coords[0]
+    best_cell = Cell(x, y, 0, polygon)
+
+    # Special case for rectangular polygons avoiding floating point error
+    bbox_cell = Cell(minx + width / 2.0, miny + height / 2, 0, polygon)
+    if bbox_cell.distance > best_cell.distance:
+        best_cell = bbox_cell
+
+    # build a regular square grid covering the polygon
+    x = minx
+    while x < maxx:
+        y = miny
+        while y < maxy:
+            heappush(cell_queue, Cell(x + h, y + h, h, polygon))
+            y += cell_size
+        x += cell_size
+
+    # minimum priority queue
+    while cell_queue:
+        cell = heappop(cell_queue)
+
+        # update the best cell if we find a better one
+        if cell.distance > best_cell.distance:
+            best_cell = cell
+
+        # continue to the next iteration if we can't find a better solution
+        # based on tolerance
+        if cell.max_distance - best_cell.distance <= tolerance:
+            continue
+
+        # split the cell into quadrants
+        h = cell.h / 2.0
+        heappush(cell_queue, Cell(cell.x - h, cell.y - h, h, polygon))
+        heappush(cell_queue, Cell(cell.x + h, cell.y - h, h, polygon))
+        heappush(cell_queue, Cell(cell.x - h, cell.y + h, h, polygon))
+        heappush(cell_queue, Cell(cell.x + h, cell.y + h, h, polygon))
+
+    return best_cell.centroid

billinglayer/python/shapely/coordinates.py

@@ -0,0 +1,198 @@
+import numpy as np
+
+from . import lib
+
+__all__ = ["transform", "count_coordinates", "get_coordinates", "set_coordinates"]
+
+
+def transform(geometry, transformation, include_z=False):
+    """Returns a copy of a geometry array with a function applied to its
+    coordinates.
+
+    With the default of ``include_z=False``, all returned geometries will be
+    two-dimensional; the third dimension will be discarded, if present.
+    When specifying ``include_z=True``, the returned geometries preserve
+    the dimensionality of the respective input geometries.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    transformation : function
+        A function that transforms a (N, 2) or (N, 3) ndarray of float64 to
+        another (N, 2) or (N, 3) ndarray of float64.
+    include_z : bool, default False
+        If True, include the third dimension in the coordinates array
+        that is passed to the ``transformation`` function. If a
+        geometry has no third dimension, the z-coordinates passed to the
+        function will be NaN.
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point
+    >>> transform(Point(0, 0), lambda x: x + 1)
+    <POINT (1 1)>
+    >>> transform(LineString([(2, 2), (4, 4)]), lambda x: x * [2, 3])
+    <LINESTRING (4 6, 8 12)>
+    >>> transform(None, lambda x: x) is None
+    True
+    >>> transform([Point(0, 0), None], lambda x: x).tolist()
+    [<POINT (0 0)>, None]
+
+    By default, the third dimension is ignored:
+
+    >>> transform(Point(0, 0, 0), lambda x: x + 1)
+    <POINT (1 1)>
+    >>> transform(Point(0, 0, 0), lambda x: x + 1, include_z=True)
+    <POINT Z (1 1 1)>
+    """
+    geometry_arr = np.array(geometry, dtype=np.object_)  # makes a copy
+    coordinates = lib.get_coordinates(geometry_arr, include_z, False)
+    new_coordinates = transformation(coordinates)
+    # check the array to yield understandable error messages
+    if not isinstance(new_coordinates, np.ndarray):
+        raise ValueError("The provided transformation did not return a numpy array")
+    if new_coordinates.dtype != np.float64:
+        raise ValueError(
+            "The provided transformation returned an array with an unexpected "
+            f"dtype ({new_coordinates.dtype})"
+        )
+    if new_coordinates.shape != coordinates.shape:
+        # if the shape is too small we will get a segfault
+        raise ValueError(
+            "The provided transformation returned an array with an unexpected "
+            f"shape ({new_coordinates.shape})"
+        )
+    geometry_arr = lib.set_coordinates(geometry_arr, new_coordinates)
+    if geometry_arr.ndim == 0 and not isinstance(geometry, np.ndarray):
+        return geometry_arr.item()
+    return geometry_arr
+
+
+def count_coordinates(geometry):
+    """Counts the number of coordinate pairs in a geometry array.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point
+    >>> count_coordinates(Point(0, 0))
+    1
+    >>> count_coordinates(LineString([(2, 2), (4, 2)]))
+    2
+    >>> count_coordinates(None)
+    0
+    >>> count_coordinates([Point(0, 0), None])
+    1
+    """
+    return lib.count_coordinates(np.asarray(geometry, dtype=np.object_))
+
+
+def get_coordinates(geometry, include_z=False, return_index=False):
+    """Gets coordinates from a geometry array as an array of floats.
+
+    The shape of the returned array is (N, 2), with N being the number of
+    coordinate pairs. With the default of ``include_z=False``, three-dimensional
+    data is ignored. When specifying ``include_z=True``, the shape of the
+    returned array is (N, 3).
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    include_z : bool, default False
+        If, True include the third dimension in the output. If a geometry
+        has no third dimension, the z-coordinates will be NaN.
+    return_index : bool, default False
+        If True, also return the index of each returned geometry as a separate
+        ndarray of integers. For multidimensional arrays, this indexes into the
+        flattened array (in C contiguous order).
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point
+    >>> get_coordinates(Point(0, 0)).tolist()
+    [[0.0, 0.0]]
+    >>> get_coordinates(LineString([(2, 2), (4, 4)])).tolist()
+    [[2.0, 2.0], [4.0, 4.0]]
+    >>> get_coordinates(None)
+    array([], shape=(0, 2), dtype=float64)
+
+    By default the third dimension is ignored:
+
+    >>> get_coordinates(Point(0, 0, 0)).tolist()
+    [[0.0, 0.0]]
+    >>> get_coordinates(Point(0, 0, 0), include_z=True).tolist()
+    [[0.0, 0.0, 0.0]]
+
+    When return_index=True, indexes are returned also:
+
+    >>> geometries = [LineString([(2, 2), (4, 4)]), Point(0, 0)]
+    >>> coordinates, index = get_coordinates(geometries, return_index=True)
+    >>> coordinates.tolist(), index.tolist()
+    ([[2.0, 2.0], [4.0, 4.0], [0.0, 0.0]], [0, 0, 1])
+    """
+    return lib.get_coordinates(
+        np.asarray(geometry, dtype=np.object_), include_z, return_index
+    )
+
+
+def set_coordinates(geometry, coordinates):
+    """Adapts the coordinates of a geometry array in-place.
+
+    If the coordinates array has shape (N, 2), all returned geometries
+    will be two-dimensional, and the third dimension will be discarded,
+    if present. If the coordinates array has shape (N, 3), the returned
+    geometries preserve the dimensionality of the input geometries.
+
+    .. warning::
+
+        The geometry array is modified in-place! If you do not want to
+        modify the original array, you can do
+        ``set_coordinates(arr.copy(), newcoords)``.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    coordinates: array_like
+
+    See Also
+    --------
+    transform : Returns a copy of a geometry array with a function applied to its
+        coordinates.
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point
+    >>> set_coordinates(Point(0, 0), [[1, 1]])
+    <POINT (1 1)>
+    >>> set_coordinates([Point(0, 0), LineString([(0, 0), (0, 0)])], [[1, 2], [3, 4], [5, 6]]).tolist()
+    [<POINT (1 2)>, <LINESTRING (3 4, 5 6)>]
+    >>> set_coordinates([None, Point(0, 0)], [[1, 2]]).tolist()
+    [None, <POINT (1 2)>]
+
+    Third dimension of input geometry is discarded if coordinates array does
+    not include one:
+
+    >>> set_coordinates(Point(0, 0, 0), [[1, 1]])
+    <POINT (1 1)>
+    >>> set_coordinates(Point(0, 0, 0), [[1, 1, 1]])
+    <POINT Z (1 1 1)>
+    """
+    geometry_arr = np.asarray(geometry, dtype=np.object_)
+    coordinates = np.atleast_2d(np.asarray(coordinates)).astype(np.float64)
+    if coordinates.ndim != 2:
+        raise ValueError(
+            "The coordinate array should have dimension of 2 "
+            f"(has {coordinates.ndim})"
+        )
+    n_coords = lib.count_coordinates(geometry_arr)
+    if (coordinates.shape[0] != n_coords) or (coordinates.shape[1] not in {2, 3}):
+        raise ValueError(
+            f"The coordinate array has an invalid shape {coordinates.shape}"
+        )
+    lib.set_coordinates(geometry_arr, coordinates)
+    if geometry_arr.ndim == 0 and not isinstance(geometry, np.ndarray):
+        return geometry_arr.item()
+    return geometry_arr

billinglayer/python/shapely/coords.py

@@ -0,0 +1,62 @@
+"""Coordinate sequence utilities
+"""
+from array import array
+
+
+class CoordinateSequence:
+    """
+    Iterative access to coordinate tuples from the parent geometry's coordinate
+    sequence.
+
+    Example:
+
+      >>> from shapely.wkt import loads
+      >>> g = loads('POINT (0.0 0.0)')
+      >>> list(g.coords)
+      [(0.0, 0.0)]
+
+    """
+
+    def __init__(self, coords):
+        self._coords = coords
+
+    def __len__(self):
+        return self._coords.shape[0]
+
+    def __iter__(self):
+        for i in range(self.__len__()):
+            yield tuple(self._coords[i].tolist())
+
+    def __getitem__(self, key):
+        m = self.__len__()
+        if isinstance(key, int):
+            if key + m < 0 or key >= m:
+                raise IndexError("index out of range")
+            if key < 0:
+                i = m + key
+            else:
+                i = key
+            return tuple(self._coords[i].tolist())
+        elif isinstance(key, slice):
+            res = []
+            start, stop, stride = key.indices(m)
+            for i in range(start, stop, stride):
+                res.append(tuple(self._coords[i].tolist()))
+            return res
+        else:
+            raise TypeError("key must be an index or slice")
+
+    def __array__(self, dtype=None):
+        return self._coords
+
+    @property
+    def xy(self):
+        """X and Y arrays"""
+        m = self.__len__()
+        x = array("d")
+        y = array("d")
+        for i in range(m):
+            xy = self._coords[i].tolist()
+            x.append(xy[0])
+            y.append(xy[1])
+        return x, y

billinglayer/python/shapely/creation.py

@@ -0,0 +1,565 @@
+import numpy as np
+
+from . import Geometry, GeometryType, lib
+from ._geometry_helpers import collections_1d, simple_geometries_1d
+from .decorators import multithreading_enabled
+from .io import from_wkt
+
+__all__ = [
+    "points",
+    "linestrings",
+    "linearrings",
+    "polygons",
+    "multipoints",
+    "multilinestrings",
+    "multipolygons",
+    "geometrycollections",
+    "box",
+    "prepare",
+    "destroy_prepared",
+    "empty",
+]
+
+
+def _xyz_to_coords(x, y, z):
+    if y is None:
+        return x
+    if z is None:
+        coords = np.broadcast_arrays(x, y)
+    else:
+        coords = np.broadcast_arrays(x, y, z)
+    return np.stack(coords, axis=-1)
+
+
+@multithreading_enabled
+def points(coords, y=None, z=None, indices=None, out=None, **kwargs):
+    """Create an array of points.
+
+    Parameters
+    ----------
+    coords : array_like
+        An array of coordinate tuples (2- or 3-dimensional) or, if ``y`` is
+        provided, an array of x coordinates.
+    y : array_like, optional
+    z : array_like, optional
+    indices : array_like, optional
+        Indices into the target array where input coordinates belong. If
+        provided, the coords should be 2D with shape (N, 2) or (N, 3) and
+        indices should be an array of shape (N,) with integers in increasing
+        order. Missing indices result in a ValueError unless ``out`` is
+        provided, in which case the original value in ``out`` is kept.
+    out : ndarray, optional
+        An array (with dtype object) to output the geometries into.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+        Ignored if ``indices`` is provided.
+
+    Examples
+    --------
+    >>> points([[0, 1], [4, 5]]).tolist()
+    [<POINT (0 1)>, <POINT (4 5)>]
+    >>> points([0, 1, 2])
+    <POINT Z (0 1 2)>
+
+    Notes
+    -----
+
+    - GEOS >=3.10 automatically converts POINT (nan nan) to POINT EMPTY.
+    - Usage of the ``y`` and ``z`` arguments will prevents lazy evaluation in ``dask``.
+      Instead provide the coordinates as an array with shape ``(..., 2)`` or ``(..., 3)`` using only the ``coords`` argument.
+    """
+    coords = _xyz_to_coords(coords, y, z)
+    if indices is None:
+        return lib.points(coords, out=out, **kwargs)
+    else:
+        return simple_geometries_1d(coords, indices, GeometryType.POINT, out=out)
+
+
+@multithreading_enabled
+def linestrings(coords, y=None, z=None, indices=None, out=None, **kwargs):
+    """Create an array of linestrings.
+
+    This function will raise an exception if a linestring contains less than
+    two points.
+
+    Parameters
+    ----------
+    coords : array_like
+        An array of lists of coordinate tuples (2- or 3-dimensional) or, if ``y``
+        is provided, an array of lists of x coordinates
+    y : array_like, optional
+    z : array_like, optional
+    indices : array_like, optional
+        Indices into the target array where input coordinates belong. If
+        provided, the coords should be 2D with shape (N, 2) or (N, 3) and
+        indices should be an array of shape (N,) with integers in increasing
+        order. Missing indices result in a ValueError unless ``out`` is
+        provided, in which case the original value in ``out`` is kept.
+    out : ndarray, optional
+        An array (with dtype object) to output the geometries into.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+        Ignored if ``indices`` is provided.
+
+    Examples
+    --------
+    >>> linestrings([[[0, 1], [4, 5]], [[2, 3], [5, 6]]]).tolist()
+    [<LINESTRING (0 1, 4 5)>, <LINESTRING (2 3, 5 6)>]
+    >>> linestrings([[0, 1], [4, 5], [2, 3], [5, 6], [7, 8]], indices=[0, 0, 1, 1, 1]).tolist()
+    [<LINESTRING (0 1, 4 5)>, <LINESTRING (2 3, 5 6, 7 8)>]
+
+    Notes
+    -----
+    - Usage of the ``y`` and ``z`` arguments will prevents lazy evaluation in ``dask``.
+      Instead provide the coordinates as a ``(..., 2)`` or ``(..., 3)`` array using only ``coords``.
+    """
+    coords = _xyz_to_coords(coords, y, z)
+    if indices is None:
+        return lib.linestrings(coords, out=out, **kwargs)
+    else:
+        return simple_geometries_1d(coords, indices, GeometryType.LINESTRING, out=out)
+
+
+@multithreading_enabled
+def linearrings(coords, y=None, z=None, indices=None, out=None, **kwargs):
+    """Create an array of linearrings.
+
+    If the provided coords do not constitute a closed linestring, or if there
+    are only 3 provided coords, the first
+    coordinate is duplicated at the end to close the ring. This function will
+    raise an exception if a linearring contains less than three points or if
+    the terminal coordinates contain NaN (not-a-number).
+
+    Parameters
+    ----------
+    coords : array_like
+        An array of lists of coordinate tuples (2- or 3-dimensional) or, if ``y``
+        is provided, an array of lists of x coordinates
+    y : array_like, optional
+    z : array_like, optional
+    indices : array_like, optional
+        Indices into the target array where input coordinates belong. If
+        provided, the coords should be 2D with shape (N, 2) or (N, 3) and
+        indices should be an array of shape (N,) with integers in increasing
+        order. Missing indices result in a ValueError unless ``out`` is
+        provided, in which case the original value in ``out`` is kept.
+    out : ndarray, optional
+        An array (with dtype object) to output the geometries into.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+        Ignored if ``indices`` is provided.
+
+    See also
+    --------
+    linestrings
+
+    Examples
+    --------
+    >>> linearrings([[0, 0], [0, 1], [1, 1], [0, 0]])
+    <LINEARRING (0 0, 0 1, 1 1, 0 0)>
+    >>> linearrings([[0, 0], [0, 1], [1, 1]])
+    <LINEARRING (0 0, 0 1, 1 1, 0 0)>
+
+    Notes
+    -----
+    - Usage of the ``y`` and ``z`` arguments will prevents lazy evaluation in ``dask``.
+      Instead provide the coordinates as a ``(..., 2)`` or ``(..., 3)`` array using only ``coords``.
+    """
+    coords = _xyz_to_coords(coords, y, z)
+    if indices is None:
+        return lib.linearrings(coords, out=out, **kwargs)
+    else:
+        return simple_geometries_1d(coords, indices, GeometryType.LINEARRING, out=out)
+
+
+@multithreading_enabled
+def polygons(geometries, holes=None, indices=None, out=None, **kwargs):
+    """Create an array of polygons.
+
+    Parameters
+    ----------
+    geometries : array_like
+        An array of linearrings or coordinates (see linearrings).
+        Unless ``indices`` are given (see description below), this
+        include the outer shells only. The ``holes`` argument should be used
+        to create polygons with holes.
+    holes : array_like, optional
+        An array of lists of linearrings that constitute holes for each shell.
+        Not to be used in combination with ``indices``.
+    indices : array_like, optional
+        Indices into the target array where input geometries belong. If
+        provided, the holes are expected to be present inside ``geometries``;
+        the first geometry for each index is the outer shell
+        and all subsequent geometries in that index are the holes.
+        Both geometries and indices should be 1D and have matching sizes.
+        Indices should be in increasing order. Missing indices result in a ValueError
+        unless ``out`` is  provided, in which case the original value in ``out`` is kept.
+    out : ndarray, optional
+        An array (with dtype object) to output the geometries into.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+        Ignored if ``indices`` is provided.
+
+    Examples
+    --------
+    Polygons are constructed from rings:
+
+    >>> ring_1 = linearrings([[0, 0], [0, 10], [10, 10], [10, 0]])
+    >>> ring_2 = linearrings([[2, 6], [2, 7], [3, 7], [3, 6]])
+    >>> polygons([ring_1, ring_2])[0]
+    <POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))>
+    >>> polygons([ring_1, ring_2])[1]
+    <POLYGON ((2 6, 2 7, 3 7, 3 6, 2 6))>
+
+    Or from coordinates directly:
+
+    >>> polygons([[0, 0], [0, 10], [10, 10], [10, 0]])
+    <POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))>
+
+    Adding holes can be done using the ``holes`` keyword argument:
+
+    >>> polygons(ring_1, holes=[ring_2])
+    <POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0), (2 6, 2 7, 3 7, 3 6, 2 6))>
+
+    Or using the ``indices`` argument:
+
+    >>> polygons([ring_1, ring_2], indices=[0, 1])[0]
+    <POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))>
+    >>> polygons([ring_1, ring_2], indices=[0, 1])[1]
+    <POLYGON ((2 6, 2 7, 3 7, 3 6, 2 6))>
+    >>> polygons([ring_1, ring_2], indices=[0, 0])[0]
+    <POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0), (2 6, 2 7, 3 7, 3 6, 2 6))>
+
+    Missing input values (``None``) are ignored and may result in an
+    empty polygon:
+
+    >>> polygons(None)
+    <POLYGON EMPTY>
+    >>> polygons(ring_1, holes=[None])
+    <POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))>
+    >>> polygons([ring_1, None], indices=[0, 0])[0]
+    <POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))>
+    """
+    geometries = np.asarray(geometries)
+    if not isinstance(geometries, Geometry) and np.issubdtype(
+        geometries.dtype, np.number
+    ):
+        geometries = linearrings(geometries)
+
+    if indices is not None:
+        if holes is not None:
+            raise TypeError("Cannot specify separate holes array when using indices.")
+        return collections_1d(geometries, indices, GeometryType.POLYGON, out=out)
+
+    if holes is None:
+        # no holes provided: initialize an empty holes array matching shells
+        shape = geometries.shape + (0,) if isinstance(geometries, np.ndarray) else (0,)
+        holes = np.empty(shape, dtype=object)
+    else:
+        holes = np.asarray(holes)
+        # convert holes coordinates into linearrings
+        if np.issubdtype(holes.dtype, np.number):
+            holes = linearrings(holes)
+
+    return lib.polygons(geometries, holes, out=out, **kwargs)
+
+
+@multithreading_enabled
+def box(xmin, ymin, xmax, ymax, ccw=True, **kwargs):
+    """Create box polygons.
+
+    Parameters
+    ----------
+    xmin : array_like
+    ymin : array_like
+    xmax : array_like
+    ymax : array_like
+    ccw : bool, default True
+        If True, box will be created in counterclockwise direction starting
+        from bottom right coordinate (xmax, ymin).
+        If False, box will be created in clockwise direction starting from
+        bottom left coordinate (xmin, ymin).
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> box(0, 0, 1, 1)
+    <POLYGON ((1 0, 1 1, 0 1, 0 0, 1 0))>
+    >>> box(0, 0, 1, 1, ccw=False)
+    <POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))>
+
+    """
+    return lib.box(xmin, ymin, xmax, ymax, ccw, **kwargs)
+
+
+@multithreading_enabled
+def multipoints(geometries, indices=None, out=None, **kwargs):
+    """Create multipoints from arrays of points
+
+    Parameters
+    ----------
+    geometries : array_like
+        An array of points or coordinates (see points).
+    indices : array_like, optional
+        Indices into the target array where input geometries belong. If
+        provided, both geometries and indices should be 1D and have matching
+        sizes. Indices should be in increasing order. Missing indices result
+        in a ValueError unless ``out`` is  provided, in which case the original
+        value in ``out`` is kept.
+    out : ndarray, optional
+        An array (with dtype object) to output the geometries into.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+        Ignored if ``indices`` is provided.
+
+    Examples
+    --------
+    Multipoints are constructed from points:
+
+    >>> point_1 = points([1, 1])
+    >>> point_2 = points([2, 2])
+    >>> multipoints([point_1, point_2])
+    <MULTIPOINT (1 1, 2 2)>
+    >>> multipoints([[point_1, point_2], [point_2, None]]).tolist()
+    [<MULTIPOINT (1 1, 2 2)>, <MULTIPOINT (2 2)>]
+
+    Or from coordinates directly:
+
+    >>> multipoints([[0, 0], [2, 2], [3, 3]])
+    <MULTIPOINT (0 0, 2 2, 3 3)>
+
+    Multiple multipoints of different sizes can be constructed efficiently using the
+    ``indices`` keyword argument:
+
+    >>> multipoints([point_1, point_2, point_2], indices=[0, 0, 1]).tolist()
+    [<MULTIPOINT (1 1, 2 2)>, <MULTIPOINT (2 2)>]
+
+    Missing input values (``None``) are ignored and may result in an
+    empty multipoint:
+
+    >>> multipoints([None])
+    <MULTIPOINT EMPTY>
+    >>> multipoints([point_1, None], indices=[0, 0]).tolist()
+    [<MULTIPOINT (1 1)>]
+    >>> multipoints([point_1, None], indices=[0, 1]).tolist()
+    [<MULTIPOINT (1 1)>, <MULTIPOINT EMPTY>]
+    """
+    typ = GeometryType.MULTIPOINT
+    geometries = np.asarray(geometries)
+    if not isinstance(geometries, Geometry) and np.issubdtype(
+        geometries.dtype, np.number
+    ):
+        geometries = points(geometries)
+    if indices is None:
+        return lib.create_collection(geometries, typ, out=out, **kwargs)
+    else:
+        return collections_1d(geometries, indices, typ, out=out)
+
+
+@multithreading_enabled
+def multilinestrings(geometries, indices=None, out=None, **kwargs):
+    """Create multilinestrings from arrays of linestrings
+
+    Parameters
+    ----------
+    geometries : array_like
+        An array of linestrings or coordinates (see linestrings).
+    indices : array_like, optional
+        Indices into the target array where input geometries belong. If
+        provided, both geometries and indices should be 1D and have matching
+        sizes. Indices should be in increasing order. Missing indices result
+        in a ValueError unless ``out`` is  provided, in which case the original
+        value in ``out`` is kept.
+    out : ndarray, optional
+        An array (with dtype object) to output the geometries into.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+        Ignored if ``indices`` is provided.
+
+    See also
+    --------
+    multipoints
+    """
+    typ = GeometryType.MULTILINESTRING
+    geometries = np.asarray(geometries)
+    if not isinstance(geometries, Geometry) and np.issubdtype(
+        geometries.dtype, np.number
+    ):
+        geometries = linestrings(geometries)
+
+    if indices is None:
+        return lib.create_collection(geometries, typ, out=out, **kwargs)
+    else:
+        return collections_1d(geometries, indices, typ, out=out)
+
+
+@multithreading_enabled
+def multipolygons(geometries, indices=None, out=None, **kwargs):
+    """Create multipolygons from arrays of polygons
+
+    Parameters
+    ----------
+    geometries : array_like
+        An array of polygons or coordinates (see polygons).
+    indices : array_like, optional
+        Indices into the target array where input geometries belong. If
+        provided, both geometries and indices should be 1D and have matching
+        sizes. Indices should be in increasing order. Missing indices result
+        in a ValueError unless ``out`` is  provided, in which case the original
+        value in ``out`` is kept.
+    out : ndarray, optional
+        An array (with dtype object) to output the geometries into.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+        Ignored if ``indices`` is provided.
+
+    See also
+    --------
+    multipoints
+    """
+    typ = GeometryType.MULTIPOLYGON
+    geometries = np.asarray(geometries)
+    if not isinstance(geometries, Geometry) and np.issubdtype(
+        geometries.dtype, np.number
+    ):
+        geometries = polygons(geometries)
+    if indices is None:
+        return lib.create_collection(geometries, typ, out=out, **kwargs)
+    else:
+        return collections_1d(geometries, indices, typ, out=out)
+
+
+@multithreading_enabled
+def geometrycollections(geometries, indices=None, out=None, **kwargs):
+    """Create geometrycollections from arrays of geometries
+
+    Parameters
+    ----------
+    geometries : array_like
+        An array of geometries
+    indices : array_like, optional
+        Indices into the target array where input geometries belong. If
+        provided, both geometries and indices should be 1D and have matching
+        sizes. Indices should be in increasing order. Missing indices result
+        in a ValueError unless ``out`` is  provided, in which case the original
+        value in ``out`` is kept.
+    out : ndarray, optional
+        An array (with dtype object) to output the geometries into.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+        Ignored if ``indices`` is provided.
+
+    See also
+    --------
+    multipoints
+    """
+    typ = GeometryType.GEOMETRYCOLLECTION
+    if indices is None:
+        return lib.create_collection(geometries, typ, out=out, **kwargs)
+    else:
+        return collections_1d(geometries, indices, typ, out=out)
+
+
+def prepare(geometry, **kwargs):
+    """Prepare a geometry, improving performance of other operations.
+
+    A prepared geometry is a normal geometry with added information such as an
+    index on the line segments. This improves the performance of the following operations:
+    contains, contains_properly, covered_by, covers, crosses, disjoint, intersects,
+    overlaps, touches, and within.
+
+    Note that if a prepared geometry is modified, the newly created Geometry object is
+    not prepared. In that case, ``prepare`` should be called again.
+
+    This function does not recompute previously prepared geometries;
+    it is efficient to call this function on an array that partially contains prepared geometries.
+
+    This function does not return any values; geometries are modified in place.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+        Geometries are changed in place
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    is_prepared : Identify whether a geometry is prepared already.
+    destroy_prepared : Destroy the prepared part of a geometry.
+
+    Examples
+    --------
+    >>> from shapely import Point, buffer, prepare, contains_properly
+    >>> poly = buffer(Point(1.0, 1.0), 1)
+    >>> prepare(poly)
+    >>> contains_properly(poly, [Point(0.0, 0.0), Point(0.5, 0.5)]).tolist()
+    [False, True]
+    """
+    lib.prepare(geometry, **kwargs)
+
+
+def destroy_prepared(geometry, **kwargs):
+    """Destroy the prepared part of a geometry, freeing up memory.
+
+    Note that the prepared geometry will always be cleaned up if the geometry itself
+    is dereferenced. This function needs only be called in very specific circumstances,
+    such as freeing up memory without losing the geometries, or benchmarking.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+        Geometries are changed inplace
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    prepare
+    """
+    lib.destroy_prepared(geometry, **kwargs)
+
+
+def empty(shape, geom_type=None, order="C"):
+    """Create a geometry array prefilled with None or with empty geometries.
+
+    Parameters
+    ----------
+    shape : int or tuple of int
+        Shape of the empty array, e.g., ``(2, 3)`` or ``2``.
+    geom_type : shapely.GeometryType, optional
+        The desired geometry type in case the array should be prefilled
+        with empty geometries. Default ``None``.
+    order : {'C', 'F'}, optional, default: 'C'
+        Whether to store multi-dimensional data in row-major
+        (C-style) or column-major (Fortran-style) order in
+        memory.
+
+    Examples
+    --------
+    >>> empty((2, 3)).tolist()
+    [[None, None, None], [None, None, None]]
+    >>> empty(2, geom_type=GeometryType.POINT).tolist()
+    [<POINT EMPTY>, <POINT EMPTY>]
+    """
+    if geom_type is None:
+        return np.empty(shape, dtype=object, order=order)
+
+    geom_type = GeometryType(geom_type)  # cast int to GeometryType
+    if geom_type is GeometryType.MISSING:
+        return np.empty(shape, dtype=object, order=order)
+
+    fill_value = from_wkt(geom_type.name + " EMPTY")
+    return np.full(shape, fill_value, dtype=object, order=order)

billinglayer/python/shapely/decorators.py

@@ -0,0 +1,82 @@
+import os
+from functools import wraps
+
+import numpy as np
+
+from . import lib
+from .errors import UnsupportedGEOSVersionError
+
+
+class requires_geos:
+    def __init__(self, version):
+        if version.count(".") != 2:
+            raise ValueError("Version must be <major>.<minor>.<patch> format")
+        self.version = tuple(int(x) for x in version.split("."))
+
+    def __call__(self, func):
+        is_compatible = lib.geos_version >= self.version
+        is_doc_build = os.environ.get("SPHINX_DOC_BUILD") == "1"  # set in docs/conf.py
+        if is_compatible and not is_doc_build:
+            return func  # return directly, do not change the docstring
+
+        msg = "'{}' requires at least GEOS {}.{}.{}.".format(
+            func.__name__, *self.version
+        )
+        if is_compatible:
+
+            @wraps(func)
+            def wrapped(*args, **kwargs):
+                return func(*args, **kwargs)
+
+        else:
+
+            @wraps(func)
+            def wrapped(*args, **kwargs):
+                raise UnsupportedGEOSVersionError(msg)
+
+        doc = wrapped.__doc__
+        if doc:
+            # Insert the message at the first double newline
+            position = doc.find("\n\n") + 2
+            # Figure out the indentation level
+            indent = 2
+            while True:
+                if doc[position + indent] == " ":
+                    indent += 1
+                else:
+                    break
+            wrapped.__doc__ = doc.replace(
+                "\n\n", "\n\n{}.. note:: {}\n\n".format(" " * indent, msg), 1
+            )
+
+        return wrapped
+
+
+def multithreading_enabled(func):
+    """Prepare multithreading by setting the writable flags of object type
+    ndarrays to False.
+
+    NB: multithreading also requires the GIL to be released, which is done in
+    the C extension (ufuncs.c)."""
+
+    @wraps(func)
+    def wrapped(*args, **kwargs):
+        array_args = [
+            arg for arg in args if isinstance(arg, np.ndarray) and arg.dtype == object
+        ] + [
+            arg
+            for name, arg in kwargs.items()
+            if name not in {"where", "out"}
+            and isinstance(arg, np.ndarray)
+            and arg.dtype == object
+        ]
+        old_flags = [arr.flags.writeable for arr in array_args]
+        try:
+            for arr in array_args:
+                arr.flags.writeable = False
+            return func(*args, **kwargs)
+        finally:
+            for arr, old_flag in zip(array_args, old_flags):
+                arr.flags.writeable = old_flag
+
+    return wrapped

billinglayer/python/shapely/errors.py

@@ -0,0 +1,81 @@
+"""Shapely errors."""
+import threading
+
+from shapely.lib import _setup_signal_checks, GEOSException, ShapelyError  # NOQA
+
+
+def setup_signal_checks(interval=10000):
+    """This enables Python signal checks in the ufunc inner loops.
+
+    Doing so allows termination (using CTRL+C) of operations on large arrays of vectors.
+
+    Parameters
+    ----------
+    interval : int, default 10000
+        Check for interrupts every x iterations. The higher the number, the slower
+        shapely will respond to a signal. However, at low values there will be a negative effect
+        on performance. The default of 10000 does not have any measureable effects on performance.
+
+    Notes
+    -----
+    For more information on signals consult the Python docs:
+
+    https://docs.python.org/3/library/signal.html
+    """
+    if interval <= 0:
+        raise ValueError("Signal checks interval must be greater than zero.")
+
+    _setup_signal_checks(interval, threading.main_thread().ident)
+
+
+class UnsupportedGEOSVersionError(ShapelyError):
+    """Raised when the GEOS library version does not support a certain operation."""
+
+
+class DimensionError(ShapelyError):
+    """An error in the number of coordinate dimensions."""
+
+
+class TopologicalError(ShapelyError):
+    """A geometry is invalid or topologically incorrect."""
+
+
+class ShapelyDeprecationWarning(FutureWarning):
+    """
+    Warning for features that will be removed or behaviour that will be
+    changed in a future release.
+    """
+
+
+class EmptyPartError(ShapelyError):
+    """An error signifying an empty part was encountered when creating a multi-part."""
+
+
+class GeometryTypeError(ShapelyError):
+    """
+    An error raised when the type of the geometry in question is
+    unrecognized or inappropriate.
+    """
+
+
+def __getattr__(name):
+    import warnings
+
+    # Alias Shapely 1.8 error classes to ShapelyError with deprecation warning
+    if name in [
+        "ReadingError",
+        "WKBReadingError",
+        "WKTReadingError",
+        "PredicateError",
+        "InvalidGeometryError",
+    ]:
+        warnings.warn(
+            f"{name} is deprecated and will be removed in a future version. "
+            "Use ShapelyError instead (functions previously raising {name} "
+            "will now raise a ShapelyError instead).",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        return ShapelyError
+
+    raise AttributeError(f"module 'shapely.errors' has no attribute '{name}'")

billinglayer/python/shapely/geometry/__init__.py

@@ -0,0 +1,28 @@
+"""Geometry classes and factories
+"""
+
+from .base import CAP_STYLE, JOIN_STYLE
+from .collection import GeometryCollection
+from .geo import box, mapping, shape
+from .linestring import LineString
+from .multilinestring import MultiLineString
+from .multipoint import MultiPoint
+from .multipolygon import MultiPolygon
+from .point import Point
+from .polygon import LinearRing, Polygon
+
+__all__ = [
+    "box",
+    "shape",
+    "mapping",
+    "Point",
+    "LineString",
+    "Polygon",
+    "MultiPoint",
+    "MultiLineString",
+    "MultiPolygon",
+    "GeometryCollection",
+    "LinearRing",
+    "CAP_STYLE",
+    "JOIN_STYLE",
+]

billinglayer/python/shapely/geometry/base.py

@@ -0,0 +1,994 @@
+"""Base geometry class and utilities
+
+Note: a third, z, coordinate value may be used when constructing
+geometry objects, but has no effect on geometric analysis. All
+operations are performed in the x-y plane. Thus, geometries with
+different z values may intersect or be equal.
+"""
+import re
+from warnings import warn
+
+import numpy as np
+
+import shapely
+from shapely._geometry_helpers import _geom_factory
+from shapely.constructive import BufferCapStyle, BufferJoinStyle
+from shapely.coords import CoordinateSequence
+from shapely.errors import GeometryTypeError, GEOSException, ShapelyDeprecationWarning
+
+GEOMETRY_TYPES = [
+    "Point",
+    "LineString",
+    "LinearRing",
+    "Polygon",
+    "MultiPoint",
+    "MultiLineString",
+    "MultiPolygon",
+    "GeometryCollection",
+]
+
+
+def geom_factory(g, parent=None):
+    """
+    Creates a Shapely geometry instance from a pointer to a GEOS geometry.
+
+    .. warning::
+        The GEOS library used to create the the GEOS geometry pointer
+        and the GEOS library used by Shapely must be exactly the same, or
+        unexpected results or segfaults may occur.
+
+    .. deprecated:: 2.0
+        Deprecated in Shapely 2.0, and will be removed in a future version.
+    """
+    warn(
+        "The 'geom_factory' function is deprecated in Shapely 2.0, and will be "
+        "removed in a future version",
+        DeprecationWarning,
+        stacklevel=2,
+    )
+    return _geom_factory(g)
+
+
+def dump_coords(geom):
+    """Dump coordinates of a geometry in the same order as data packing"""
+    if not isinstance(geom, BaseGeometry):
+        raise ValueError(
+            "Must be instance of a geometry class; found " + geom.__class__.__name__
+        )
+    elif geom.geom_type in ("Point", "LineString", "LinearRing"):
+        return geom.coords[:]
+    elif geom.geom_type == "Polygon":
+        return geom.exterior.coords[:] + [i.coords[:] for i in geom.interiors]
+    elif geom.geom_type.startswith("Multi") or geom.geom_type == "GeometryCollection":
+        # Recursive call
+        return [dump_coords(part) for part in geom.geoms]
+    else:
+        raise GeometryTypeError("Unhandled geometry type: " + repr(geom.geom_type))
+
+
+def _maybe_unpack(result):
+    if result.ndim == 0:
+        # convert numpy 0-d array / scalar to python scalar
+        return result.item()
+    else:
+        # >=1 dim array
+        return result
+
+
+class CAP_STYLE:
+    round = BufferCapStyle.round
+    flat = BufferCapStyle.flat
+    square = BufferCapStyle.square
+
+
+class JOIN_STYLE:
+    round = BufferJoinStyle.round
+    mitre = BufferJoinStyle.mitre
+    bevel = BufferJoinStyle.bevel
+
+
+class BaseGeometry(shapely.Geometry):
+    """
+    Provides GEOS spatial predicates and topological operations.
+
+    """
+
+    __slots__ = []
+
+    def __new__(self):
+        warn(
+            "Directly calling the base class 'BaseGeometry()' is deprecated, and "
+            "will raise an error in the future. To create an empty geometry, "
+            "use one of the subclasses instead, for example 'GeometryCollection()'.",
+            ShapelyDeprecationWarning,
+            stacklevel=2,
+        )
+        return shapely.from_wkt("GEOMETRYCOLLECTION EMPTY")
+
+    @property
+    def _ndim(self):
+        return shapely.get_coordinate_dimension(self)
+
+    def __bool__(self):
+        return self.is_empty is False
+
+    def __nonzero__(self):
+        return self.__bool__()
+
+    def __format__(self, format_spec):
+        """Format a geometry using a format specification."""
+        # bypass regexp for simple cases
+        if format_spec == "":
+            return shapely.to_wkt(self, rounding_precision=-1)
+        elif format_spec == "x":
+            return shapely.to_wkb(self, hex=True).lower()
+        elif format_spec == "X":
+            return shapely.to_wkb(self, hex=True)
+
+        # fmt: off
+        format_spec_regexp = (
+            "(?:0?\\.(?P<prec>[0-9]+))?"
+            "(?P<fmt_code>[fFgGxX]?)"
+        )
+        # fmt: on
+        match = re.fullmatch(format_spec_regexp, format_spec)
+        if match is None:
+            raise ValueError(f"invalid format specifier: {format_spec}")
+
+        prec, fmt_code = match.groups()
+
+        if prec:
+            prec = int(prec)
+        else:
+            # GEOS has a default rounding_precision -1
+            prec = -1
+
+        if not fmt_code:
+            fmt_code = "g"
+
+        if fmt_code in ("g", "G"):
+            res = shapely.to_wkt(self, rounding_precision=prec, trim=True)
+        elif fmt_code in ("f", "F"):
+            res = shapely.to_wkt(self, rounding_precision=prec, trim=False)
+        elif fmt_code in ("x", "X"):
+            raise ValueError("hex representation does not specify precision")
+        else:
+            raise NotImplementedError(f"unhandled fmt_code: {fmt_code}")
+
+        if fmt_code.isupper():
+            return res.upper()
+        else:
+            return res
+
+    def __repr__(self):
+        try:
+            wkt = super().__str__()
+        except (GEOSException, ValueError):
+            # we never want a repr() to fail; that can be very confusing
+            return "<shapely.{} Exception in WKT writer>".format(
+                self.__class__.__name__
+            )
+
+        # the total length is limited to 80 characters including brackets
+        max_length = 78
+        if len(wkt) > max_length:
+            return f"<{wkt[: max_length - 3]}...>"
+
+        return f"<{wkt}>"
+
+    def __str__(self):
+        return self.wkt
+
+    def __reduce__(self):
+        return (shapely.from_wkb, (shapely.to_wkb(self, include_srid=True),))
+
+    # Operators
+    # ---------
+
+    def __and__(self, other):
+        return self.intersection(other)
+
+    def __or__(self, other):
+        return self.union(other)
+
+    def __sub__(self, other):
+        return self.difference(other)
+
+    def __xor__(self, other):
+        return self.symmetric_difference(other)
+
+    # Coordinate access
+    # -----------------
+
+    @property
+    def coords(self):
+        """Access to geometry's coordinates (CoordinateSequence)"""
+        coords_array = shapely.get_coordinates(self, include_z=self.has_z)
+        return CoordinateSequence(coords_array)
+
+    @property
+    def xy(self):
+        """Separate arrays of X and Y coordinate values"""
+        raise NotImplementedError
+
+    # Python feature protocol
+
+    @property
+    def __geo_interface__(self):
+        """Dictionary representation of the geometry"""
+        raise NotImplementedError
+
+    # Type of geometry and its representations
+    # ----------------------------------------
+
+    def geometryType(self):
+        warn(
+            "The 'GeometryType()' method is deprecated, and will be removed in "
+            "the future. You can use the 'geom_type' attribute instead.",
+            ShapelyDeprecationWarning,
+            stacklevel=2,
+        )
+        return self.geom_type
+
+    @property
+    def type(self):
+        warn(
+            "The 'type' attribute is deprecated, and will be removed in "
+            "the future. You can use the 'geom_type' attribute instead.",
+            ShapelyDeprecationWarning,
+            stacklevel=2,
+        )
+        return self.geom_type
+
+    @property
+    def wkt(self):
+        """WKT representation of the geometry"""
+        # TODO(shapely-2.0) keep default of not trimming?
+        return shapely.to_wkt(self, rounding_precision=-1)
+
+    @property
+    def wkb(self):
+        """WKB representation of the geometry"""
+        return shapely.to_wkb(self)
+
+    @property
+    def wkb_hex(self):
+        """WKB hex representation of the geometry"""
+        return shapely.to_wkb(self, hex=True)
+
+    def svg(self, scale_factor=1.0, **kwargs):
+        """Raises NotImplementedError"""
+        raise NotImplementedError
+
+    def _repr_svg_(self):
+        """SVG representation for iPython notebook"""
+        svg_top = (
+            '<svg xmlns="http://www.w3.org/2000/svg" '
+            'xmlns:xlink="http://www.w3.org/1999/xlink" '
+        )
+        if self.is_empty:
+            return svg_top + "/>"
+        else:
+            # Establish SVG canvas that will fit all the data + small space
+            xmin, ymin, xmax, ymax = self.bounds
+            if xmin == xmax and ymin == ymax:
+                # This is a point; buffer using an arbitrary size
+                xmin, ymin, xmax, ymax = self.buffer(1).bounds
+            else:
+                # Expand bounds by a fraction of the data ranges
+                expand = 0.04  # or 4%, same as R plots
+                widest_part = max([xmax - xmin, ymax - ymin])
+                expand_amount = widest_part * expand
+                xmin -= expand_amount
+                ymin -= expand_amount
+                xmax += expand_amount
+                ymax += expand_amount
+            dx = xmax - xmin
+            dy = ymax - ymin
+            width = min([max([100.0, dx]), 300])
+            height = min([max([100.0, dy]), 300])
+            try:
+                scale_factor = max([dx, dy]) / max([width, height])
+            except ZeroDivisionError:
+                scale_factor = 1.0
+            view_box = f"{xmin} {ymin} {dx} {dy}"
+            transform = f"matrix(1,0,0,-1,0,{ymax + ymin})"
+            return svg_top + (
+                'width="{1}" height="{2}" viewBox="{0}" '
+                'preserveAspectRatio="xMinYMin meet">'
+                '<g transform="{3}">{4}</g></svg>'
+            ).format(view_box, width, height, transform, self.svg(scale_factor))
+
+    @property
+    def geom_type(self):
+        """Name of the geometry's type, such as 'Point'"""
+        return GEOMETRY_TYPES[shapely.get_type_id(self)]
+
+    # Real-valued properties and methods
+    # ----------------------------------
+
+    @property
+    def area(self):
+        """Unitless area of the geometry (float)"""
+        return float(shapely.area(self))
+
+    def distance(self, other):
+        """Unitless distance to other geometry (float)"""
+        return _maybe_unpack(shapely.distance(self, other))
+
+    def hausdorff_distance(self, other):
+        """Unitless hausdorff distance to other geometry (float)"""
+        return _maybe_unpack(shapely.hausdorff_distance(self, other))
+
+    @property
+    def length(self):
+        """Unitless length of the geometry (float)"""
+        return float(shapely.length(self))
+
+    @property
+    def minimum_clearance(self):
+        """Unitless distance by which a node could be moved to produce an invalid geometry (float)"""
+        return float(shapely.minimum_clearance(self))
+
+    # Topological properties
+    # ----------------------
+
+    @property
+    def boundary(self):
+        """Returns a lower dimension geometry that bounds the object
+
+        The boundary of a polygon is a line, the boundary of a line is a
+        collection of points. The boundary of a point is an empty (null)
+        collection.
+        """
+        return shapely.boundary(self)
+
+    @property
+    def bounds(self):
+        """Returns minimum bounding region (minx, miny, maxx, maxy)"""
+        return tuple(shapely.bounds(self).tolist())
+
+    @property
+    def centroid(self):
+        """Returns the geometric center of the object"""
+        return shapely.centroid(self)
+
+    def point_on_surface(self):
+        """Returns a point guaranteed to be within the object, cheaply.
+
+        Alias of `representative_point`.
+        """
+        return shapely.point_on_surface(self)
+
+    def representative_point(self):
+        """Returns a point guaranteed to be within the object, cheaply.
+
+        Alias of `point_on_surface`.
+        """
+        return shapely.point_on_surface(self)
+
+    @property
+    def convex_hull(self):
+        """Imagine an elastic band stretched around the geometry: that's a
+        convex hull, more or less
+
+        The convex hull of a three member multipoint, for example, is a
+        triangular polygon.
+        """
+        return shapely.convex_hull(self)
+
+    @property
+    def envelope(self):
+        """A figure that envelopes the geometry"""
+        return shapely.envelope(self)
+
+    @property
+    def oriented_envelope(self):
+        """
+        Returns the oriented envelope (minimum rotated rectangle) that
+        encloses the geometry.
+
+        Unlike envelope this rectangle is not constrained to be parallel to the
+        coordinate axes. If the convex hull of the object is a degenerate (line
+        or point) this degenerate is returned.
+
+        Alias of `minimum_rotated_rectangle`.
+        """
+        return shapely.oriented_envelope(self)
+
+    @property
+    def minimum_rotated_rectangle(self):
+        """
+        Returns the oriented envelope (minimum rotated rectangle) that
+        encloses the geometry.
+
+        Unlike `envelope` this rectangle is not constrained to be parallel to the
+        coordinate axes. If the convex hull of the object is a degenerate (line
+        or point) this degenerate is returned.
+
+        Alias of `oriented_envelope`.
+        """
+        return shapely.oriented_envelope(self)
+
+    def buffer(
+        self,
+        distance,
+        quad_segs=16,
+        cap_style="round",
+        join_style="round",
+        mitre_limit=5.0,
+        single_sided=False,
+        **kwargs,
+    ):
+        """Get a geometry that represents all points within a distance
+        of this geometry.
+
+        A positive distance produces a dilation, a negative distance an
+        erosion. A very small or zero distance may sometimes be used to
+        "tidy" a polygon.
+
+        Parameters
+        ----------
+        distance : float
+            The distance to buffer around the object.
+        resolution : int, optional
+            The resolution of the buffer around each vertex of the
+            object.
+        quad_segs : int, optional
+            Sets the number of line segments used to approximate an
+            angle fillet.
+        cap_style : shapely.BufferCapStyle or {'round', 'square', 'flat'}, default 'round'
+            Specifies the shape of buffered line endings. BufferCapStyle.round ('round')
+            results in circular line endings (see ``quad_segs``). Both BufferCapStyle.square
+            ('square') and BufferCapStyle.flat ('flat') result in rectangular line endings,
+            only BufferCapStyle.flat ('flat') will end at the original vertex,
+            while BufferCapStyle.square ('square') involves adding the buffer width.
+        join_style : shapely.BufferJoinStyle or {'round', 'mitre', 'bevel'}, default 'round'
+            Specifies the shape of buffered line midpoints. BufferJoinStyle.ROUND ('round')
+            results in rounded shapes. BufferJoinStyle.bevel ('bevel') results in a beveled
+            edge that touches the original vertex. BufferJoinStyle.mitre ('mitre') results
+            in a single vertex that is beveled depending on the ``mitre_limit`` parameter.
+        mitre_limit : float, optional
+            The mitre limit ratio is used for very sharp corners. The
+            mitre ratio is the ratio of the distance from the corner to
+            the end of the mitred offset corner. When two line segments
+            meet at a sharp angle, a miter join will extend the original
+            geometry. To prevent unreasonable geometry, the mitre limit
+            allows controlling the maximum length of the join corner.
+            Corners with a ratio which exceed the limit will be beveled.
+        single_side : bool, optional
+            The side used is determined by the sign of the buffer
+            distance:
+
+                a positive distance indicates the left-hand side
+                a negative distance indicates the right-hand side
+
+            The single-sided buffer of point geometries is the same as
+            the regular buffer.  The End Cap Style for single-sided
+            buffers is always ignored, and forced to the equivalent of
+            CAP_FLAT.
+        quadsegs : int, optional
+            Deprecated alias for `quad_segs`.
+
+        Returns
+        -------
+        Geometry
+
+        Notes
+        -----
+        The return value is a strictly two-dimensional geometry. All
+        Z coordinates of the original geometry will be ignored.
+
+        Examples
+        --------
+        >>> from shapely.wkt import loads
+        >>> g = loads('POINT (0.0 0.0)')
+
+        16-gon approx of a unit radius circle:
+
+        >>> g.buffer(1.0).area  # doctest: +ELLIPSIS
+        3.1365484905459...
+
+        128-gon approximation:
+
+        >>> g.buffer(1.0, 128).area  # doctest: +ELLIPSIS
+        3.141513801144...
+
+        triangle approximation:
+
+        >>> g.buffer(1.0, 3).area
+        3.0
+        >>> list(g.buffer(1.0, cap_style=BufferCapStyle.square).exterior.coords)
+        [(1.0, 1.0), (1.0, -1.0), (-1.0, -1.0), (-1.0, 1.0), (1.0, 1.0)]
+        >>> g.buffer(1.0, cap_style=BufferCapStyle.square).area
+        4.0
+
+        """
+        quadsegs = kwargs.pop("quadsegs", None)
+        if quadsegs is not None:
+            warn(
+                "The `quadsegs` argument is deprecated. Use `quad_segs` instead.",
+                FutureWarning,
+            )
+            quad_segs = quadsegs
+
+        # TODO deprecate `resolution` keyword for shapely 2.1
+        resolution = kwargs.pop("resolution", None)
+        if resolution is not None:
+            quad_segs = resolution
+        if kwargs:
+            kwarg = list(kwargs.keys())[0]  # noqa
+            raise TypeError(f"buffer() got an unexpected keyword argument '{kwarg}'")
+
+        if mitre_limit == 0.0:
+            raise ValueError("Cannot compute offset from zero-length line segment")
+        elif not np.isfinite(distance).all():
+            raise ValueError("buffer distance must be finite")
+
+        return shapely.buffer(
+            self,
+            distance,
+            quad_segs=quad_segs,
+            cap_style=cap_style,
+            join_style=join_style,
+            mitre_limit=mitre_limit,
+            single_sided=single_sided,
+        )
+
+    def simplify(self, tolerance, preserve_topology=True):
+        """Returns a simplified geometry produced by the Douglas-Peucker
+        algorithm
+
+        Coordinates of the simplified geometry will be no more than the
+        tolerance distance from the original. Unless the topology preserving
+        option is used, the algorithm may produce self-intersecting or
+        otherwise invalid geometries.
+        """
+        return shapely.simplify(self, tolerance, preserve_topology=preserve_topology)
+
+    def normalize(self):
+        """Converts geometry to normal form (or canonical form).
+
+        This method orders the coordinates, rings of a polygon and parts of
+        multi geometries consistently. Typically useful for testing purposes
+        (for example in combination with `equals_exact`).
+
+        Examples
+        --------
+        >>> from shapely import MultiLineString
+        >>> line = MultiLineString([[(0, 0), (1, 1)], [(3, 3), (2, 2)]])
+        >>> line.normalize()
+        <MULTILINESTRING ((2 2, 3 3), (0 0, 1 1))>
+        """
+        return shapely.normalize(self)
+
+    # Overlay operations
+    # ---------------------------
+
+    def difference(self, other, grid_size=None):
+        """
+        Returns the difference of the geometries.
+
+        Refer to `shapely.difference` for full documentation.
+        """
+        return shapely.difference(self, other, grid_size=grid_size)
+
+    def intersection(self, other, grid_size=None):
+        """
+        Returns the intersection of the geometries.
+
+        Refer to `shapely.intersection` for full documentation.
+        """
+        return shapely.intersection(self, other, grid_size=grid_size)
+
+    def symmetric_difference(self, other, grid_size=None):
+        """
+        Returns the symmetric difference of the geometries.
+
+        Refer to `shapely.symmetric_difference` for full documentation.
+        """
+        return shapely.symmetric_difference(self, other, grid_size=grid_size)
+
+    def union(self, other, grid_size=None):
+        """
+        Returns the union of the geometries.
+
+        Refer to `shapely.union` for full documentation.
+        """
+        return shapely.union(self, other, grid_size=grid_size)
+
+    # Unary predicates
+    # ----------------
+
+    @property
+    def has_z(self):
+        """True if the geometry's coordinate sequence(s) have z values (are
+        3-dimensional)"""
+        return bool(shapely.has_z(self))
+
+    @property
+    def is_empty(self):
+        """True if the set of points in this geometry is empty, else False"""
+        return bool(shapely.is_empty(self))
+
+    @property
+    def is_ring(self):
+        """True if the geometry is a closed ring, else False"""
+        return bool(shapely.is_ring(self))
+
+    @property
+    def is_closed(self):
+        """True if the geometry is closed, else False
+
+        Applicable only to 1-D geometries."""
+        if self.geom_type == "LinearRing":
+            return True
+        return bool(shapely.is_closed(self))
+
+    @property
+    def is_simple(self):
+        """True if the geometry is simple, meaning that any self-intersections
+        are only at boundary points, else False"""
+        return bool(shapely.is_simple(self))
+
+    @property
+    def is_valid(self):
+        """True if the geometry is valid (definition depends on sub-class),
+        else False"""
+        return bool(shapely.is_valid(self))
+
+    # Binary predicates
+    # -----------------
+
+    def relate(self, other):
+        """Returns the DE-9IM intersection matrix for the two geometries
+        (string)"""
+        return shapely.relate(self, other)
+
+    def covers(self, other):
+        """Returns True if the geometry covers the other, else False"""
+        return _maybe_unpack(shapely.covers(self, other))
+
+    def covered_by(self, other):
+        """Returns True if the geometry is covered by the other, else False"""
+        return _maybe_unpack(shapely.covered_by(self, other))
+
+    def contains(self, other):
+        """Returns True if the geometry contains the other, else False"""
+        return _maybe_unpack(shapely.contains(self, other))
+
+    def contains_properly(self, other):
+        """
+        Returns True if the geometry completely contains the other, with no
+        common boundary points, else False
+
+        Refer to `shapely.contains_properly` for full documentation.
+        """
+        return _maybe_unpack(shapely.contains_properly(self, other))
+
+    def crosses(self, other):
+        """Returns True if the geometries cross, else False"""
+        return _maybe_unpack(shapely.crosses(self, other))
+
+    def disjoint(self, other):
+        """Returns True if geometries are disjoint, else False"""
+        return _maybe_unpack(shapely.disjoint(self, other))
+
+    def equals(self, other):
+        """Returns True if geometries are equal, else False.
+
+        This method considers point-set equality (or topological
+        equality), and is equivalent to (self.within(other) &
+        self.contains(other)).
+
+        Examples
+        --------
+        >>> LineString(
+        ...     [(0, 0), (2, 2)]
+        ... ).equals(
+        ...     LineString([(0, 0), (1, 1), (2, 2)])
+        ... )
+        True
+
+        Returns
+        -------
+        bool
+
+        """
+        return _maybe_unpack(shapely.equals(self, other))
+
+    def intersects(self, other):
+        """Returns True if geometries intersect, else False"""
+        return _maybe_unpack(shapely.intersects(self, other))
+
+    def overlaps(self, other):
+        """Returns True if geometries overlap, else False"""
+        return _maybe_unpack(shapely.overlaps(self, other))
+
+    def touches(self, other):
+        """Returns True if geometries touch, else False"""
+        return _maybe_unpack(shapely.touches(self, other))
+
+    def within(self, other):
+        """Returns True if geometry is within the other, else False"""
+        return _maybe_unpack(shapely.within(self, other))
+
+    def dwithin(self, other, distance):
+        """
+        Returns True if geometry is within a given distance from the other, else False.
+
+        Refer to `shapely.dwithin` for full documentation.
+        """
+        return _maybe_unpack(shapely.dwithin(self, other, distance))
+
+    def equals_exact(self, other, tolerance):
+        """True if geometries are equal to within a specified
+        tolerance.
+
+        Parameters
+        ----------
+        other : BaseGeometry
+            The other geometry object in this comparison.
+        tolerance : float
+            Absolute tolerance in the same units as coordinates.
+
+        This method considers coordinate equality, which requires
+        coordinates to be equal and in the same order for all components
+        of a geometry.
+
+        Because of this it is possible for "equals()" to be True for two
+        geometries and "equals_exact()" to be False.
+
+        Examples
+        --------
+        >>> LineString(
+        ...     [(0, 0), (2, 2)]
+        ... ).equals_exact(
+        ...     LineString([(0, 0), (1, 1), (2, 2)]),
+        ...     1e-6
+        ... )
+        False
+
+        Returns
+        -------
+        bool
+
+        """
+        return _maybe_unpack(shapely.equals_exact(self, other, tolerance))
+
+    def almost_equals(self, other, decimal=6):
+        """True if geometries are equal at all coordinates to a
+        specified decimal place.
+
+        .. deprecated:: 1.8.0
+            The 'almost_equals()' method is deprecated
+            and will be removed in Shapely 2.1 because the name is
+            confusing. The 'equals_exact()' method should be used
+            instead.
+
+        Refers to approximate coordinate equality, which requires
+        coordinates to be approximately equal and in the same order for
+        all components of a geometry.
+
+        Because of this it is possible for "equals()" to be True for two
+        geometries and "almost_equals()" to be False.
+
+        Examples
+        --------
+        >>> LineString(
+        ...     [(0, 0), (2, 2)]
+        ... ).equals_exact(
+        ...     LineString([(0, 0), (1, 1), (2, 2)]),
+        ...     1e-6
+        ... )
+        False
+
+        Returns
+        -------
+        bool
+
+        """
+        warn(
+            "The 'almost_equals()' method is deprecated and will be "
+            "removed in Shapely 2.1; use 'equals_exact()' instead",
+            ShapelyDeprecationWarning,
+            stacklevel=2,
+        )
+        return self.equals_exact(other, 0.5 * 10 ** (-decimal))
+
+    def relate_pattern(self, other, pattern):
+        """Returns True if the DE-9IM string code for the relationship between
+        the geometries satisfies the pattern, else False"""
+        return _maybe_unpack(shapely.relate_pattern(self, other, pattern))
+
+    # Linear referencing
+    # ------------------
+
+    def line_locate_point(self, other, normalized=False):
+        """Returns the distance along this geometry to a point nearest the
+        specified point
+
+        If the normalized arg is True, return the distance normalized to the
+        length of the linear geometry.
+
+        Alias of `project`.
+        """
+        return shapely.line_locate_point(self, other, normalized=normalized)
+
+    def project(self, other, normalized=False):
+        """Returns the distance along this geometry to a point nearest the
+        specified point
+
+        If the normalized arg is True, return the distance normalized to the
+        length of the linear geometry.
+
+        Alias of `line_locate_point`.
+        """
+        return shapely.line_locate_point(self, other, normalized=normalized)
+
+    def line_interpolate_point(self, distance, normalized=False):
+        """Return a point at the specified distance along a linear geometry
+
+        Negative length values are taken as measured in the reverse
+        direction from the end of the geometry. Out-of-range index
+        values are handled by clamping them to the valid range of values.
+        If the normalized arg is True, the distance will be interpreted as a
+        fraction of the geometry's length.
+
+        Alias of `interpolate`.
+        """
+        return shapely.line_interpolate_point(self, distance, normalized=normalized)
+
+    def interpolate(self, distance, normalized=False):
+        """Return a point at the specified distance along a linear geometry
+
+        Negative length values are taken as measured in the reverse
+        direction from the end of the geometry. Out-of-range index
+        values are handled by clamping them to the valid range of values.
+        If the normalized arg is True, the distance will be interpreted as a
+        fraction of the geometry's length.
+
+        Alias of `line_interpolate_point`.
+        """
+        return shapely.line_interpolate_point(self, distance, normalized=normalized)
+
+    def segmentize(self, max_segment_length):
+        """Adds vertices to line segments based on maximum segment length.
+
+        Additional vertices will be added to every line segment in an input geometry
+        so that segments are no longer than the provided maximum segment length. New
+        vertices will evenly subdivide each segment.
+
+        Only linear components of input geometries are densified; other geometries
+        are returned unmodified.
+
+        Parameters
+        ----------
+        max_segment_length : float or array_like
+            Additional vertices will be added so that all line segments are no
+            longer this value.  Must be greater than 0.
+
+        Examples
+        --------
+        >>> from shapely import LineString, Polygon
+        >>> LineString([(0, 0), (0, 10)]).segmentize(max_segment_length=5)
+        <LINESTRING (0 0, 0 5, 0 10)>
+        >>> Polygon([(0, 0), (10, 0), (10, 10), (0, 10), (0, 0)]).segmentize(max_segment_length=5)
+        <POLYGON ((0 0, 5 0, 10 0, 10 5, 10 10, 5 10, 0 10, 0 5, 0 0))>
+        """
+        return shapely.segmentize(self, max_segment_length)
+
+    def reverse(self):
+        """Returns a copy of this geometry with the order of coordinates reversed.
+
+        If the geometry is a polygon with interior rings, the interior rings are also
+        reversed.
+
+        Points are unchanged.
+
+        See also
+        --------
+        is_ccw : Checks if a geometry is clockwise.
+
+        Examples
+        --------
+        >>> from shapely import LineString, Polygon
+        >>> LineString([(0, 0), (1, 2)]).reverse()
+        <LINESTRING (1 2, 0 0)>
+        >>> Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]).reverse()
+        <POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))>
+        """
+        return shapely.reverse(self)
+
+
+class BaseMultipartGeometry(BaseGeometry):
+
+    __slots__ = []
+
+    @property
+    def coords(self):
+        raise NotImplementedError(
+            "Sub-geometries may have coordinate sequences, "
+            "but multi-part geometries do not"
+        )
+
+    @property
+    def geoms(self):
+        return GeometrySequence(self)
+
+    def __bool__(self):
+        return self.is_empty is False
+
+    def svg(self, scale_factor=1.0, color=None):
+        """Returns a group of SVG elements for the multipart geometry.
+
+        Parameters
+        ==========
+        scale_factor : float
+            Multiplication factor for the SVG stroke-width.  Default is 1.
+        color : str, optional
+            Hex string for stroke or fill color. Default is to use "#66cc99"
+            if geometry is valid, and "#ff3333" if invalid.
+        """
+        if self.is_empty:
+            return "<g />"
+        if color is None:
+            color = "#66cc99" if self.is_valid else "#ff3333"
+        return "<g>" + "".join(p.svg(scale_factor, color) for p in self.geoms) + "</g>"
+
+
+class GeometrySequence:
+    """
+    Iterative access to members of a homogeneous multipart geometry.
+    """
+
+    # Attributes
+    # ----------
+    # _parent : object
+    #     Parent (Shapely) geometry
+    _parent = None
+
+    def __init__(self, parent):
+        self._parent = parent
+
+    def _get_geom_item(self, i):
+        return shapely.get_geometry(self._parent, i)
+
+    def __iter__(self):
+        for i in range(self.__len__()):
+            yield self._get_geom_item(i)
+
+    def __len__(self):
+        return shapely.get_num_geometries(self._parent)
+
+    def __getitem__(self, key):
+        m = self.__len__()
+        if isinstance(key, (int, np.integer)):
+            if key + m < 0 or key >= m:
+                raise IndexError("index out of range")
+            if key < 0:
+                i = m + key
+            else:
+                i = key
+            return self._get_geom_item(i)
+        elif isinstance(key, slice):
+            res = []
+            start, stop, stride = key.indices(m)
+            for i in range(start, stop, stride):
+                res.append(self._get_geom_item(i))
+            return type(self._parent)(res or None)
+        else:
+            raise TypeError("key must be an index or slice")
+
+
+class EmptyGeometry(BaseGeometry):
+    def __new__(self):
+        """Create an empty geometry."""
+        warn(
+            "The 'EmptyGeometry()' constructor to create an empty geometry is "
+            "deprecated, and will raise an error in the future. Use one of the "
+            "geometry subclasses instead, for example 'GeometryCollection()'.",
+            ShapelyDeprecationWarning,
+            stacklevel=2,
+        )
+        return shapely.from_wkt("GEOMETRYCOLLECTION EMPTY")

billinglayer/python/shapely/geometry/collection.py

@@ -0,0 +1,58 @@
+"""Multi-part collections of geometries
+"""
+
+import shapely
+from shapely.geometry.base import BaseGeometry, BaseMultipartGeometry
+
+
+class GeometryCollection(BaseMultipartGeometry):
+    """
+    A collection of one or more geometries that may contain more than one type
+    of geometry.
+
+    Parameters
+    ----------
+    geoms : list
+        A list of shapely geometry instances, which may be of varying
+        geometry types.
+
+    Attributes
+    ----------
+    geoms : sequence
+        A sequence of Shapely geometry instances
+
+    Examples
+    --------
+    Create a GeometryCollection with a Point and a LineString
+
+    >>> from shapely import LineString, Point
+    >>> p = Point(51, -1)
+    >>> l = LineString([(52, -1), (49, 2)])
+    >>> gc = GeometryCollection([p, l])
+    """
+
+    __slots__ = []
+
+    def __new__(self, geoms=None):
+        if not geoms:
+            # TODO better empty constructor
+            return shapely.from_wkt("GEOMETRYCOLLECTION EMPTY")
+        if isinstance(geoms, BaseGeometry):
+            # TODO(shapely-2.0) do we actually want to split Multi-part geometries?
+            # this is needed for the split() tests
+            if hasattr(geoms, "geoms"):
+                geoms = geoms.geoms
+            else:
+                geoms = [geoms]
+
+        return shapely.geometrycollections(geoms)
+
+    @property
+    def __geo_interface__(self):
+        geometries = []
+        for geom in self.geoms:
+            geometries.append(geom.__geo_interface__)
+        return dict(type="GeometryCollection", geometries=geometries)
+
+
+shapely.lib.registry[7] = GeometryCollection

billinglayer/python/shapely/geometry/conftest.py

@@ -0,0 +1,10 @@
+"""Autouse fixtures for doctests."""
+
+import pytest
+
+from .linestring import LineString
+
+
+@pytest.fixture(autouse=True)
+def add_linestring(doctest_namespace):
+    doctest_namespace["LineString"] = LineString

billinglayer/python/shapely/geometry/geo.py

@@ -0,0 +1,136 @@
+"""
+Geometry factories based on the geo interface
+"""
+import numpy as np
+
+from shapely.errors import GeometryTypeError
+
+from .collection import GeometryCollection
+from .linestring import LineString
+from .multilinestring import MultiLineString
+from .multipoint import MultiPoint
+from .multipolygon import MultiPolygon
+from .point import Point
+from .polygon import LinearRing, Polygon
+
+
+def _is_coordinates_empty(coordinates):
+    """Helper to identify if coordinates or subset of coordinates are empty"""
+
+    if coordinates is None:
+        return True
+
+    if isinstance(coordinates, (list, tuple, np.ndarray)):
+        if len(coordinates) == 0:
+            return True
+        return all(map(_is_coordinates_empty, coordinates))
+    else:
+        return False
+
+
+def _empty_shape_for_no_coordinates(geom_type):
+    """Return empty counterpart for geom_type"""
+    if geom_type == "point":
+        return Point()
+    elif geom_type == "multipoint":
+        return MultiPoint()
+    elif geom_type == "linestring":
+        return LineString()
+    elif geom_type == "multilinestring":
+        return MultiLineString()
+    elif geom_type == "polygon":
+        return Polygon()
+    elif geom_type == "multipolygon":
+        return MultiPolygon()
+    else:
+        raise GeometryTypeError(f"Unknown geometry type: {geom_type!r}")
+
+
+def box(minx, miny, maxx, maxy, ccw=True):
+    """Returns a rectangular polygon with configurable normal vector"""
+    coords = [(maxx, miny), (maxx, maxy), (minx, maxy), (minx, miny)]
+    if not ccw:
+        coords = coords[::-1]
+    return Polygon(coords)
+
+
+def shape(context):
+    """
+    Returns a new, independent geometry with coordinates *copied* from the
+    context. Changes to the original context will not be reflected in the
+    geometry object.
+
+    Parameters
+    ----------
+    context :
+        a GeoJSON-like dict, which provides a "type" member describing the type
+        of the geometry and "coordinates" member providing a list of coordinates,
+        or an object which implements __geo_interface__.
+
+    Returns
+    -------
+    Geometry object
+
+    Examples
+    --------
+    Create a Point from GeoJSON, and then create a copy using __geo_interface__.
+
+    >>> context = {'type': 'Point', 'coordinates': [0, 1]}
+    >>> geom = shape(context)
+    >>> geom.geom_type == 'Point'
+    True
+    >>> geom.wkt
+    'POINT (0 1)'
+    >>> geom2 = shape(geom)
+    >>> geom == geom2
+    True
+    """
+    if hasattr(context, "__geo_interface__"):
+        ob = context.__geo_interface__
+    else:
+        ob = context
+    geom_type = ob.get("type").lower()
+    if "coordinates" in ob and _is_coordinates_empty(ob["coordinates"]):
+        return _empty_shape_for_no_coordinates(geom_type)
+    elif geom_type == "point":
+        return Point(ob["coordinates"])
+    elif geom_type == "linestring":
+        return LineString(ob["coordinates"])
+    elif geom_type == "linearring":
+        return LinearRing(ob["coordinates"])
+    elif geom_type == "polygon":
+        return Polygon(ob["coordinates"][0], ob["coordinates"][1:])
+    elif geom_type == "multipoint":
+        return MultiPoint(ob["coordinates"])
+    elif geom_type == "multilinestring":
+        return MultiLineString(ob["coordinates"])
+    elif geom_type == "multipolygon":
+        return MultiPolygon([[c[0], c[1:]] for c in ob["coordinates"]])
+    elif geom_type == "geometrycollection":
+        geoms = [shape(g) for g in ob.get("geometries", [])]
+        return GeometryCollection(geoms)
+    else:
+        raise GeometryTypeError(f"Unknown geometry type: {geom_type!r}")
+
+
+def mapping(ob):
+    """
+    Returns a GeoJSON-like mapping from a Geometry or any
+    object which implements __geo_interface__
+
+    Parameters
+    ----------
+    ob :
+        An object which implements __geo_interface__.
+
+    Returns
+    -------
+    dict
+
+    Examples
+    --------
+    >>> pt = Point(0, 0)
+    >>> mapping(pt)
+    {'type': 'Point', 'coordinates': (0.0, 0.0)}
+    """
+    return ob.__geo_interface__

billinglayer/python/shapely/geometry/linestring.py

@@ -0,0 +1,188 @@
+"""Line strings and related utilities
+"""
+import numpy as np
+
+import shapely
+from shapely.geometry.base import BaseGeometry, JOIN_STYLE
+from shapely.geometry.point import Point
+
+__all__ = ["LineString"]
+
+
+class LineString(BaseGeometry):
+    """
+    A geometry type composed of one or more line segments.
+
+    A LineString is a one-dimensional feature and has a non-zero length but
+    zero area. It may approximate a curve and need not be straight. Unlike a
+    LinearRing, a LineString is not closed.
+
+    Parameters
+    ----------
+    coordinates : sequence
+        A sequence of (x, y, [,z]) numeric coordinate pairs or triples, or
+        an array-like with shape (N, 2) or (N, 3).
+        Also can be a sequence of Point objects.
+
+    Examples
+    --------
+    Create a LineString with two segments
+
+    >>> a = LineString([[0, 0], [1, 0], [1, 1]])
+    >>> a.length
+    2.0
+    """
+
+    __slots__ = []
+
+    def __new__(self, coordinates=None):
+        if coordinates is None:
+            # empty geometry
+            # TODO better constructor
+            return shapely.from_wkt("LINESTRING EMPTY")
+        elif isinstance(coordinates, LineString):
+            if type(coordinates) == LineString:
+                # return original objects since geometries are immutable
+                return coordinates
+            else:
+                # LinearRing
+                # TODO convert LinearRing to LineString more directly
+                coordinates = coordinates.coords
+        else:
+            if hasattr(coordinates, "__array__"):
+                coordinates = np.asarray(coordinates)
+            if isinstance(coordinates, np.ndarray) and np.issubdtype(
+                coordinates.dtype, np.number
+            ):
+                pass
+            else:
+                # check coordinates on points
+                def _coords(o):
+                    if isinstance(o, Point):
+                        return o.coords[0]
+                    else:
+                        return [float(c) for c in o]
+
+                coordinates = [_coords(o) for o in coordinates]
+
+        if len(coordinates) == 0:
+            # empty geometry
+            # TODO better constructor + should shapely.linestrings handle this?
+            return shapely.from_wkt("LINESTRING EMPTY")
+
+        geom = shapely.linestrings(coordinates)
+        if not isinstance(geom, LineString):
+            raise ValueError("Invalid values passed to LineString constructor")
+        return geom
+
+    @property
+    def __geo_interface__(self):
+        return {"type": "LineString", "coordinates": tuple(self.coords)}
+
+    def svg(self, scale_factor=1.0, stroke_color=None, opacity=None):
+        """Returns SVG polyline element for the LineString geometry.
+
+        Parameters
+        ==========
+        scale_factor : float
+            Multiplication factor for the SVG stroke-width.  Default is 1.
+        stroke_color : str, optional
+            Hex string for stroke color. Default is to use "#66cc99" if
+            geometry is valid, and "#ff3333" if invalid.
+        opacity : float
+            Float number between 0 and 1 for color opacity. Default value is 0.8
+        """
+        if self.is_empty:
+            return "<g />"
+        if stroke_color is None:
+            stroke_color = "#66cc99" if self.is_valid else "#ff3333"
+        if opacity is None:
+            opacity = 0.8
+        pnt_format = " ".join(["{},{}".format(*c) for c in self.coords])
+        return (
+            '<polyline fill="none" stroke="{2}" stroke-width="{1}" '
+            'points="{0}" opacity="{3}" />'
+        ).format(pnt_format, 2.0 * scale_factor, stroke_color, opacity)
+
+    @property
+    def xy(self):
+        """Separate arrays of X and Y coordinate values
+
+        Example:
+
+          >>> x, y = LineString([(0, 0), (1, 1)]).xy
+          >>> list(x)
+          [0.0, 1.0]
+          >>> list(y)
+          [0.0, 1.0]
+        """
+        return self.coords.xy
+
+    def offset_curve(
+        self,
+        distance,
+        quad_segs=16,
+        join_style=JOIN_STYLE.round,
+        mitre_limit=5.0,
+    ):
+        """Returns a LineString or MultiLineString geometry at a distance from
+        the object on its right or its left side.
+
+        The side is determined by the sign of the `distance` parameter
+        (negative for right side offset, positive for left side offset). The
+        resolution of the buffer around each vertex of the object increases
+        by increasing the `quad_segs` keyword parameter.
+
+        The join style is for outside corners between line segments. Accepted
+        values are JOIN_STYLE.round (1), JOIN_STYLE.mitre (2), and
+        JOIN_STYLE.bevel (3).
+
+        The mitre ratio limit is used for very sharp corners. It is the ratio
+        of the distance from the corner to the end of the mitred offset corner.
+        When two line segments meet at a sharp angle, a miter join will extend
+        far beyond the original geometry. To prevent unreasonable geometry, the
+        mitre limit allows controlling the maximum length of the join corner.
+        Corners with a ratio which exceed the limit will be beveled.
+
+        Note: the behaviour regarding orientation of the resulting line
+        depends on the GEOS version. With GEOS < 3.11, the line retains the
+        same direction for a left offset (positive distance) or has reverse
+        direction for a right offset (negative distance), and this behaviour
+        was documented as such in previous Shapely versions. Starting with
+        GEOS 3.11, the function tries to preserve the orientation of the
+        original line.
+        """
+        if mitre_limit == 0.0:
+            raise ValueError("Cannot compute offset from zero-length line segment")
+        elif not np.isfinite(distance):
+            raise ValueError("offset_curve distance must be finite")
+        return shapely.offset_curve(self, distance, quad_segs, join_style, mitre_limit)
+
+    def parallel_offset(
+        self,
+        distance,
+        side="right",
+        resolution=16,
+        join_style=JOIN_STYLE.round,
+        mitre_limit=5.0,
+    ):
+        """
+        Alternative method to :meth:`offset_curve` method.
+
+        Older alternative method to the :meth:`offset_curve` method, but uses
+        ``resolution`` instead of ``quad_segs`` and a ``side`` keyword
+        ('left' or 'right') instead of sign of the distance. This method is
+        kept for backwards compatibility for now, but is is recommended to
+        use :meth:`offset_curve` instead.
+        """
+        if side == "right":
+            distance *= -1
+        return self.offset_curve(
+            distance,
+            quad_segs=resolution,
+            join_style=join_style,
+            mitre_limit=mitre_limit,
+        )
+
+
+shapely.lib.registry[1] = LineString

billinglayer/python/shapely/geometry/multilinestring.py

@@ -0,0 +1,93 @@
+"""Collections of linestrings and related utilities
+"""
+
+import shapely
+from shapely.errors import EmptyPartError
+from shapely.geometry import linestring
+from shapely.geometry.base import BaseMultipartGeometry
+
+__all__ = ["MultiLineString"]
+
+
+class MultiLineString(BaseMultipartGeometry):
+    """
+    A collection of one or more LineStrings.
+
+    A MultiLineString has non-zero length and zero area.
+
+    Parameters
+    ----------
+    lines : sequence
+        A sequence LineStrings, or a sequence of line-like coordinate
+        sequences or array-likes (see accepted input for LineString).
+
+    Attributes
+    ----------
+    geoms : sequence
+        A sequence of LineStrings
+
+    Examples
+    --------
+    Construct a MultiLineString containing two LineStrings.
+
+    >>> lines = MultiLineString([[[0, 0], [1, 2]], [[4, 4], [5, 6]]])
+    """
+
+    __slots__ = []
+
+    def __new__(self, lines=None):
+        if not lines:
+            # allow creation of empty multilinestrings, to support unpickling
+            # TODO better empty constructor
+            return shapely.from_wkt("MULTILINESTRING EMPTY")
+        elif isinstance(lines, MultiLineString):
+            return lines
+
+        lines = getattr(lines, "geoms", lines)
+        m = len(lines)
+        subs = []
+        for i in range(m):
+            line = linestring.LineString(lines[i])
+            if line.is_empty:
+                raise EmptyPartError(
+                    "Can't create MultiLineString with empty component"
+                )
+            subs.append(line)
+
+        if len(lines) == 0:
+            return shapely.from_wkt("MULTILINESTRING EMPTY")
+
+        return shapely.multilinestrings(subs)
+
+    @property
+    def __geo_interface__(self):
+        return {
+            "type": "MultiLineString",
+            "coordinates": tuple(tuple(c for c in g.coords) for g in self.geoms),
+        }
+
+    def svg(self, scale_factor=1.0, stroke_color=None, opacity=None):
+        """Returns a group of SVG polyline elements for the LineString geometry.
+
+        Parameters
+        ==========
+        scale_factor : float
+            Multiplication factor for the SVG stroke-width.  Default is 1.
+        stroke_color : str, optional
+            Hex string for stroke color. Default is to use "#66cc99" if
+            geometry is valid, and "#ff3333" if invalid.
+        opacity : float
+            Float number between 0 and 1 for color opacity. Default value is 0.8
+        """
+        if self.is_empty:
+            return "<g />"
+        if stroke_color is None:
+            stroke_color = "#66cc99" if self.is_valid else "#ff3333"
+        return (
+            "<g>"
+            + "".join(p.svg(scale_factor, stroke_color, opacity) for p in self.geoms)
+            + "</g>"
+        )
+
+
+shapely.lib.registry[5] = MultiLineString

billinglayer/python/shapely/geometry/multipoint.py

@@ -0,0 +1,95 @@
+"""Collections of points and related utilities
+"""
+
+import shapely
+from shapely.errors import EmptyPartError
+from shapely.geometry import point
+from shapely.geometry.base import BaseMultipartGeometry
+
+__all__ = ["MultiPoint"]
+
+
+class MultiPoint(BaseMultipartGeometry):
+    """
+    A collection of one or more Points.
+
+    A MultiPoint has zero area and zero length.
+
+    Parameters
+    ----------
+    points : sequence
+        A sequence of Points, or a sequence of (x, y [,z]) numeric coordinate
+        pairs or triples, or an array-like of shape (N, 2) or (N, 3).
+
+    Attributes
+    ----------
+    geoms : sequence
+        A sequence of Points
+
+    Examples
+    --------
+    Construct a MultiPoint containing two Points
+
+    >>> from shapely import Point
+    >>> ob = MultiPoint([[0.0, 0.0], [1.0, 2.0]])
+    >>> len(ob.geoms)
+    2
+    >>> type(ob.geoms[0]) == Point
+    True
+    """
+
+    __slots__ = []
+
+    def __new__(self, points=None):
+        if points is None:
+            # allow creation of empty multipoints, to support unpickling
+            # TODO better empty constructor
+            return shapely.from_wkt("MULTIPOINT EMPTY")
+        elif isinstance(points, MultiPoint):
+            return points
+
+        m = len(points)
+        subs = []
+        for i in range(m):
+            p = point.Point(points[i])
+            if p.is_empty:
+                raise EmptyPartError("Can't create MultiPoint with empty component")
+            subs.append(p)
+
+        if len(points) == 0:
+            return shapely.from_wkt("MULTIPOINT EMPTY")
+
+        return shapely.multipoints(subs)
+
+    @property
+    def __geo_interface__(self):
+        return {
+            "type": "MultiPoint",
+            "coordinates": tuple(g.coords[0] for g in self.geoms),
+        }
+
+    def svg(self, scale_factor=1.0, fill_color=None, opacity=None):
+        """Returns a group of SVG circle elements for the MultiPoint geometry.
+
+        Parameters
+        ==========
+        scale_factor : float
+            Multiplication factor for the SVG circle diameters.  Default is 1.
+        fill_color : str, optional
+            Hex string for fill color. Default is to use "#66cc99" if
+            geometry is valid, and "#ff3333" if invalid.
+        opacity : float
+            Float number between 0 and 1 for color opacity. Default value is 0.6
+        """
+        if self.is_empty:
+            return "<g />"
+        if fill_color is None:
+            fill_color = "#66cc99" if self.is_valid else "#ff3333"
+        return (
+            "<g>"
+            + "".join(p.svg(scale_factor, fill_color, opacity) for p in self.geoms)
+            + "</g>"
+        )
+
+
+shapely.lib.registry[4] = MultiPoint

billinglayer/python/shapely/geometry/multipolygon.py

@@ -0,0 +1,123 @@
+"""Collections of polygons and related utilities
+"""
+
+import shapely
+from shapely.geometry import polygon
+from shapely.geometry.base import BaseMultipartGeometry
+
+__all__ = ["MultiPolygon"]
+
+
+class MultiPolygon(BaseMultipartGeometry):
+    """
+    A collection of one or more Polygons.
+
+    If component polygons overlap the collection is invalid and some
+    operations on it may fail.
+
+    Parameters
+    ----------
+    polygons : sequence
+        A sequence of Polygons, or a sequence of (shell, holes) tuples
+        where shell is the sequence representation of a linear ring
+        (see LinearRing) and holes is a sequence of such linear rings.
+
+    Attributes
+    ----------
+    geoms : sequence
+        A sequence of `Polygon` instances
+
+    Examples
+    --------
+    Construct a MultiPolygon from a sequence of coordinate tuples
+
+    >>> from shapely import Polygon
+    >>> ob = MultiPolygon([
+    ...     (
+    ...     ((0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0)),
+    ...     [((0.1,0.1), (0.1,0.2), (0.2,0.2), (0.2,0.1))]
+    ...     )
+    ... ])
+    >>> len(ob.geoms)
+    1
+    >>> type(ob.geoms[0]) == Polygon
+    True
+    """
+
+    __slots__ = []
+
+    def __new__(self, polygons=None):
+        if not polygons:
+            # allow creation of empty multipolygons, to support unpickling
+            # TODO better empty constructor
+            return shapely.from_wkt("MULTIPOLYGON EMPTY")
+        elif isinstance(polygons, MultiPolygon):
+            return polygons
+
+        polygons = getattr(polygons, "geoms", polygons)
+        polygons = [
+            p
+            for p in polygons
+            if p and not (isinstance(p, polygon.Polygon) and p.is_empty)
+        ]
+
+        L = len(polygons)
+
+        # Bail immediately if we have no input points.
+        if L == 0:
+            return shapely.from_wkt("MULTIPOLYGON EMPTY")
+
+        # This function does not accept sequences of MultiPolygons: there is
+        # no implicit flattening.
+        if isinstance(polygons[0], MultiPolygon):
+            raise ValueError("Sequences of multi-polygons are not valid arguments")
+
+        subs = []
+        for i in range(L):
+            ob = polygons[i]
+            if not isinstance(ob, polygon.Polygon):
+                shell = ob[0]
+                holes = ob[1]
+                p = polygon.Polygon(shell, holes)
+            else:
+                p = polygon.Polygon(ob)
+            subs.append(p)
+
+        return shapely.multipolygons(subs)
+
+    @property
+    def __geo_interface__(self):
+        allcoords = []
+        for geom in self.geoms:
+            coords = []
+            coords.append(tuple(geom.exterior.coords))
+            for hole in geom.interiors:
+                coords.append(tuple(hole.coords))
+            allcoords.append(tuple(coords))
+        return {"type": "MultiPolygon", "coordinates": allcoords}
+
+    def svg(self, scale_factor=1.0, fill_color=None, opacity=None):
+        """Returns group of SVG path elements for the MultiPolygon geometry.
+
+        Parameters
+        ==========
+        scale_factor : float
+            Multiplication factor for the SVG stroke-width.  Default is 1.
+        fill_color : str, optional
+            Hex string for fill color. Default is to use "#66cc99" if
+            geometry is valid, and "#ff3333" if invalid.
+        opacity : float
+            Float number between 0 and 1 for color opacity. Default value is 0.6
+        """
+        if self.is_empty:
+            return "<g />"
+        if fill_color is None:
+            fill_color = "#66cc99" if self.is_valid else "#ff3333"
+        return (
+            "<g>"
+            + "".join(p.svg(scale_factor, fill_color, opacity) for p in self.geoms)
+            + "</g>"
+        )
+
+
+shapely.lib.registry[6] = MultiPolygon

billinglayer/python/shapely/geometry/point.py

@@ -0,0 +1,145 @@
+"""Points and related utilities
+"""
+import numpy as np
+
+import shapely
+from shapely.errors import DimensionError
+from shapely.geometry.base import BaseGeometry
+
+__all__ = ["Point"]
+
+
+class Point(BaseGeometry):
+    """
+    A geometry type that represents a single coordinate with
+    x,y and possibly z values.
+
+    A point is a zero-dimensional feature and has zero length and zero area.
+
+    Parameters
+    ----------
+    args : float, or sequence of floats
+        The coordinates can either be passed as a single parameter, or as
+        individual float values using multiple parameters:
+
+        1) 1 parameter: a sequence or array-like of with 2 or 3 values.
+        2) 2 or 3 parameters (float): x, y, and possibly z.
+
+    Attributes
+    ----------
+    x, y, z : float
+        Coordinate values
+
+    Examples
+    --------
+    Constructing the Point using separate parameters for x and y:
+
+    >>> p = Point(1.0, -1.0)
+
+    Constructing the Point using a list of x, y coordinates:
+
+    >>> p = Point([1.0, -1.0])
+    >>> print(p)
+    POINT (1 -1)
+    >>> p.y
+    -1.0
+    >>> p.x
+    1.0
+    """
+
+    __slots__ = []
+
+    def __new__(self, *args):
+        if len(args) == 0:
+            # empty geometry
+            # TODO better constructor
+            return shapely.from_wkt("POINT EMPTY")
+        elif len(args) > 3:
+            raise TypeError(f"Point() takes at most 3 arguments ({len(args)} given)")
+        elif len(args) == 1:
+            coords = args[0]
+            if isinstance(coords, Point):
+                return coords
+
+            # Accept either (x, y) or [(x, y)]
+            if not hasattr(coords, "__getitem__"):  # generators
+                coords = list(coords)
+            coords = np.asarray(coords).squeeze()
+        else:
+            # 2 or 3 args
+            coords = np.array(args).squeeze()
+
+        if coords.ndim > 1:
+            raise ValueError(
+                f"Point() takes only scalar or 1-size vector arguments, got {args}"
+            )
+        if not np.issubdtype(coords.dtype, np.number):
+            coords = [float(c) for c in coords]
+        geom = shapely.points(coords)
+        if not isinstance(geom, Point):
+            raise ValueError("Invalid values passed to Point constructor")
+        return geom
+
+    # Coordinate getters and setters
+
+    @property
+    def x(self):
+        """Return x coordinate."""
+        return shapely.get_x(self)
+
+    @property
+    def y(self):
+        """Return y coordinate."""
+        return shapely.get_y(self)
+
+    @property
+    def z(self):
+        """Return z coordinate."""
+        if not shapely.has_z(self):
+            raise DimensionError("This point has no z coordinate.")
+        # return shapely.get_z(self) -> get_z only supported for GEOS 3.7+
+        return self.coords[0][2]
+
+    @property
+    def __geo_interface__(self):
+        return {"type": "Point", "coordinates": self.coords[0]}
+
+    def svg(self, scale_factor=1.0, fill_color=None, opacity=None):
+        """Returns SVG circle element for the Point geometry.
+
+        Parameters
+        ==========
+        scale_factor : float
+            Multiplication factor for the SVG circle diameter.  Default is 1.
+        fill_color : str, optional
+            Hex string for fill color. Default is to use "#66cc99" if
+            geometry is valid, and "#ff3333" if invalid.
+        opacity : float
+            Float number between 0 and 1 for color opacity. Default value is 0.6
+        """
+        if self.is_empty:
+            return "<g />"
+        if fill_color is None:
+            fill_color = "#66cc99" if self.is_valid else "#ff3333"
+        if opacity is None:
+            opacity = 0.6
+        return (
+            '<circle cx="{0.x}" cy="{0.y}" r="{1}" '
+            'stroke="#555555" stroke-width="{2}" fill="{3}" opacity="{4}" />'
+        ).format(self, 3.0 * scale_factor, 1.0 * scale_factor, fill_color, opacity)
+
+    @property
+    def xy(self):
+        """Separate arrays of X and Y coordinate values
+
+        Example:
+          >>> x, y = Point(0, 0).xy
+          >>> list(x)
+          [0.0]
+          >>> list(y)
+          [0.0]
+        """
+        return self.coords.xy
+
+
+shapely.lib.registry[0] = Point

billinglayer/python/shapely/geometry/polygon.py

@@ -0,0 +1,326 @@
+"""Polygons and their linear ring components
+"""
+
+import numpy as np
+
+import shapely
+from shapely.algorithms.cga import is_ccw_impl, signed_area
+from shapely.errors import TopologicalError
+from shapely.geometry.base import BaseGeometry
+from shapely.geometry.linestring import LineString
+from shapely.geometry.point import Point
+
+__all__ = ["Polygon", "LinearRing"]
+
+
+def _unpickle_linearring(wkb):
+    linestring = shapely.from_wkb(wkb)
+    srid = shapely.get_srid(linestring)
+    linearring = shapely.linearrings(shapely.get_coordinates(linestring))
+    if srid:
+        linearring = shapely.set_srid(linearring, srid)
+    return linearring
+
+
+class LinearRing(LineString):
+    """
+    A geometry type composed of one or more line segments
+    that forms a closed loop.
+
+    A LinearRing is a closed, one-dimensional feature.
+    A LinearRing that crosses itself or touches itself at a single point is
+    invalid and operations on it may fail.
+
+    Parameters
+    ----------
+    coordinates : sequence
+        A sequence of (x, y [,z]) numeric coordinate pairs or triples, or
+        an array-like with shape (N, 2) or (N, 3).
+        Also can be a sequence of Point objects.
+
+    Notes
+    -----
+    Rings are automatically closed. There is no need to specify a final
+    coordinate pair identical to the first.
+
+    Examples
+    --------
+    Construct a square ring.
+
+    >>> ring = LinearRing( ((0, 0), (0, 1), (1 ,1 ), (1 , 0)) )
+    >>> ring.is_closed
+    True
+    >>> list(ring.coords)
+    [(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)]
+    >>> ring.length
+    4.0
+
+    """
+
+    __slots__ = []
+
+    def __new__(self, coordinates=None):
+        if coordinates is None:
+            # empty geometry
+            # TODO better way?
+            return shapely.from_wkt("LINEARRING EMPTY")
+        elif isinstance(coordinates, LineString):
+            if type(coordinates) == LinearRing:
+                # return original objects since geometries are immutable
+                return coordinates
+            elif not coordinates.is_valid:
+                raise TopologicalError("An input LineString must be valid.")
+            else:
+                # LineString
+                # TODO convert LineString to LinearRing more directly?
+                coordinates = coordinates.coords
+
+        else:
+            if hasattr(coordinates, "__array__"):
+                coordinates = np.asarray(coordinates)
+            if isinstance(coordinates, np.ndarray) and np.issubdtype(
+                coordinates.dtype, np.number
+            ):
+                pass
+            else:
+                # check coordinates on points
+                def _coords(o):
+                    if isinstance(o, Point):
+                        return o.coords[0]
+                    else:
+                        return [float(c) for c in o]
+
+                coordinates = np.array([_coords(o) for o in coordinates])
+                if not np.issubdtype(coordinates.dtype, np.number):
+                    # conversion of coords to 2D array failed, this might be due
+                    # to inconsistent coordinate dimensionality
+                    raise ValueError("Inconsistent coordinate dimensionality")
+
+        if len(coordinates) == 0:
+            # empty geometry
+            # TODO better constructor + should shapely.linearrings handle this?
+            return shapely.from_wkt("LINEARRING EMPTY")
+
+        geom = shapely.linearrings(coordinates)
+        if not isinstance(geom, LinearRing):
+            raise ValueError("Invalid values passed to LinearRing constructor")
+        return geom
+
+    @property
+    def __geo_interface__(self):
+        return {"type": "LinearRing", "coordinates": tuple(self.coords)}
+
+    def __reduce__(self):
+        """WKB doesn't differentiate between LineString and LinearRing so we
+        need to move the coordinate sequence into the correct geometry type"""
+        return (_unpickle_linearring, (shapely.to_wkb(self, include_srid=True),))
+
+    @property
+    def is_ccw(self):
+        """True is the ring is oriented counter clock-wise"""
+        return bool(is_ccw_impl()(self))
+
+    @property
+    def is_simple(self):
+        """True if the geometry is simple, meaning that any self-intersections
+        are only at boundary points, else False"""
+        return bool(shapely.is_simple(self))
+
+
+shapely.lib.registry[2] = LinearRing
+
+
+class InteriorRingSequence:
+
+    _parent = None
+    _ndim = None
+    _index = 0
+    _length = 0
+
+    def __init__(self, parent):
+        self._parent = parent
+        self._ndim = parent._ndim
+
+    def __iter__(self):
+        self._index = 0
+        self._length = self.__len__()
+        return self
+
+    def __next__(self):
+        if self._index < self._length:
+            ring = self._get_ring(self._index)
+            self._index += 1
+            return ring
+        else:
+            raise StopIteration
+
+    def __len__(self):
+        return shapely.get_num_interior_rings(self._parent)
+
+    def __getitem__(self, key):
+        m = self.__len__()
+        if isinstance(key, int):
+            if key + m < 0 or key >= m:
+                raise IndexError("index out of range")
+            if key < 0:
+                i = m + key
+            else:
+                i = key
+            return self._get_ring(i)
+        elif isinstance(key, slice):
+            res = []
+            start, stop, stride = key.indices(m)
+            for i in range(start, stop, stride):
+                res.append(self._get_ring(i))
+            return res
+        else:
+            raise TypeError("key must be an index or slice")
+
+    def _get_ring(self, i):
+        return shapely.get_interior_ring(self._parent, i)
+
+
+class Polygon(BaseGeometry):
+    """
+    A geometry type representing an area that is enclosed by a linear ring.
+
+    A polygon is a two-dimensional feature and has a non-zero area. It may
+    have one or more negative-space "holes" which are also bounded by linear
+    rings. If any rings cross each other, the feature is invalid and
+    operations on it may fail.
+
+    Parameters
+    ----------
+    shell : sequence
+        A sequence of (x, y [,z]) numeric coordinate pairs or triples, or
+        an array-like with shape (N, 2) or (N, 3).
+        Also can be a sequence of Point objects.
+    holes : sequence
+        A sequence of objects which satisfy the same requirements as the
+        shell parameters above
+
+    Attributes
+    ----------
+    exterior : LinearRing
+        The ring which bounds the positive space of the polygon.
+    interiors : sequence
+        A sequence of rings which bound all existing holes.
+
+    Examples
+    --------
+    Create a square polygon with no holes
+
+    >>> coords = ((0., 0.), (0., 1.), (1., 1.), (1., 0.), (0., 0.))
+    >>> polygon = Polygon(coords)
+    >>> polygon.area
+    1.0
+    """
+
+    __slots__ = []
+
+    def __new__(self, shell=None, holes=None):
+        if shell is None:
+            # empty geometry
+            # TODO better way?
+            return shapely.from_wkt("POLYGON EMPTY")
+        elif isinstance(shell, Polygon):
+            # return original objects since geometries are immutable
+            return shell
+        else:
+            shell = LinearRing(shell)
+
+        if holes is not None:
+            if len(holes) == 0:
+                # shapely constructor cannot handle holes=[]
+                holes = None
+            else:
+                holes = [LinearRing(ring) for ring in holes]
+
+        geom = shapely.polygons(shell, holes=holes)
+        if not isinstance(geom, Polygon):
+            raise ValueError("Invalid values passed to Polygon constructor")
+        return geom
+
+    @property
+    def exterior(self):
+        return shapely.get_exterior_ring(self)
+
+    @property
+    def interiors(self):
+        if self.is_empty:
+            return []
+        return InteriorRingSequence(self)
+
+    @property
+    def coords(self):
+        raise NotImplementedError(
+            "Component rings have coordinate sequences, but the polygon does not"
+        )
+
+    @property
+    def __geo_interface__(self):
+        if self.exterior == LinearRing():
+            coords = []
+        else:
+            coords = [tuple(self.exterior.coords)]
+            for hole in self.interiors:
+                coords.append(tuple(hole.coords))
+        return {"type": "Polygon", "coordinates": tuple(coords)}
+
+    def svg(self, scale_factor=1.0, fill_color=None, opacity=None):
+        """Returns SVG path element for the Polygon geometry.
+
+        Parameters
+        ==========
+        scale_factor : float
+            Multiplication factor for the SVG stroke-width.  Default is 1.
+        fill_color : str, optional
+            Hex string for fill color. Default is to use "#66cc99" if
+            geometry is valid, and "#ff3333" if invalid.
+        opacity : float
+            Float number between 0 and 1 for color opacity. Default value is 0.6
+        """
+        if self.is_empty:
+            return "<g />"
+        if fill_color is None:
+            fill_color = "#66cc99" if self.is_valid else "#ff3333"
+        if opacity is None:
+            opacity = 0.6
+        exterior_coords = [["{},{}".format(*c) for c in self.exterior.coords]]
+        interior_coords = [
+            ["{},{}".format(*c) for c in interior.coords] for interior in self.interiors
+        ]
+        path = " ".join(
+            [
+                "M {} L {} z".format(coords[0], " L ".join(coords[1:]))
+                for coords in exterior_coords + interior_coords
+            ]
+        )
+        return (
+            '<path fill-rule="evenodd" fill="{2}" stroke="#555555" '
+            'stroke-width="{0}" opacity="{3}" d="{1}" />'
+        ).format(2.0 * scale_factor, path, fill_color, opacity)
+
+    @classmethod
+    def from_bounds(cls, xmin, ymin, xmax, ymax):
+        """Construct a `Polygon()` from spatial bounds."""
+        return cls([(xmin, ymin), (xmin, ymax), (xmax, ymax), (xmax, ymin)])
+
+
+shapely.lib.registry[3] = Polygon
+
+
+def orient(polygon, sign=1.0):
+    s = float(sign)
+    rings = []
+    ring = polygon.exterior
+    if signed_area(ring) / s >= 0.0:
+        rings.append(ring)
+    else:
+        rings.append(list(ring.coords)[::-1])
+    for ring in polygon.interiors:
+        if signed_area(ring) / s <= 0.0:
+            rings.append(ring)
+        else:
+            rings.append(list(ring.coords)[::-1])
+    return Polygon(rings[0], rings[1:])

billinglayer/python/shapely/geos.py

@@ -0,0 +1,8 @@
+"""
+Proxies for libgeos, GEOS-specific exceptions, and utilities
+"""
+import shapely
+
+geos_version_string = shapely.geos_capi_version_string
+geos_version = shapely.geos_version
+geos_capi_version = shapely.geos_capi_version

billinglayer/python/shapely/io.py

@@ -0,0 +1,378 @@
+import numpy as np
+
+from . import lib
+from ._enum import ParamEnum
+
+# include ragged array functions here for reference documentation purpose
+from ._ragged_array import from_ragged_array, to_ragged_array
+from .decorators import requires_geos
+from .errors import UnsupportedGEOSVersionError
+
+__all__ = [
+    "from_geojson",
+    "from_ragged_array",
+    "from_wkb",
+    "from_wkt",
+    "to_geojson",
+    "to_ragged_array",
+    "to_wkb",
+    "to_wkt",
+]
+
+
+# Allowed options for handling WKB/WKT decoding errors
+# Note: cannot use standard constructor since "raise" is a keyword
+DecodingErrorOptions = ParamEnum(
+    "DecodingErrorOptions", {"ignore": 0, "warn": 1, "raise": 2}
+)
+
+WKBFlavorOptions = ParamEnum("WKBFlavorOptions", {"extended": 1, "iso": 2})
+
+
+def to_wkt(
+    geometry,
+    rounding_precision=6,
+    trim=True,
+    output_dimension=3,
+    old_3d=False,
+    **kwargs,
+):
+    """
+    Converts to the Well-Known Text (WKT) representation of a Geometry.
+
+    The Well-known Text format is defined in the `OGC Simple Features
+    Specification for SQL <https://www.opengeospatial.org/standards/sfs>`__.
+
+    The following limitations apply to WKT serialization:
+
+    - for GEOS <= 3.8 a multipoint with an empty sub-geometry will raise an exception
+    - for GEOS <= 3.8 empty geometries are always serialized to 2D
+    - for GEOS >= 3.9 only simple empty geometries can be 3D, collections are still
+      always 2D
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    rounding_precision : int, default 6
+        The rounding precision when writing the WKT string. Set to a value of
+        -1 to indicate the full precision.
+    trim : bool, default True
+        If True, trim unnecessary decimals (trailing zeros).
+    output_dimension : int, default 3
+        The output dimension for the WKT string. Supported values are 2 and 3.
+        Specifying 3 means that up to 3 dimensions will be written but 2D
+        geometries will still be represented as 2D in the WKT string.
+    old_3d : bool, default False
+        Enable old style 3D/4D WKT generation. By default, new style 3D/4D WKT
+        (ie. "POINT Z (10 20 30)") is returned, but with ``old_3d=True``
+        the WKT will be formatted in the style "POINT (10 20 30)".
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import Point
+    >>> to_wkt(Point(0, 0))
+    'POINT (0 0)'
+    >>> to_wkt(Point(0, 0), rounding_precision=3, trim=False)
+    'POINT (0.000 0.000)'
+    >>> to_wkt(Point(0, 0), rounding_precision=-1, trim=False)
+    'POINT (0.0000000000000000 0.0000000000000000)'
+    >>> to_wkt(Point(1, 2, 3), trim=True)
+    'POINT Z (1 2 3)'
+    >>> to_wkt(Point(1, 2, 3), trim=True, output_dimension=2)
+    'POINT (1 2)'
+    >>> to_wkt(Point(1, 2, 3), trim=True, old_3d=True)
+    'POINT (1 2 3)'
+
+    Notes
+    -----
+    The defaults differ from the default of the GEOS library. To mimic this,
+    use::
+
+        to_wkt(geometry, rounding_precision=-1, trim=False, output_dimension=2)
+
+    """
+    if not np.isscalar(rounding_precision):
+        raise TypeError("rounding_precision only accepts scalar values")
+    if not np.isscalar(trim):
+        raise TypeError("trim only accepts scalar values")
+    if not np.isscalar(output_dimension):
+        raise TypeError("output_dimension only accepts scalar values")
+    if not np.isscalar(old_3d):
+        raise TypeError("old_3d only accepts scalar values")
+
+    return lib.to_wkt(
+        geometry,
+        np.intc(rounding_precision),
+        np.bool_(trim),
+        np.intc(output_dimension),
+        np.bool_(old_3d),
+        **kwargs,
+    )
+
+
+def to_wkb(
+    geometry,
+    hex=False,
+    output_dimension=3,
+    byte_order=-1,
+    include_srid=False,
+    flavor="extended",
+    **kwargs,
+):
+    r"""
+    Converts to the Well-Known Binary (WKB) representation of a Geometry.
+
+    The Well-Known Binary format is defined in the `OGC Simple Features
+    Specification for SQL <https://www.opengeospatial.org/standards/sfs>`__.
+
+    The following limitations apply to WKB serialization:
+
+    - linearrings will be converted to linestrings
+    - a point with only NaN coordinates is converted to an empty point
+    - for GEOS <= 3.7, empty points are always serialized to 3D if
+      output_dimension=3, and to 2D if output_dimension=2
+    - for GEOS == 3.8, empty points are always serialized to 2D
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    hex : bool, default False
+        If true, export the WKB as a hexidecimal string. The default is to
+        return a binary bytes object.
+    output_dimension : int, default 3
+        The output dimension for the WKB. Supported values are 2 and 3.
+        Specifying 3 means that up to 3 dimensions will be written but 2D
+        geometries will still be represented as 2D in the WKB represenation.
+    byte_order : int, default -1
+        Defaults to native machine byte order (-1). Use 0 to force big endian
+        and 1 for little endian.
+    include_srid : bool, default False
+        If True, the SRID is be included in WKB (this is an extension
+        to the OGC WKB specification). Not allowed when flavor is "iso".
+    flavor : {"iso", "extended"}, default "extended"
+        Which flavor of WKB will be returned. The flavor determines how
+        extra dimensionality is encoded with the type number, and whether
+        SRID can be included in the WKB. ISO flavor is "more standard" for
+        3D output, and does not support SRID embedding.
+        Both flavors are equivalent when ``output_dimension=2`` (or with 2D
+        geometries) and ``include_srid=False``.
+        The `from_wkb` function can read both flavors.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import Point
+    >>> point = Point(1, 1)
+    >>> to_wkb(point, byte_order=1)
+    b'\x01\x01\x00\x00\x00\x00\x00\x00\x00\x00\x00\xf0?\x00\x00\x00\x00\x00\x00\xf0?'
+    >>> to_wkb(point, hex=True, byte_order=1)
+    '0101000000000000000000F03F000000000000F03F'
+    """
+    if not np.isscalar(hex):
+        raise TypeError("hex only accepts scalar values")
+    if not np.isscalar(output_dimension):
+        raise TypeError("output_dimension only accepts scalar values")
+    if not np.isscalar(byte_order):
+        raise TypeError("byte_order only accepts scalar values")
+    if not np.isscalar(include_srid):
+        raise TypeError("include_srid only accepts scalar values")
+    if not np.isscalar(flavor):
+        raise TypeError("flavor only accepts scalar values")
+    if lib.geos_version < (3, 10, 0) and flavor == "iso":
+        raise UnsupportedGEOSVersionError(
+            'The "iso" option requires at least GEOS 3.10.0'
+        )
+    if flavor == "iso" and include_srid:
+        raise ValueError('flavor="iso" and include_srid=True cannot be used together')
+    flavor = WKBFlavorOptions.get_value(flavor)
+
+    return lib.to_wkb(
+        geometry,
+        np.bool_(hex),
+        np.intc(output_dimension),
+        np.intc(byte_order),
+        np.bool_(include_srid),
+        np.intc(flavor),
+        **kwargs,
+    )
+
+
+@requires_geos("3.10.0")
+def to_geojson(geometry, indent=None, **kwargs):
+    """Converts to the GeoJSON representation of a Geometry.
+
+    The GeoJSON format is defined in the `RFC 7946 <https://geojson.org/>`__.
+    NaN (not-a-number) coordinates will be written as 'null'.
+
+    The following are currently unsupported:
+
+    - Geometries of type LINEARRING: these are output as 'null'.
+    - Three-dimensional geometries: the third dimension is ignored.
+
+    Parameters
+    ----------
+    geometry : str, bytes or array_like
+    indent : int, optional
+        If indent is a non-negative integer, then GeoJSON will be formatted.
+        An indent level of 0 will only insert newlines. None (the default)
+        selects the most compact representation.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import Point
+    >>> point = Point(1, 1)
+    >>> to_geojson(point)
+    '{"type":"Point","coordinates":[1.0,1.0]}'
+    >>> print(to_geojson(point, indent=2))
+    {
+      "type": "Point",
+      "coordinates": [
+          1.0,
+          1.0
+      ]
+    }
+    """
+    # GEOS Tickets:
+    # - handle linearrings: https://trac.osgeo.org/geos/ticket/1140
+    # - support 3D: https://trac.osgeo.org/geos/ticket/1141
+    if indent is None:
+        indent = -1
+    elif not np.isscalar(indent):
+        raise TypeError("indent only accepts scalar values")
+    elif indent < 0:
+        raise ValueError("indent cannot be negative")
+
+    return lib.to_geojson(geometry, np.intc(indent), **kwargs)
+
+
+def from_wkt(geometry, on_invalid="raise", **kwargs):
+    """
+    Creates geometries from the Well-Known Text (WKT) representation.
+
+    The Well-known Text format is defined in the `OGC Simple Features
+    Specification for SQL <https://www.opengeospatial.org/standards/sfs>`__.
+
+    Parameters
+    ----------
+    geometry : str or array_like
+        The WKT string(s) to convert.
+    on_invalid : {"raise", "warn", "ignore"}, default "raise"
+        - raise: an exception will be raised if WKT input geometries are invalid.
+        - warn: a warning will be raised and invalid WKT geometries will be
+          returned as ``None``.
+        - ignore: invalid WKT geometries will be returned as ``None`` without a warning.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from_wkt('POINT (0 0)')
+    <POINT (0 0)>
+    """
+    if not np.isscalar(on_invalid):
+        raise TypeError("on_invalid only accepts scalar values")
+
+    invalid_handler = np.uint8(DecodingErrorOptions.get_value(on_invalid))
+
+    return lib.from_wkt(geometry, invalid_handler, **kwargs)
+
+
+def from_wkb(geometry, on_invalid="raise", **kwargs):
+    r"""
+    Creates geometries from the Well-Known Binary (WKB) representation.
+
+    The Well-Known Binary format is defined in the `OGC Simple Features
+    Specification for SQL <https://www.opengeospatial.org/standards/sfs>`__.
+
+
+    Parameters
+    ----------
+    geometry : str or array_like
+        The WKB byte object(s) to convert.
+    on_invalid : {"raise", "warn", "ignore"}, default "raise"
+        - raise: an exception will be raised if a WKB input geometry is invalid.
+        - warn: a warning will be raised and invalid WKB geometries will be
+          returned as ``None``.
+        - ignore: invalid WKB geometries will be returned as ``None`` without a warning.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from_wkb(b'\x01\x01\x00\x00\x00\x00\x00\x00\x00\x00\x00\xf0?\x00\x00\x00\x00\x00\x00\xf0?')
+    <POINT (1 1)>
+    """
+
+    if not np.isscalar(on_invalid):
+        raise TypeError("on_invalid only accepts scalar values")
+
+    invalid_handler = np.uint8(DecodingErrorOptions.get_value(on_invalid))
+
+    # ensure the input has object dtype, to avoid numpy inferring it as a
+    # fixed-length string dtype (which removes trailing null bytes upon access
+    # of array elements)
+    geometry = np.asarray(geometry, dtype=object)
+    return lib.from_wkb(geometry, invalid_handler, **kwargs)
+
+
+@requires_geos("3.10.1")
+def from_geojson(geometry, on_invalid="raise", **kwargs):
+    """Creates geometries from GeoJSON representations (strings).
+
+    If a GeoJSON is a FeatureCollection, it is read as a single geometry
+    (with type GEOMETRYCOLLECTION). This may be unpacked using the ``pygeos.get_parts``.
+    Properties are not read.
+
+    The GeoJSON format is defined in `RFC 7946 <https://geojson.org/>`__.
+
+    The following are currently unsupported:
+
+    - Three-dimensional geometries: the third dimension is ignored.
+    - Geometries having 'null' in the coordinates.
+
+    Parameters
+    ----------
+    geometry : str, bytes or array_like
+        The GeoJSON string or byte object(s) to convert.
+    on_invalid : {"raise", "warn", "ignore"}, default "raise"
+        - raise: an exception will be raised if an input GeoJSON is invalid.
+        - warn: a warning will be raised and invalid input geometries will be
+          returned as ``None``.
+        - ignore: invalid input geometries will be returned as ``None`` without a warning.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    get_parts
+
+    Examples
+    --------
+    >>> from_geojson('{"type": "Point","coordinates": [1, 2]}')
+    <POINT (1 2)>
+    """
+    # GEOS Tickets:
+    # - support 3D: https://trac.osgeo.org/geos/ticket/1141
+    # - handle null coordinates: https://trac.osgeo.org/geos/ticket/1142
+    if not np.isscalar(on_invalid):
+        raise TypeError("on_invalid only accepts scalar values")
+
+    invalid_handler = np.uint8(DecodingErrorOptions.get_value(on_invalid))
+
+    # ensure the input has object dtype, to avoid numpy inferring it as a
+    # fixed-length string dtype (which removes trailing null bytes upon access
+    # of array elements)
+    geometry = np.asarray(geometry, dtype=object)
+
+    return lib.from_geojson(geometry, invalid_handler, **kwargs)

billinglayer/python/shapely/linear.py

@@ -0,0 +1,208 @@
+from . import lib
+from .decorators import multithreading_enabled
+from .errors import UnsupportedGEOSVersionError
+
+__all__ = [
+    "line_interpolate_point",
+    "line_locate_point",
+    "line_merge",
+    "shared_paths",
+    "shortest_line",
+]
+
+
+@multithreading_enabled
+def line_interpolate_point(line, distance, normalized=False, **kwargs):
+    """Returns a point interpolated at given distance on a line.
+
+    Parameters
+    ----------
+    line : Geometry or array_like
+        For multilinestrings or geometrycollections, the first geometry is taken
+        and the rest is ignored. This function raises a TypeError for non-linear
+        geometries. For empty linear geometries, empty points are returned.
+    distance : float or array_like
+        Negative values measure distance from the end of the line. Out-of-range
+        values will be clipped to the line endings.
+    normalized : bool, default False
+        If True, the distance is a fraction of the total
+        line length instead of the absolute distance.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point
+    >>> line = LineString([(0, 2), (0, 10)])
+    >>> line_interpolate_point(line, 2)
+    <POINT (0 4)>
+    >>> line_interpolate_point(line, 100)
+    <POINT (0 10)>
+    >>> line_interpolate_point(line, -2)
+    <POINT (0 8)>
+    >>> line_interpolate_point(line, [0.25, -0.25], normalized=True).tolist()
+    [<POINT (0 4)>, <POINT (0 8)>]
+    >>> line_interpolate_point(LineString(), 1)
+    <POINT EMPTY>
+    """
+    if normalized:
+        return lib.line_interpolate_point_normalized(line, distance)
+    else:
+        return lib.line_interpolate_point(line, distance)
+
+
+@multithreading_enabled
+def line_locate_point(line, other, normalized=False, **kwargs):
+    """Returns the distance to the line origin of given point.
+
+    If given point does not intersect with the line, the point will first be
+    projected onto the line after which the distance is taken.
+
+    Parameters
+    ----------
+    line : Geometry or array_like
+    point : Geometry or array_like
+    normalized : bool, default False
+        If True, the distance is a fraction of the total
+        line length instead of the absolute distance.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point
+    >>> line = LineString([(0, 2), (0, 10)])
+    >>> point = Point(4, 4)
+    >>> line_locate_point(line, point)
+    2.0
+    >>> line_locate_point(line, point, normalized=True)
+    0.25
+    >>> line_locate_point(line, Point(0, 18))
+    8.0
+    >>> line_locate_point(LineString(), point)
+    nan
+    """
+    if normalized:
+        return lib.line_locate_point_normalized(line, other)
+    else:
+        return lib.line_locate_point(line, other)
+
+
+@multithreading_enabled
+def line_merge(line, directed=False, **kwargs):
+    """Returns (Multi)LineStrings formed by combining the lines in a
+    MultiLineString.
+
+    Lines are joined together at their endpoints in case two lines are
+    intersecting. Lines are not joined when 3 or more lines are intersecting at
+    the endpoints. Line elements that cannot be joined are kept as is in the
+    resulting MultiLineString.
+
+    The direction of each merged LineString will be that of the majority of the
+    LineStrings from which it was derived. Except if ``directed=True`` is
+    specified, then the operation will not change the order of points within
+    lines and so only lines which can be joined with no change in direction
+    are merged.
+
+    Parameters
+    ----------
+    line : Geometry or array_like
+    directed : bool, default False
+        Only combine lines if possible without changing point order.
+        Requires GEOS >= 3.11.0
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import MultiLineString
+    >>> line_merge(MultiLineString([[(0, 2), (0, 10)], [(0, 10), (5, 10)]]))
+    <LINESTRING (0 2, 0 10, 5 10)>
+    >>> line_merge(MultiLineString([[(0, 2), (0, 10)], [(0, 11), (5, 10)]]))
+    <MULTILINESTRING ((0 2, 0 10), (0 11, 5 10))>
+    >>> line_merge(MultiLineString())
+    <GEOMETRYCOLLECTION EMPTY>
+    >>> line_merge(MultiLineString([[(0, 0), (1, 0)], [(0, 0), (3, 0)]]))
+    <LINESTRING (1 0, 0 0, 3 0)>
+    >>> line_merge(MultiLineString([[(0, 0), (1, 0)], [(0, 0), (3, 0)]]), directed=True)
+    <MULTILINESTRING ((0 0, 1 0), (0 0, 3 0))>
+    """
+    if directed:
+        if lib.geos_version < (3, 11, 0):
+            raise UnsupportedGEOSVersionError(
+                "'{}' requires at least GEOS {}.{}.{}.".format(
+                    "line_merge", *(3, 11, 0)
+                )
+            )
+        return lib.line_merge_directed(line, **kwargs)
+    return lib.line_merge(line, **kwargs)
+
+
+@multithreading_enabled
+def shared_paths(a, b, **kwargs):
+    """Returns the shared paths between geom1 and geom2.
+
+    Both geometries should be linestrings or arrays of linestrings.
+    A geometrycollection or array of geometrycollections is returned
+    with two elements in each geometrycollection. The first element is a
+    multilinestring containing shared paths with the same direction
+    for both inputs. The second element is a multilinestring containing
+    shared paths with the opposite direction for the two inputs.
+
+    Parameters
+    ----------
+    a : Geometry or array_like
+    b : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import LineString
+    >>> line1 = LineString([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
+    >>> line2 = LineString([(1, 0), (2, 0), (2, 1), (1, 1), (1, 0)])
+    >>> shared_paths(line1, line2).wkt
+    'GEOMETRYCOLLECTION (MULTILINESTRING EMPTY, MULTILINESTRING ((1 0, 1 1)))'
+    >>> line3 = LineString([(1, 1), (0, 1)])
+    >>> shared_paths(line1, line3).wkt
+    'GEOMETRYCOLLECTION (MULTILINESTRING ((1 1, 0 1)), MULTILINESTRING EMPTY)'
+    """
+    return lib.shared_paths(a, b, **kwargs)
+
+
+@multithreading_enabled
+def shortest_line(a, b, **kwargs):
+    """
+    Returns the shortest line between two geometries.
+
+    The resulting line consists of two points, representing the nearest
+    points between the geometry pair. The line always starts in the first
+    geometry `a` and ends in he second geometry `b`. The endpoints of the
+    line will not necessarily be existing vertices of the input geometries
+    `a` and `b`, but can also be a point along a line segment.
+
+    Parameters
+    ----------
+    a : Geometry or array_like
+    b : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    prepare : improve performance by preparing ``a`` (the first argument) (for GEOS>=3.9)
+
+    Examples
+    --------
+    >>> from shapely import LineString
+    >>> line1 = LineString([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
+    >>> line2 = LineString([(0, 3), (3, 0), (5, 3)])
+    >>> shortest_line(line1, line2)
+    <LINESTRING (1 1, 1.5 1.5)>
+    """
+    return lib.shortest_line(a, b, **kwargs)

billinglayer/python/shapely/measurement.py

@@ -0,0 +1,335 @@
+import warnings
+
+import numpy as np
+
+from . import lib
+from .decorators import multithreading_enabled, requires_geos
+
+__all__ = [
+    "area",
+    "distance",
+    "bounds",
+    "total_bounds",
+    "length",
+    "hausdorff_distance",
+    "frechet_distance",
+    "minimum_clearance",
+    "minimum_bounding_radius",
+]
+
+
+@multithreading_enabled
+def area(geometry, **kwargs):
+    """Computes the area of a (multi)polygon.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import MultiPolygon, Polygon
+    >>> polygon = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)])
+    >>> area(polygon)
+    100.0
+    >>> area(MultiPolygon([polygon, Polygon([(10, 10), (10, 20), (20, 20), (20, 10), (10, 10)])]))
+    200.0
+    >>> area(Polygon())
+    0.0
+    >>> area(None)
+    nan
+    """
+    return lib.area(geometry, **kwargs)
+
+
+@multithreading_enabled
+def distance(a, b, **kwargs):
+    """Computes the Cartesian distance between two geometries.
+
+    Parameters
+    ----------
+    a, b : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point, Polygon
+    >>> point = Point(0, 0)
+    >>> distance(Point(10, 0), point)
+    10.0
+    >>> distance(LineString([(1, 1), (1, -1)]), point)
+    1.0
+    >>> distance(Polygon([(3, 0), (5, 0), (5, 5), (3, 5), (3, 0)]), point)
+    3.0
+    >>> distance(Point(), point)
+    nan
+    >>> distance(None, point)
+    nan
+    """
+    return lib.distance(a, b, **kwargs)
+
+
+@multithreading_enabled
+def bounds(geometry, **kwargs):
+    """Computes the bounds (extent) of a geometry.
+
+    For each geometry these 4 numbers are returned: min x, min y, max x, max y.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point, Polygon
+    >>> bounds(Point(2, 3)).tolist()
+    [2.0, 3.0, 2.0, 3.0]
+    >>> bounds(LineString([(0, 0), (0, 2), (3, 2)])).tolist()
+    [0.0, 0.0, 3.0, 2.0]
+    >>> bounds(Polygon()).tolist()
+    [nan, nan, nan, nan]
+    >>> bounds(None).tolist()
+    [nan, nan, nan, nan]
+    """
+    # We need to provide the `out` argument here for compatibility with
+    # numpy < 1.16. See https://github.com/numpy/numpy/issues/14949
+    geometry_arr = np.asarray(geometry, dtype=np.object_)
+    out = np.empty(geometry_arr.shape + (4,), dtype="float64")
+    return lib.bounds(geometry_arr, out=out, **kwargs)
+
+
+def total_bounds(geometry, **kwargs):
+    """Computes the total bounds (extent) of the geometry.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Returns
+    -------
+    numpy ndarray of [xmin, ymin, xmax, ymax]
+
+    Examples
+    --------
+    >>> from shapely import LineString, Point, Polygon
+    >>> total_bounds(Point(2, 3)).tolist()
+    [2.0, 3.0, 2.0, 3.0]
+    >>> total_bounds([Point(2, 3), Point(4, 5)]).tolist()
+    [2.0, 3.0, 4.0, 5.0]
+    >>> total_bounds([
+    ...     LineString([(0, 1), (0, 2), (3, 2)]),
+    ...     LineString([(4, 4), (4, 6), (6, 7)])
+    ... ]).tolist()
+    [0.0, 1.0, 6.0, 7.0]
+    >>> total_bounds(Polygon()).tolist()
+    [nan, nan, nan, nan]
+    >>> total_bounds([Polygon(), Point(2, 3)]).tolist()
+    [2.0, 3.0, 2.0, 3.0]
+    >>> total_bounds(None).tolist()
+    [nan, nan, nan, nan]
+    """
+    b = bounds(geometry, **kwargs)
+    if b.ndim == 1:
+        return b
+
+    with warnings.catch_warnings():
+        # ignore 'All-NaN slice encountered' warnings
+        warnings.simplefilter("ignore", RuntimeWarning)
+        return np.array(
+            [
+                np.nanmin(b[..., 0]),
+                np.nanmin(b[..., 1]),
+                np.nanmax(b[..., 2]),
+                np.nanmax(b[..., 3]),
+            ]
+        )
+
+
+@multithreading_enabled
+def length(geometry, **kwargs):
+    """Computes the length of a (multi)linestring or polygon perimeter.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import LineString, MultiLineString, Polygon
+    >>> length(LineString([(0, 0), (0, 2), (3, 2)]))
+    5.0
+    >>> length(MultiLineString([
+    ...     LineString([(0, 0), (1, 0)]),
+    ...     LineString([(1, 0), (2, 0)])
+    ... ]))
+    2.0
+    >>> length(Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)]))
+    40.0
+    >>> length(LineString())
+    0.0
+    >>> length(None)
+    nan
+    """
+    return lib.length(geometry, **kwargs)
+
+
+@multithreading_enabled
+def hausdorff_distance(a, b, densify=None, **kwargs):
+    """Compute the discrete Hausdorff distance between two geometries.
+
+    The Hausdorff distance is a measure of similarity: it is the greatest
+    distance between any point in A and the closest point in B. The discrete
+    distance is an approximation of this metric: only vertices are considered.
+    The parameter 'densify' makes this approximation less coarse by splitting
+    the line segments between vertices before computing the distance.
+
+    Parameters
+    ----------
+    a, b : Geometry or array_like
+    densify : float or array_like, optional
+        The value of densify is required to be between 0 and 1.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import LineString
+    >>> line1 = LineString([(130, 0), (0, 0), (0, 150)])
+    >>> line2 = LineString([(10, 10), (10, 150), (130, 10)])
+    >>> hausdorff_distance(line1, line2)  # doctest: +ELLIPSIS
+    14.14...
+    >>> hausdorff_distance(line1, line2, densify=0.5)
+    70.0
+    >>> hausdorff_distance(line1, LineString())
+    nan
+    >>> hausdorff_distance(line1, None)
+    nan
+    """
+    if densify is None:
+        return lib.hausdorff_distance(a, b, **kwargs)
+    else:
+        return lib.hausdorff_distance_densify(a, b, densify, **kwargs)
+
+
+@requires_geos("3.7.0")
+@multithreading_enabled
+def frechet_distance(a, b, densify=None, **kwargs):
+    """Compute the discrete Fréchet distance between two geometries.
+
+    The Fréchet distance is a measure of similarity: it is the greatest
+    distance between any point in A and the closest point in B. The discrete
+    distance is an approximation of this metric: only vertices are considered.
+    The parameter 'densify' makes this approximation less coarse by splitting
+    the line segments between vertices before computing the distance.
+
+    Fréchet distance sweep continuously along their respective curves
+    and the direction of curves is significant. This makes it a better measure
+    of similarity than Hausdorff distance for curve or surface matching.
+
+    Parameters
+    ----------
+    a, b : Geometry or array_like
+    densify : float or array_like, optional
+        The value of densify is required to be between 0 and 1.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import LineString
+    >>> line1 = LineString([(0, 0), (100, 0)])
+    >>> line2 = LineString([(0, 0), (50, 50), (100, 0)])
+    >>> frechet_distance(line1, line2)  # doctest: +ELLIPSIS
+    70.71...
+    >>> frechet_distance(line1, line2, densify=0.5)
+    50.0
+    >>> frechet_distance(line1, LineString())
+    nan
+    >>> frechet_distance(line1, None)
+    nan
+    """
+    if densify is None:
+        return lib.frechet_distance(a, b, **kwargs)
+    return lib.frechet_distance_densify(a, b, densify, **kwargs)
+
+
+@requires_geos("3.6.0")
+@multithreading_enabled
+def minimum_clearance(geometry, **kwargs):
+    """Computes the Minimum Clearance distance.
+
+    A geometry's "minimum clearance" is the smallest distance by which
+    a vertex of the geometry could be moved to produce an invalid geometry.
+
+    If no minimum clearance exists for a geometry (for example, a single
+    point, or an empty geometry), infinity is returned.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    Examples
+    --------
+    >>> from shapely import Polygon
+    >>> polygon = Polygon([(0, 0), (0, 10), (5, 6), (10, 10), (10, 0), (5, 4), (0, 0)])
+    >>> minimum_clearance(polygon)
+    2.0
+    >>> minimum_clearance(Polygon())
+    inf
+    >>> minimum_clearance(None)
+    nan
+    """
+    return lib.minimum_clearance(geometry, **kwargs)
+
+
+@requires_geos("3.8.0")
+@multithreading_enabled
+def minimum_bounding_radius(geometry, **kwargs):
+    """Computes the radius of the minimum bounding circle that encloses an input geometry.
+
+    Parameters
+    ----------
+    geometry : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+
+    Examples
+    --------
+    >>> from shapely import GeometryCollection, LineString, MultiPoint, Point, Polygon
+    >>> minimum_bounding_radius(Polygon([(0, 5), (5, 10), (10, 5), (5, 0), (0, 5)]))
+    5.0
+    >>> minimum_bounding_radius(LineString([(1, 1), (1, 10)]))
+    4.5
+    >>> minimum_bounding_radius(MultiPoint([(2, 2), (4, 2)]))
+    1.0
+    >>> minimum_bounding_radius(Point(0, 1))
+    0.0
+    >>> minimum_bounding_radius(GeometryCollection())
+    0.0
+
+    See also
+    --------
+    minimum_bounding_circle
+    """
+    return lib.minimum_bounding_radius(geometry, **kwargs)

billinglayer/python/shapely/ops.py

@@ -0,0 +1,746 @@
+"""Support for various GEOS geometry operations
+"""
+
+from warnings import warn
+
+import shapely
+from shapely.algorithms.polylabel import polylabel  # noqa
+from shapely.errors import GeometryTypeError, ShapelyDeprecationWarning
+from shapely.geometry import (
+    GeometryCollection,
+    LineString,
+    MultiLineString,
+    MultiPoint,
+    Point,
+    Polygon,
+    shape,
+)
+from shapely.geometry.base import BaseGeometry, BaseMultipartGeometry
+from shapely.geometry.polygon import orient as orient_
+from shapely.prepared import prep
+
+__all__ = [
+    "cascaded_union",
+    "linemerge",
+    "operator",
+    "polygonize",
+    "polygonize_full",
+    "transform",
+    "unary_union",
+    "triangulate",
+    "voronoi_diagram",
+    "split",
+    "nearest_points",
+    "validate",
+    "snap",
+    "shared_paths",
+    "clip_by_rect",
+    "orient",
+    "substring",
+]
+
+
+class CollectionOperator:
+    def shapeup(self, ob):
+        if isinstance(ob, BaseGeometry):
+            return ob
+        else:
+            try:
+                return shape(ob)
+            except (ValueError, AttributeError):
+                return LineString(ob)
+
+    def polygonize(self, lines):
+        """Creates polygons from a source of lines
+
+        The source may be a MultiLineString, a sequence of LineString objects,
+        or a sequence of objects than can be adapted to LineStrings.
+        """
+        source = getattr(lines, "geoms", None) or lines
+        try:
+            source = iter(source)
+        except TypeError:
+            source = [source]
+        finally:
+            obs = [self.shapeup(line) for line in source]
+        collection = shapely.polygonize(obs)
+        return collection.geoms
+
+    def polygonize_full(self, lines):
+        """Creates polygons from a source of lines, returning the polygons
+        and leftover geometries.
+
+        The source may be a MultiLineString, a sequence of LineString objects,
+        or a sequence of objects than can be adapted to LineStrings.
+
+        Returns a tuple of objects: (polygons, cut edges, dangles, invalid ring
+        lines). Each are a geometry collection.
+
+        Dangles are edges which have one or both ends which are not incident on
+        another edge endpoint. Cut edges are connected at both ends but do not
+        form part of polygon. Invalid ring lines form rings which are invalid
+        (bowties, etc).
+        """
+        source = getattr(lines, "geoms", None) or lines
+        try:
+            source = iter(source)
+        except TypeError:
+            source = [source]
+        finally:
+            obs = [self.shapeup(line) for line in source]
+        return shapely.polygonize_full(obs)
+
+    def linemerge(self, lines, directed=False):
+        """Merges all connected lines from a source
+
+        The source may be a MultiLineString, a sequence of LineString objects,
+        or a sequence of objects than can be adapted to LineStrings.  Returns a
+        LineString or MultiLineString when lines are not contiguous.
+        """
+        source = None
+        if getattr(lines, "geom_type", None) == "MultiLineString":
+            source = lines
+        elif hasattr(lines, "geoms"):
+            # other Multi geometries
+            source = MultiLineString([ls.coords for ls in lines.geoms])
+        elif hasattr(lines, "__iter__"):
+            try:
+                source = MultiLineString([ls.coords for ls in lines])
+            except AttributeError:
+                source = MultiLineString(lines)
+        if source is None:
+            raise ValueError(f"Cannot linemerge {lines}")
+        return shapely.line_merge(source, directed=directed)
+
+    def cascaded_union(self, geoms):
+        """Returns the union of a sequence of geometries
+
+        .. deprecated:: 1.8
+            This function was superseded by :meth:`unary_union`.
+        """
+        warn(
+            "The 'cascaded_union()' function is deprecated. "
+            "Use 'unary_union()' instead.",
+            ShapelyDeprecationWarning,
+            stacklevel=2,
+        )
+        return shapely.union_all(geoms, axis=None)
+
+    def unary_union(self, geoms):
+        """Returns the union of a sequence of geometries
+
+        Usually used to convert a collection into the smallest set of polygons
+        that cover the same area.
+        """
+        return shapely.union_all(geoms, axis=None)
+
+
+operator = CollectionOperator()
+polygonize = operator.polygonize
+polygonize_full = operator.polygonize_full
+linemerge = operator.linemerge
+cascaded_union = operator.cascaded_union
+unary_union = operator.unary_union
+
+
+def triangulate(geom, tolerance=0.0, edges=False):
+    """Creates the Delaunay triangulation and returns a list of geometries
+
+    The source may be any geometry type. All vertices of the geometry will be
+    used as the points of the triangulation.
+
+    From the GEOS documentation:
+    tolerance is the snapping tolerance used to improve the robustness of
+    the triangulation computation. A tolerance of 0.0 specifies that no
+    snapping will take place.
+
+    If edges is False, a list of Polygons (triangles) will be returned.
+    Otherwise the list of LineString edges is returned.
+
+    """
+    collection = shapely.delaunay_triangles(geom, tolerance=tolerance, only_edges=edges)
+    return [g for g in collection.geoms]
+
+
+def voronoi_diagram(geom, envelope=None, tolerance=0.0, edges=False):
+    """
+    Constructs a Voronoi Diagram [1] from the given geometry.
+    Returns a list of geometries.
+
+    Parameters
+    ----------
+    geom: geometry
+        the input geometry whose vertices will be used to calculate
+        the final diagram.
+    envelope: geometry, None
+        clipping envelope for the returned diagram, automatically
+        determined if None. The diagram will be clipped to the larger
+        of this envelope or an envelope surrounding the sites.
+    tolerance: float, 0.0
+        sets the snapping tolerance used to improve the robustness
+        of the computation. A tolerance of 0.0 specifies that no
+        snapping will take place.
+    edges: bool, False
+        If False, return regions as polygons. Else, return only
+        edges e.g. LineStrings.
+
+    GEOS documentation can be found at [2]
+
+    Returns
+    -------
+    GeometryCollection
+        geometries representing the Voronoi regions.
+
+    Notes
+    -----
+    The tolerance `argument` can be finicky and is known to cause the
+    algorithm to fail in several cases. If you're using `tolerance`
+    and getting a failure, try removing it. The test cases in
+    tests/test_voronoi_diagram.py show more details.
+
+
+    References
+    ----------
+    [1] https://en.wikipedia.org/wiki/Voronoi_diagram
+    [2] https://geos.osgeo.org/doxygen/geos__c_8h_source.html  (line 730)
+    """
+    try:
+        result = shapely.voronoi_polygons(
+            geom, tolerance=tolerance, extend_to=envelope, only_edges=edges
+        )
+    except shapely.GEOSException as err:
+        errstr = "Could not create Voronoi Diagram with the specified inputs "
+        errstr += f"({err!s})."
+        if tolerance:
+            errstr += " Try running again with default tolerance value."
+        raise ValueError(errstr) from err
+
+    if result.geom_type != "GeometryCollection":
+        return GeometryCollection([result])
+    return result
+
+
+def validate(geom):
+    return shapely.is_valid_reason(geom)
+
+
+def transform(func, geom):
+    """Applies `func` to all coordinates of `geom` and returns a new
+    geometry of the same type from the transformed coordinates.
+
+    `func` maps x, y, and optionally z to output xp, yp, zp. The input
+    parameters may iterable types like lists or arrays or single values.
+    The output shall be of the same type. Scalars in, scalars out.
+    Lists in, lists out.
+
+    For example, here is an identity function applicable to both types
+    of input.
+
+      def id_func(x, y, z=None):
+          return tuple(filter(None, [x, y, z]))
+
+      g2 = transform(id_func, g1)
+
+    Using pyproj >= 2.1, this example will accurately project Shapely geometries:
+
+      import pyproj
+
+      wgs84 = pyproj.CRS('EPSG:4326')
+      utm = pyproj.CRS('EPSG:32618')
+
+      project = pyproj.Transformer.from_crs(wgs84, utm, always_xy=True).transform
+
+      g2 = transform(project, g1)
+
+    Note that the always_xy kwarg is required here as Shapely geometries only support
+    X,Y coordinate ordering.
+
+    Lambda expressions such as the one in
+
+      g2 = transform(lambda x, y, z=None: (x+1.0, y+1.0), g1)
+
+    also satisfy the requirements for `func`.
+    """
+    if geom.is_empty:
+        return geom
+    if geom.geom_type in ("Point", "LineString", "LinearRing", "Polygon"):
+
+        # First we try to apply func to x, y, z sequences. When func is
+        # optimized for sequences, this is the fastest, though zipping
+        # the results up to go back into the geometry constructors adds
+        # extra cost.
+        try:
+            if geom.geom_type in ("Point", "LineString", "LinearRing"):
+                return type(geom)(zip(*func(*zip(*geom.coords))))
+            elif geom.geom_type == "Polygon":
+                shell = type(geom.exterior)(zip(*func(*zip(*geom.exterior.coords))))
+                holes = list(
+                    type(ring)(zip(*func(*zip(*ring.coords))))
+                    for ring in geom.interiors
+                )
+                return type(geom)(shell, holes)
+
+        # A func that assumes x, y, z are single values will likely raise a
+        # TypeError, in which case we'll try again.
+        except TypeError:
+            if geom.geom_type in ("Point", "LineString", "LinearRing"):
+                return type(geom)([func(*c) for c in geom.coords])
+            elif geom.geom_type == "Polygon":
+                shell = type(geom.exterior)([func(*c) for c in geom.exterior.coords])
+                holes = list(
+                    type(ring)([func(*c) for c in ring.coords])
+                    for ring in geom.interiors
+                )
+                return type(geom)(shell, holes)
+
+    elif geom.geom_type.startswith("Multi") or geom.geom_type == "GeometryCollection":
+        return type(geom)([transform(func, part) for part in geom.geoms])
+    else:
+        raise GeometryTypeError(f"Type {geom.geom_type!r} not recognized")
+
+
+def nearest_points(g1, g2):
+    """Returns the calculated nearest points in the input geometries
+
+    The points are returned in the same order as the input geometries.
+    """
+    seq = shapely.shortest_line(g1, g2)
+    if seq is None:
+        if g1.is_empty:
+            raise ValueError("The first input geometry is empty")
+        else:
+            raise ValueError("The second input geometry is empty")
+
+    p1 = shapely.get_point(seq, 0)
+    p2 = shapely.get_point(seq, 1)
+    return (p1, p2)
+
+
+def snap(g1, g2, tolerance):
+    """Snap one geometry to another with a given tolerance
+
+    Vertices of the first geometry are snapped to vertices of the second
+    geometry. The resulting snapped geometry is returned. The input geometries
+    are not modified.
+
+    Parameters
+    ----------
+    g1 : geometry
+        The first geometry
+    g2 : geometry
+        The second geometry
+    tolerance : float
+        The snapping tolerance
+
+    Example
+    -------
+    >>> square = Polygon([(1,1), (2, 1), (2, 2), (1, 2), (1, 1)])
+    >>> line = LineString([(0,0), (0.8, 0.8), (1.8, 0.95), (2.6, 0.5)])
+    >>> result = snap(line, square, 0.5)
+    >>> result.wkt
+    'LINESTRING (0 0, 1 1, 2 1, 2.6 0.5)'
+    """
+    return shapely.snap(g1, g2, tolerance)
+
+
+def shared_paths(g1, g2):
+    """Find paths shared between the two given lineal geometries
+
+    Returns a GeometryCollection with two elements:
+     - First element is a MultiLineString containing shared paths with the
+       same direction for both inputs.
+     - Second element is a MultiLineString containing shared paths with the
+       opposite direction for the two inputs.
+
+    Parameters
+    ----------
+    g1 : geometry
+        The first geometry
+    g2 : geometry
+        The second geometry
+    """
+    if not isinstance(g1, LineString):
+        raise GeometryTypeError("First geometry must be a LineString")
+    if not isinstance(g2, LineString):
+        raise GeometryTypeError("Second geometry must be a LineString")
+    return shapely.shared_paths(g1, g2)
+
+
+class SplitOp:
+    @staticmethod
+    def _split_polygon_with_line(poly, splitter):
+        """Split a Polygon with a LineString"""
+        if not isinstance(poly, Polygon):
+            raise GeometryTypeError("First argument must be a Polygon")
+        if not isinstance(splitter, LineString):
+            raise GeometryTypeError("Second argument must be a LineString")
+
+        union = poly.boundary.union(splitter)
+
+        # greatly improves split performance for big geometries with many
+        # holes (the following contains checks) with minimal overhead
+        # for common cases
+        poly = prep(poly)
+
+        # some polygonized geometries may be holes, we do not want them
+        # that's why we test if the original polygon (poly) contains
+        # an inner point of polygonized geometry (pg)
+        return [
+            pg for pg in polygonize(union) if poly.contains(pg.representative_point())
+        ]
+
+    @staticmethod
+    def _split_line_with_line(line, splitter):
+        """Split a LineString with another (Multi)LineString or (Multi)Polygon"""
+
+        # if splitter is a polygon, pick it's boundary
+        if splitter.geom_type in ("Polygon", "MultiPolygon"):
+            splitter = splitter.boundary
+
+        if not isinstance(line, LineString):
+            raise GeometryTypeError("First argument must be a LineString")
+        if not isinstance(splitter, LineString) and not isinstance(
+            splitter, MultiLineString
+        ):
+            raise GeometryTypeError(
+                "Second argument must be either a LineString or a MultiLineString"
+            )
+
+        # |    s\l   | Interior | Boundary | Exterior |
+        # |----------|----------|----------|----------|
+        # | Interior |  0 or F  |    *     |    *     |   At least one of these two must be 0
+        # | Boundary |  0 or F  |    *     |    *     |   So either '0********' or '[0F]**0*****'
+        # | Exterior |    *     |    *     |    *     |   No overlapping interiors ('1********')
+        relation = splitter.relate(line)
+        if relation[0] == "1":
+            # The lines overlap at some segment (linear intersection of interiors)
+            raise ValueError("Input geometry segment overlaps with the splitter.")
+        elif relation[0] == "0" or relation[3] == "0":
+            # The splitter crosses or touches the line's interior --> return multilinestring from the split
+            return line.difference(splitter)
+        else:
+            # The splitter does not cross or touch the line's interior --> return collection with identity line
+            return [line]
+
+    @staticmethod
+    def _split_line_with_point(line, splitter):
+        """Split a LineString with a Point"""
+        if not isinstance(line, LineString):
+            raise GeometryTypeError("First argument must be a LineString")
+        if not isinstance(splitter, Point):
+            raise GeometryTypeError("Second argument must be a Point")
+
+        # check if point is in the interior of the line
+        if not line.relate_pattern(splitter, "0********"):
+            # point not on line interior --> return collection with single identity line
+            # (REASONING: Returning a list with the input line reference and creating a
+            # GeometryCollection at the general split function prevents unnecessary copying
+            # of linestrings in multipoint splitting function)
+            return [line]
+        elif line.coords[0] == splitter.coords[0]:
+            # if line is a closed ring the previous test doesn't behave as desired
+            return [line]
+
+        # point is on line, get the distance from the first point on line
+        distance_on_line = line.project(splitter)
+        coords = list(line.coords)
+        # split the line at the point and create two new lines
+        current_position = 0.0
+        for i in range(len(coords) - 1):
+            point1 = coords[i]
+            point2 = coords[i + 1]
+            dx = point1[0] - point2[0]
+            dy = point1[1] - point2[1]
+            segment_length = (dx**2 + dy**2) ** 0.5
+            current_position += segment_length
+            if distance_on_line == current_position:
+                # splitter is exactly on a vertex
+                return [LineString(coords[: i + 2]), LineString(coords[i + 1 :])]
+            elif distance_on_line < current_position:
+                # splitter is between two vertices
+                return [
+                    LineString(coords[: i + 1] + [splitter.coords[0]]),
+                    LineString([splitter.coords[0]] + coords[i + 1 :]),
+                ]
+        return [line]
+
+    @staticmethod
+    def _split_line_with_multipoint(line, splitter):
+        """Split a LineString with a MultiPoint"""
+
+        if not isinstance(line, LineString):
+            raise GeometryTypeError("First argument must be a LineString")
+        if not isinstance(splitter, MultiPoint):
+            raise GeometryTypeError("Second argument must be a MultiPoint")
+
+        chunks = [line]
+        for pt in splitter.geoms:
+            new_chunks = []
+            for chunk in filter(lambda x: not x.is_empty, chunks):
+                # add the newly split 2 lines or the same line if not split
+                new_chunks.extend(SplitOp._split_line_with_point(chunk, pt))
+            chunks = new_chunks
+
+        return chunks
+
+    @staticmethod
+    def split(geom, splitter):
+        """
+        Splits a geometry by another geometry and returns a collection of geometries. This function is the theoretical
+        opposite of the union of the split geometry parts. If the splitter does not split the geometry, a collection
+        with a single geometry equal to the input geometry is returned.
+        The function supports:
+          - Splitting a (Multi)LineString by a (Multi)Point or (Multi)LineString or (Multi)Polygon
+          - Splitting a (Multi)Polygon by a LineString
+
+        It may be convenient to snap the splitter with low tolerance to the geometry. For example in the case
+        of splitting a line by a point, the point must be exactly on the line, for the line to be correctly split.
+        When splitting a line by a polygon, the boundary of the polygon is used for the operation.
+        When splitting a line by another line, a ValueError is raised if the two overlap at some segment.
+
+        Parameters
+        ----------
+        geom : geometry
+            The geometry to be split
+        splitter : geometry
+            The geometry that will split the input geom
+
+        Example
+        -------
+        >>> pt = Point((1, 1))
+        >>> line = LineString([(0,0), (2,2)])
+        >>> result = split(line, pt)
+        >>> result.wkt
+        'GEOMETRYCOLLECTION (LINESTRING (0 0, 1 1), LINESTRING (1 1, 2 2))'
+        """
+
+        if geom.geom_type in ("MultiLineString", "MultiPolygon"):
+            return GeometryCollection(
+                [i for part in geom.geoms for i in SplitOp.split(part, splitter).geoms]
+            )
+
+        elif geom.geom_type == "LineString":
+            if splitter.geom_type in (
+                "LineString",
+                "MultiLineString",
+                "Polygon",
+                "MultiPolygon",
+            ):
+                split_func = SplitOp._split_line_with_line
+            elif splitter.geom_type == "Point":
+                split_func = SplitOp._split_line_with_point
+            elif splitter.geom_type == "MultiPoint":
+                split_func = SplitOp._split_line_with_multipoint
+            else:
+                raise GeometryTypeError(
+                    f"Splitting a LineString with a {splitter.geom_type} is not supported"
+                )
+
+        elif geom.geom_type == "Polygon":
+            if splitter.geom_type == "LineString":
+                split_func = SplitOp._split_polygon_with_line
+            else:
+                raise GeometryTypeError(
+                    f"Splitting a Polygon with a {splitter.geom_type} is not supported"
+                )
+
+        else:
+            raise GeometryTypeError(
+                f"Splitting {geom.geom_type} geometry is not supported"
+            )
+
+        return GeometryCollection(split_func(geom, splitter))
+
+
+split = SplitOp.split
+
+
+def substring(geom, start_dist, end_dist, normalized=False):
+    """Return a line segment between specified distances along a LineString
+
+    Negative distance values are taken as measured in the reverse
+    direction from the end of the geometry. Out-of-range index
+    values are handled by clamping them to the valid range of values.
+
+    If the start distance equals the end distance, a Point is returned.
+
+    If the start distance is actually beyond the end distance, then the
+    reversed substring is returned such that the start distance is
+    at the first coordinate.
+
+    Parameters
+    ----------
+    geom : LineString
+        The geometry to get a substring of.
+    start_dist : float
+        The distance along `geom` of the start of the substring.
+    end_dist : float
+        The distance along `geom` of the end of the substring.
+    normalized : bool, False
+        Whether the distance parameters are interpreted as a
+        fraction of the geometry's length.
+
+    Returns
+    -------
+    Union[Point, LineString]
+        The substring between `start_dist` and `end_dist` or a Point
+        if they are at the same location.
+
+    Raises
+    ------
+    TypeError
+        If `geom` is not a LineString.
+
+    Examples
+    --------
+    >>> from shapely.geometry import LineString
+    >>> from shapely.ops import substring
+    >>> ls = LineString((i, 0) for i in range(6))
+    >>> ls.wkt
+    'LINESTRING (0 0, 1 0, 2 0, 3 0, 4 0, 5 0)'
+    >>> substring(ls, start_dist=1, end_dist=3).wkt
+    'LINESTRING (1 0, 2 0, 3 0)'
+    >>> substring(ls, start_dist=3, end_dist=1).wkt
+    'LINESTRING (3 0, 2 0, 1 0)'
+    >>> substring(ls, start_dist=1, end_dist=-3).wkt
+    'LINESTRING (1 0, 2 0)'
+    >>> substring(ls, start_dist=0.2, end_dist=-0.6, normalized=True).wkt
+    'LINESTRING (1 0, 2 0)'
+
+    Returning a `Point` when `start_dist` and `end_dist` are at the
+    same location.
+
+    >>> substring(ls, 2.5, -2.5).wkt
+    'POINT (2.5 0)'
+    """
+
+    if not isinstance(geom, LineString):
+        raise GeometryTypeError(
+            "Can only calculate a substring of LineString geometries. "
+            f"A {geom.geom_type} was provided."
+        )
+
+    # Filter out cases in which to return a point
+    if start_dist == end_dist:
+        return geom.interpolate(start_dist, normalized)
+    elif not normalized and start_dist >= geom.length and end_dist >= geom.length:
+        return geom.interpolate(geom.length, normalized)
+    elif not normalized and -start_dist >= geom.length and -end_dist >= geom.length:
+        return geom.interpolate(0, normalized)
+    elif normalized and start_dist >= 1 and end_dist >= 1:
+        return geom.interpolate(1, normalized)
+    elif normalized and -start_dist >= 1 and -end_dist >= 1:
+        return geom.interpolate(0, normalized)
+
+    if normalized:
+        start_dist *= geom.length
+        end_dist *= geom.length
+
+    # Filter out cases where distances meet at a middle point from opposite ends.
+    if start_dist < 0 < end_dist and abs(start_dist) + end_dist == geom.length:
+        return geom.interpolate(end_dist)
+    elif end_dist < 0 < start_dist and abs(end_dist) + start_dist == geom.length:
+        return geom.interpolate(start_dist)
+
+    start_point = geom.interpolate(start_dist)
+    end_point = geom.interpolate(end_dist)
+
+    if start_dist < 0:
+        start_dist = geom.length + start_dist  # Values may still be negative,
+    if end_dist < 0:  # but only in the out-of-range
+        end_dist = geom.length + end_dist  # sense, not the wrap-around sense.
+
+    reverse = start_dist > end_dist
+    if reverse:
+        start_dist, end_dist = end_dist, start_dist
+
+    if start_dist < 0:
+        start_dist = 0  # to avoid duplicating the first vertex
+
+    if reverse:
+        vertex_list = [(end_point.x, end_point.y)]
+    else:
+        vertex_list = [(start_point.x, start_point.y)]
+
+    coords = list(geom.coords)
+    current_distance = 0
+    for p1, p2 in zip(coords, coords[1:]):
+        if start_dist < current_distance < end_dist:
+            vertex_list.append(p1)
+        elif current_distance >= end_dist:
+            break
+
+        current_distance += ((p2[0] - p1[0]) ** 2 + (p2[1] - p1[1]) ** 2) ** 0.5
+
+    if reverse:
+        vertex_list.append((start_point.x, start_point.y))
+        # reverse direction result
+        vertex_list = reversed(vertex_list)
+    else:
+        vertex_list.append((end_point.x, end_point.y))
+
+    return LineString(vertex_list)
+
+
+def clip_by_rect(geom, xmin, ymin, xmax, ymax):
+    """Returns the portion of a geometry within a rectangle
+
+    The geometry is clipped in a fast but possibly dirty way. The output is
+    not guaranteed to be valid. No exceptions will be raised for topological
+    errors.
+
+    Parameters
+    ----------
+    geom : geometry
+        The geometry to be clipped
+    xmin : float
+        Minimum x value of the rectangle
+    ymin : float
+        Minimum y value of the rectangle
+    xmax : float
+        Maximum x value of the rectangle
+    ymax : float
+        Maximum y value of the rectangle
+
+    Notes
+    -----
+    Requires GEOS >= 3.5.0
+    New in 1.7.
+    """
+    if geom.is_empty:
+        return geom
+    return shapely.clip_by_rect(geom, xmin, ymin, xmax, ymax)
+
+
+def orient(geom, sign=1.0):
+    """A properly oriented copy of the given geometry.
+
+    The signed area of the result will have the given sign. A sign of
+    1.0 means that the coordinates of the product's exterior rings will
+    be oriented counter-clockwise.
+
+    Parameters
+    ----------
+    geom : Geometry
+        The original geometry. May be a Polygon, MultiPolygon, or
+        GeometryCollection.
+    sign : float, optional.
+        The sign of the result's signed area.
+
+    Returns
+    -------
+    Geometry
+
+    """
+    if isinstance(geom, BaseMultipartGeometry):
+        return geom.__class__(
+            list(
+                map(
+                    lambda geom: orient(geom, sign),
+                    geom.geoms,
+                )
+            )
+        )
+    if isinstance(geom, (Polygon,)):
+        return orient_(geom, sign)
+    return geom

billinglayer/python/shapely/plotting.py

@@ -0,0 +1,217 @@
+"""
+Plot single geometries using Matplotlib.
+
+Note: this module is experimental, and mainly targetting (interactive)
+exploration, debugging and illustration purposes.
+
+"""
+
+import numpy as np
+
+import shapely
+
+
+def _default_ax():
+    import matplotlib.pyplot as plt
+
+    ax = plt.gca()
+    ax.grid(True)
+    ax.set_aspect("equal")
+    return ax
+
+
+def _path_from_polygon(polygon):
+    from matplotlib.path import Path
+
+    if isinstance(polygon, shapely.MultiPolygon):
+        return Path.make_compound_path(
+            *[_path_from_polygon(poly) for poly in polygon.geoms]
+        )
+    else:
+        return Path.make_compound_path(
+            Path(np.asarray(polygon.exterior.coords)[:, :2]),
+            *[Path(np.asarray(ring.coords)[:, :2]) for ring in polygon.interiors],
+        )
+
+
+def patch_from_polygon(polygon, **kwargs):
+    """
+    Gets a Matplotlib patch from a (Multi)Polygon.
+
+    Note: this function is experimental, and mainly targetting (interactive)
+    exploration, debugging and illustration purposes.
+
+    Parameters
+    ----------
+    polygon : shapely.Polygon or shapely.MultiPolygon
+    **kwargs
+        Additional keyword arguments passed to the matplotlib Patch.
+
+    Returns
+    -------
+    Matplotlib artist (PathPatch)
+    """
+    from matplotlib.patches import PathPatch
+
+    return PathPatch(_path_from_polygon(polygon), **kwargs)
+
+
+def plot_polygon(
+    polygon,
+    ax=None,
+    add_points=True,
+    color=None,
+    facecolor=None,
+    edgecolor=None,
+    linewidth=None,
+    **kwargs
+):
+    """
+    Plot a (Multi)Polygon.
+
+    Note: this function is experimental, and mainly targetting (interactive)
+    exploration, debugging and illustration purposes.
+
+    Parameters
+    ----------
+    polygon : shapely.Polygon or shapely.MultiPolygon
+    ax : matplotlib Axes, default None
+        The axes on which to draw the plot. If not specified, will get the
+        current active axes or create a new figure.
+    add_points : bool, default True
+        If True, also plot the coordinates (vertices) as points.
+    color : matplotlib color specification
+        Color for both the polygon fill (face) and boundary (edge). By default,
+        the fill is using an alpha of 0.3. You can specify `facecolor` and
+        `edgecolor` separately for greater control.
+    facecolor : matplotlib color specification
+        Color for the polygon fill.
+    edgecolor : matplotlib color specification
+        Color for the polygon boundary.
+    linewidth : float
+        The line width for the polygon boundary.
+    **kwargs
+        Additional keyword arguments passed to the matplotlib Patch.
+
+    Returns
+    -------
+    Matplotlib artist (PathPatch), if `add_points` is false.
+    A tuple of Matplotlib artists (PathPatch, Line2D), if `add_points` is true.
+    """
+    from matplotlib import colors
+
+    if ax is None:
+        ax = _default_ax()
+
+    if color is None:
+        color = "C0"
+    color = colors.to_rgba(color)
+
+    if facecolor is None:
+        facecolor = list(color)
+        facecolor[-1] = 0.3
+        facecolor = tuple(facecolor)
+
+    if edgecolor is None:
+        edgecolor = color
+
+    patch = patch_from_polygon(
+        polygon, facecolor=facecolor, edgecolor=edgecolor, linewidth=linewidth, **kwargs
+    )
+    ax.add_patch(patch)
+    ax.autoscale_view()
+
+    if add_points:
+        line = plot_points(polygon, ax=ax, color=color)
+        return patch, line
+
+    return patch
+
+
+def plot_line(line, ax=None, add_points=True, color=None, linewidth=2, **kwargs):
+    """
+    Plot a (Multi)LineString/LinearRing.
+
+    Note: this function is experimental, and mainly targetting (interactive)
+    exploration, debugging and illustration purposes.
+
+    Parameters
+    ----------
+    line : shapely.LineString or shapely.LinearRing
+    ax : matplotlib Axes, default None
+        The axes on which to draw the plot. If not specified, will get the
+        current active axes or create a new figure.
+    add_points : bool, default True
+        If True, also plot the coordinates (vertices) as points.
+    color : matplotlib color specification
+        Color for the line (edgecolor under the hood) and pointes.
+    linewidth : float, default 2
+        The line width for the polygon boundary.
+    **kwargs
+        Additional keyword arguments passed to the matplotlib Patch.
+
+    Returns
+    -------
+    Matplotlib artist (PathPatch)
+    """
+    from matplotlib.patches import PathPatch
+    from matplotlib.path import Path
+
+    if ax is None:
+        ax = _default_ax()
+
+    if color is None:
+        color = "C0"
+
+    if isinstance(line, shapely.MultiLineString):
+        path = Path.make_compound_path(
+            *[Path(np.asarray(mline.coords)[:, :2]) for mline in line.geoms]
+        )
+    else:
+        path = Path(np.asarray(line.coords)[:, :2])
+
+    patch = PathPatch(
+        path, facecolor="none", edgecolor=color, linewidth=linewidth, **kwargs
+    )
+    ax.add_patch(patch)
+    ax.autoscale_view()
+
+    if add_points:
+        line = plot_points(line, ax=ax, color=color)
+        return patch, line
+
+    return patch
+
+
+def plot_points(geom, ax=None, color=None, marker="o", **kwargs):
+    """
+    Plot a Point/MultiPoint or the vertices of any other geometry type.
+
+    Parameters
+    ----------
+    geom : shapely.Geometry
+        Any shapely Geometry object, from which all vertices are extracted
+        and plotted.
+    ax : matplotlib Axes, default None
+        The axes on which to draw the plot. If not specified, will get the
+        current active axes or create a new figure.
+    color : matplotlib color specification
+        Color for the filled points. You can use `markeredgecolor` and
+        `markeredgecolor` to have different edge and fill colors.
+    marker : str, default "o"
+        The matplotlib marker for the points.
+    **kwargs
+        Additional keyword arguments passed to matplotlib `plot` (Line2D).
+
+    Returns
+    -------
+    Matplotlib artist (Line2D)
+    """
+    if ax is None:
+        ax = _default_ax()
+
+    coords = shapely.get_coordinates(geom)
+    (line,) = ax.plot(
+        coords[:, 0], coords[:, 1], linestyle="", marker=marker, color=color, **kwargs
+    )
+    return line

billinglayer/python/shapely/prepared.py

@@ -0,0 +1,77 @@
+"""
+Support for GEOS prepared geometry operations.
+"""
+from pickle import PicklingError
+
+import shapely
+
+
+class PreparedGeometry:
+    """
+    A geometry prepared for efficient comparison to a set of other geometries.
+
+    Example:
+
+      >>> from shapely.geometry import Point, Polygon
+      >>> triangle = Polygon([(0.0, 0.0), (1.0, 1.0), (1.0, -1.0)])
+      >>> p = prep(triangle)
+      >>> p.intersects(Point(0.5, 0.5))
+      True
+    """
+
+    def __init__(self, context):
+        if isinstance(context, PreparedGeometry):
+            self.context = context.context
+        else:
+            shapely.prepare(context)
+            self.context = context
+        self.prepared = True
+
+    def contains(self, other):
+        """Returns True if the geometry contains the other, else False"""
+        return self.context.contains(other)
+
+    def contains_properly(self, other):
+        """Returns True if the geometry properly contains the other, else False"""
+        # TODO temporary hack until shapely exposes contains properly as predicate function
+        from shapely import STRtree
+
+        tree = STRtree([other])
+        idx = tree.query(self.context, predicate="contains_properly")
+        return bool(len(idx))
+
+    def covers(self, other):
+        """Returns True if the geometry covers the other, else False"""
+        return self.context.covers(other)
+
+    def crosses(self, other):
+        """Returns True if the geometries cross, else False"""
+        return self.context.crosses(other)
+
+    def disjoint(self, other):
+        """Returns True if geometries are disjoint, else False"""
+        return self.context.disjoint(other)
+
+    def intersects(self, other):
+        """Returns True if geometries intersect, else False"""
+        return self.context.intersects(other)
+
+    def overlaps(self, other):
+        """Returns True if geometries overlap, else False"""
+        return self.context.overlaps(other)
+
+    def touches(self, other):
+        """Returns True if geometries touch, else False"""
+        return self.context.touches(other)
+
+    def within(self, other):
+        """Returns True if geometry is within the other, else False"""
+        return self.context.within(other)
+
+    def __reduce__(self):
+        raise PicklingError("Prepared geometries cannot be pickled.")
+
+
+def prep(ob):
+    """Creates and returns a prepared geometric object."""
+    return PreparedGeometry(ob)

billinglayer/python/shapely/set_operations.py

@@ -0,0 +1,514 @@
+import numpy as np
+
+from . import GeometryType, lib
+from .decorators import multithreading_enabled, requires_geos
+from .errors import UnsupportedGEOSVersionError
+
+__all__ = [
+    "difference",
+    "intersection",
+    "intersection_all",
+    "symmetric_difference",
+    "symmetric_difference_all",
+    "unary_union",
+    "union",
+    "union_all",
+    "coverage_union",
+    "coverage_union_all",
+]
+
+
+@multithreading_enabled
+def difference(a, b, grid_size=None, **kwargs):
+    """Returns the part of geometry A that does not intersect with geometry B.
+
+    If grid_size is nonzero, input coordinates will be snapped to a precision
+    grid of that size and resulting coordinates will be snapped to that same
+    grid.  If 0, this operation will use double precision coordinates.  If None,
+    the highest precision of the inputs will be used, which may be previously
+    set using set_precision.  Note: returned geometry does not have precision
+    set unless specified previously by set_precision.
+
+    Parameters
+    ----------
+    a : Geometry or array_like
+    b : Geometry or array_like
+    grid_size : float, optional
+        Precision grid size; requires GEOS >= 3.9.0.  Will use the highest
+        precision of the inputs by default.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    set_precision
+
+    Examples
+    --------
+    >>> from shapely import box, LineString, normalize, Polygon
+    >>> line = LineString([(0, 0), (2, 2)])
+    >>> difference(line, LineString([(1, 1), (3, 3)]))
+    <LINESTRING (0 0, 1 1)>
+    >>> difference(line, LineString())
+    <LINESTRING (0 0, 2 2)>
+    >>> difference(line, None) is None
+    True
+    >>> box1 = box(0, 0, 2, 2)
+    >>> box2 = box(1, 1, 3, 3)
+    >>> normalize(difference(box1, box2))
+    <POLYGON ((0 0, 0 2, 1 2, 1 1, 2 1, 2 0, 0 0))>
+    >>> box1 = box(0.1, 0.2, 2.1, 2.1)
+    >>> difference(box1, box2, grid_size=1)
+    <POLYGON ((2 0, 0 0, 0 2, 1 2, 1 1, 2 1, 2 0))>
+    """
+
+    if grid_size is not None:
+        if lib.geos_version < (3, 9, 0):
+            raise UnsupportedGEOSVersionError(
+                "grid_size parameter requires GEOS >= 3.9.0"
+            )
+
+        if not np.isscalar(grid_size):
+            raise ValueError("grid_size parameter only accepts scalar values")
+
+        return lib.difference_prec(a, b, grid_size, **kwargs)
+
+    return lib.difference(a, b, **kwargs)
+
+
+@multithreading_enabled
+def intersection(a, b, grid_size=None, **kwargs):
+    """Returns the geometry that is shared between input geometries.
+
+    If grid_size is nonzero, input coordinates will be snapped to a precision
+    grid of that size and resulting coordinates will be snapped to that same
+    grid.  If 0, this operation will use double precision coordinates.  If None,
+    the highest precision of the inputs will be used, which may be previously
+    set using set_precision.  Note: returned geometry does not have precision
+    set unless specified previously by set_precision.
+
+    Parameters
+    ----------
+    a : Geometry or array_like
+    b : Geometry or array_like
+    grid_size : float, optional
+        Precision grid size; requires GEOS >= 3.9.0.  Will use the highest
+        precision of the inputs by default.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    intersection_all
+    set_precision
+
+    Examples
+    --------
+    >>> from shapely import box, LineString, normalize, Polygon
+    >>> line = LineString([(0, 0), (2, 2)])
+    >>> intersection(line, LineString([(1, 1), (3, 3)]))
+    <LINESTRING (1 1, 2 2)>
+    >>> box1 = box(0, 0, 2, 2)
+    >>> box2 = box(1, 1, 3, 3)
+    >>> normalize(intersection(box1, box2))
+    <POLYGON ((1 1, 1 2, 2 2, 2 1, 1 1))>
+    >>> box1 = box(0.1, 0.2, 2.1, 2.1)
+    >>> intersection(box1, box2, grid_size=1)
+    <POLYGON ((2 2, 2 1, 1 1, 1 2, 2 2))>
+    """
+
+    if grid_size is not None:
+        if lib.geos_version < (3, 9, 0):
+            raise UnsupportedGEOSVersionError(
+                "grid_size parameter requires GEOS >= 3.9.0"
+            )
+
+        if not np.isscalar(grid_size):
+            raise ValueError("grid_size parameter only accepts scalar values")
+
+        return lib.intersection_prec(a, b, grid_size, **kwargs)
+
+    return lib.intersection(a, b, **kwargs)
+
+
+@multithreading_enabled
+def intersection_all(geometries, axis=None, **kwargs):
+    """Returns the intersection of multiple geometries.
+
+    This function ignores None values when other Geometry elements are present.
+    If all elements of the given axis are None, an empty GeometryCollection is
+    returned.
+
+    Parameters
+    ----------
+    geometries : array_like
+    axis : int, optional
+        Axis along which the operation is performed. The default (None)
+        performs the operation over all axes, returning a scalar value.
+        Axis may be negative, in which case it counts from the last to the
+        first axis.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc.reduce docs <https://numpy.org/doc/stable/reference/generated/numpy.ufunc.reduce.html#numpy.ufunc.reduce>`_.
+
+    See also
+    --------
+    intersection
+
+    Examples
+    --------
+    >>> from shapely import LineString
+    >>> line1 = LineString([(0, 0), (2, 2)])
+    >>> line2 = LineString([(1, 1), (3, 3)])
+    >>> intersection_all([line1, line2])
+    <LINESTRING (1 1, 2 2)>
+    >>> intersection_all([[line1, line2, None]], axis=1).tolist()
+    [<LINESTRING (1 1, 2 2)>]
+    >>> intersection_all([line1, None])
+    <LINESTRING (0 0, 2 2)>
+    """
+    geometries = np.asarray(geometries)
+    if axis is None:
+        geometries = geometries.ravel()
+    else:
+        geometries = np.rollaxis(geometries, axis=axis, start=geometries.ndim)
+
+    return lib.intersection_all(geometries, **kwargs)
+
+
+@multithreading_enabled
+def symmetric_difference(a, b, grid_size=None, **kwargs):
+    """Returns the geometry that represents the portions of input geometries
+    that do not intersect.
+
+    If grid_size is nonzero, input coordinates will be snapped to a precision
+    grid of that size and resulting coordinates will be snapped to that same
+    grid.  If 0, this operation will use double precision coordinates.  If None,
+    the highest precision of the inputs will be used, which may be previously
+    set using set_precision.  Note: returned geometry does not have precision
+    set unless specified previously by set_precision.
+
+    Parameters
+    ----------
+    a : Geometry or array_like
+    b : Geometry or array_like
+    grid_size : float, optional
+        Precision grid size; requires GEOS >= 3.9.0.  Will use the highest
+        precision of the inputs by default.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    symmetric_difference_all
+    set_precision
+
+    Examples
+    --------
+    >>> from shapely import box, LineString, normalize
+    >>> line = LineString([(0, 0), (2, 2)])
+    >>> symmetric_difference(line, LineString([(1, 1), (3, 3)]))
+    <MULTILINESTRING ((0 0, 1 1), (2 2, 3 3))>
+    >>> box1 = box(0, 0, 2, 2)
+    >>> box2 = box(1, 1, 3, 3)
+    >>> normalize(symmetric_difference(box1, box2))
+    <MULTIPOLYGON (((1 2, 1 3, 3 3, 3 1, 2 1, 2 2, 1 2)), ((0 0, 0 2, 1 2, 1 1, ...>
+    >>> box1 = box(0.1, 0.2, 2.1, 2.1)
+    >>> symmetric_difference(box1, box2, grid_size=1)
+    <MULTIPOLYGON (((2 0, 0 0, 0 2, 1 2, 1 1, 2 1, 2 0)), ((2 2, 1 2, 1 3, 3 3, ...>
+    """
+
+    if grid_size is not None:
+        if lib.geos_version < (3, 9, 0):
+            raise UnsupportedGEOSVersionError(
+                "grid_size parameter requires GEOS >= 3.9.0"
+            )
+
+        if not np.isscalar(grid_size):
+            raise ValueError("grid_size parameter only accepts scalar values")
+
+        return lib.symmetric_difference_prec(a, b, grid_size, **kwargs)
+
+    return lib.symmetric_difference(a, b, **kwargs)
+
+
+@multithreading_enabled
+def symmetric_difference_all(geometries, axis=None, **kwargs):
+    """Returns the symmetric difference of multiple geometries.
+
+    This function ignores None values when other Geometry elements are present.
+    If all elements of the given axis are None an empty GeometryCollection is
+    returned.
+
+    Parameters
+    ----------
+    geometries : array_like
+    axis : int, optional
+        Axis along which the operation is performed. The default (None)
+        performs the operation over all axes, returning a scalar value.
+        Axis may be negative, in which case it counts from the last to the
+        first axis.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc.reduce docs <https://numpy.org/doc/stable/reference/generated/numpy.ufunc.reduce.html#numpy.ufunc.reduce>`_.
+
+    See also
+    --------
+    symmetric_difference
+
+    Examples
+    --------
+    >>> from shapely import LineString
+    >>> line1 = LineString([(0, 0), (2, 2)])
+    >>> line2 = LineString([(1, 1), (3, 3)])
+    >>> symmetric_difference_all([line1, line2])
+    <MULTILINESTRING ((0 0, 1 1), (2 2, 3 3))>
+    >>> symmetric_difference_all([[line1, line2, None]], axis=1).tolist()
+    [<MULTILINESTRING ((0 0, 1 1), (2 2, 3 3))>]
+    >>> symmetric_difference_all([line1, None])
+    <LINESTRING (0 0, 2 2)>
+    >>> symmetric_difference_all([None, None])
+    <GEOMETRYCOLLECTION EMPTY>
+    """
+    geometries = np.asarray(geometries)
+    if axis is None:
+        geometries = geometries.ravel()
+    else:
+        geometries = np.rollaxis(geometries, axis=axis, start=geometries.ndim)
+
+    return lib.symmetric_difference_all(geometries, **kwargs)
+
+
+@multithreading_enabled
+def union(a, b, grid_size=None, **kwargs):
+    """Merges geometries into one.
+
+    If grid_size is nonzero, input coordinates will be snapped to a precision
+    grid of that size and resulting coordinates will be snapped to that same
+    grid.  If 0, this operation will use double precision coordinates.  If None,
+    the highest precision of the inputs will be used, which may be previously
+    set using set_precision.  Note: returned geometry does not have precision
+    set unless specified previously by set_precision.
+
+    Parameters
+    ----------
+    a : Geometry or array_like
+    b : Geometry or array_like
+    grid_size : float, optional
+        Precision grid size; requires GEOS >= 3.9.0.  Will use the highest
+        precision of the inputs by default.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    union_all
+    set_precision
+
+    Examples
+    --------
+    >>> from shapely import box, LineString, normalize
+    >>> line = LineString([(0, 0), (2, 2)])
+    >>> union(line, LineString([(2, 2), (3, 3)]))
+    <MULTILINESTRING ((0 0, 2 2), (2 2, 3 3))>
+    >>> union(line, None) is None
+    True
+    >>> box1 = box(0, 0, 2, 2)
+    >>> box2 = box(1, 1, 3, 3)
+    >>> normalize(union(box1, box2))
+    <POLYGON ((0 0, 0 2, 1 2, 1 3, 3 3, 3 1, 2 1, 2 0, 0 0))>
+    >>> box1 = box(0.1, 0.2, 2.1, 2.1)
+    >>> union(box1, box2, grid_size=1)
+    <POLYGON ((2 0, 0 0, 0 2, 1 2, 1 3, 3 3, 3 1, 2 1, 2 0))>
+    """
+
+    if grid_size is not None:
+        if lib.geos_version < (3, 9, 0):
+            raise UnsupportedGEOSVersionError(
+                "grid_size parameter requires GEOS >= 3.9.0"
+            )
+
+        if not np.isscalar(grid_size):
+            raise ValueError("grid_size parameter only accepts scalar values")
+
+        return lib.union_prec(a, b, grid_size, **kwargs)
+
+    return lib.union(a, b, **kwargs)
+
+
+@multithreading_enabled
+def union_all(geometries, grid_size=None, axis=None, **kwargs):
+    """Returns the union of multiple geometries.
+
+    This function ignores None values when other Geometry elements are present.
+    If all elements of the given axis are None an empty GeometryCollection is
+    returned.
+
+    If grid_size is nonzero, input coordinates will be snapped to a precision
+    grid of that size and resulting coordinates will be snapped to that same
+    grid.  If 0, this operation will use double precision coordinates.  If None,
+    the highest precision of the inputs will be used, which may be previously
+    set using set_precision.  Note: returned geometry does not have precision
+    set unless specified previously by set_precision.
+
+    `unary_union` is an alias of `union_all`.
+
+    Parameters
+    ----------
+    geometries : array_like
+    grid_size : float, optional
+        Precision grid size; requires GEOS >= 3.9.0.  Will use the highest
+        precision of the inputs by default.
+    axis : int, optional
+        Axis along which the operation is performed. The default (None)
+        performs the operation over all axes, returning a scalar value.
+        Axis may be negative, in which case it counts from the last to the
+        first axis.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    union
+    set_precision
+
+    Examples
+    --------
+    >>> from shapely import box, LineString, normalize, Point
+    >>> line1 = LineString([(0, 0), (2, 2)])
+    >>> line2 = LineString([(2, 2), (3, 3)])
+    >>> union_all([line1, line2])
+    <MULTILINESTRING ((0 0, 2 2), (2 2, 3 3))>
+    >>> union_all([[line1, line2, None]], axis=1).tolist()
+    [<MULTILINESTRING ((0 0, 2 2), (2 2, 3 3))>]
+    >>> box1 = box(0, 0, 2, 2)
+    >>> box2 = box(1, 1, 3, 3)
+    >>> normalize(union_all([box1, box2]))
+    <POLYGON ((0 0, 0 2, 1 2, 1 3, 3 3, 3 1, 2 1, 2 0, 0 0))>
+    >>> box1 = box(0.1, 0.2, 2.1, 2.1)
+    >>> union_all([box1, box2], grid_size=1)
+    <POLYGON ((2 0, 0 0, 0 2, 1 2, 1 3, 3 3, 3 1, 2 1, 2 0))>
+    >>> union_all([None, Point(0, 1)])
+    <POINT (0 1)>
+    >>> union_all([None, None])
+    <GEOMETRYCOLLECTION EMPTY>
+    >>> union_all([])
+    <GEOMETRYCOLLECTION EMPTY>
+    """
+    # for union_all, GEOS provides an efficient route through first creating
+    # GeometryCollections
+    # first roll the aggregation axis backwards
+    geometries = np.asarray(geometries)
+    if axis is None:
+        geometries = geometries.ravel()
+    else:
+        geometries = np.rollaxis(geometries, axis=axis, start=geometries.ndim)
+
+    # create_collection acts on the inner axis
+    collections = lib.create_collection(geometries, GeometryType.GEOMETRYCOLLECTION)
+
+    if grid_size is not None:
+        if lib.geos_version < (3, 9, 0):
+            raise UnsupportedGEOSVersionError(
+                "grid_size parameter requires GEOS >= 3.9.0"
+            )
+
+        if not np.isscalar(grid_size):
+            raise ValueError("grid_size parameter only accepts scalar values")
+
+        return lib.unary_union_prec(collections, grid_size, **kwargs)
+
+    return lib.unary_union(collections, **kwargs)
+
+
+unary_union = union_all
+
+
+@requires_geos("3.8.0")
+@multithreading_enabled
+def coverage_union(a, b, **kwargs):
+    """Merges multiple polygons into one. This is an optimized version of
+    union which assumes the polygons to be non-overlapping.
+
+    Parameters
+    ----------
+    a : Geometry or array_like
+    b : Geometry or array_like
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    coverage_union_all
+
+    Examples
+    --------
+    >>> from shapely import normalize, Polygon
+    >>> polygon = Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
+    >>> normalize(coverage_union(polygon, Polygon([(1, 0), (1, 1), (2, 1), (2, 0), (1, 0)])))
+    <POLYGON ((0 0, 0 1, 1 1, 2 1, 2 0, 1 0, 0 0))>
+
+    Union with None returns same polygon
+    >>> normalize(coverage_union(polygon, None))
+    <POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))>
+    """
+    return coverage_union_all([a, b], **kwargs)
+
+
+@requires_geos("3.8.0")
+@multithreading_enabled
+def coverage_union_all(geometries, axis=None, **kwargs):
+    """Returns the union of multiple polygons of a geometry collection.
+    This is an optimized version of union which assumes the polygons
+    to be non-overlapping.
+
+    This function ignores None values when other Geometry elements are present.
+    If all elements of the given axis are None, an empty MultiPolygon is
+    returned.
+
+    Parameters
+    ----------
+    geometries : array_like
+    axis : int, optional
+        Axis along which the operation is performed. The default (None)
+        performs the operation over all axes, returning a scalar value.
+        Axis may be negative, in which case it counts from the last to the
+        first axis.
+    **kwargs
+        For other keyword-only arguments, see the
+        `NumPy ufunc docs <https://numpy.org/doc/stable/reference/ufuncs.html#ufuncs-kwargs>`_.
+
+    See also
+    --------
+    coverage_union
+
+    Examples
+    --------
+    >>> from shapely import normalize, Polygon
+    >>> polygon_1 = Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
+    >>> polygon_2 = Polygon([(1, 0), (1, 1), (2, 1), (2, 0), (1, 0)])
+    >>> normalize(coverage_union_all([polygon_1, polygon_2]))
+    <POLYGON ((0 0, 0 1, 1 1, 2 1, 2 0, 1 0, 0 0))>
+    >>> normalize(coverage_union_all([polygon_1, None]))
+    <POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))>
+    >>> normalize(coverage_union_all([None, None]))
+    <MULTIPOLYGON EMPTY>
+    """
+    # coverage union in GEOS works over GeometryCollections
+    # first roll the aggregation axis backwards
+    geometries = np.asarray(geometries)
+    if axis is None:
+        geometries = geometries.ravel()
+    else:
+        geometries = np.rollaxis(
+            np.asarray(geometries), axis=axis, start=geometries.ndim
+        )
+    # create_collection acts on the inner axis
+    collections = lib.create_collection(geometries, GeometryType.GEOMETRYCOLLECTION)
+    return lib.coverage_union(collections, **kwargs)

billinglayer/python/shapely/speedups.py

@@ -0,0 +1,36 @@
+import warnings
+
+__all__ = ["available", "enable", "disable", "enabled"]
+
+
+available = True
+enabled = True
+
+
+_MSG = (
+    "This function has no longer any effect, and will be removed in a "
+    "future release. Starting with Shapely 2.0, equivalent speedups are "
+    "always available"
+)
+
+
+def enable():
+    """
+    This function has no longer any effect, and will be removed in a future
+    release.
+
+    Previously, this function enabled cython-based speedups. Starting with
+    Shapely 2.0, equivalent speedups are available in every installation.
+    """
+    warnings.warn(_MSG, DeprecationWarning, stacklevel=2)
+
+
+def disable():
+    """
+    This function has no longer any effect, and will be removed in a future
+    release.
+
+    Previously, this function enabled cython-based speedups. Starting with
+    Shapely 2.0, equivalent speedups are available in every installation.
+    """
+    warnings.warn(_MSG, DeprecationWarning, stacklevel=2)

billinglayer/python/shapely/strtree.py

@@ -0,0 +1,544 @@
+from typing import Any, Iterable, Union
+
+import numpy as np
+
+from . import lib
+from ._enum import ParamEnum
+from .decorators import requires_geos, UnsupportedGEOSVersionError
+from .geometry.base import BaseGeometry
+from .predicates import is_empty, is_missing
+
+__all__ = ["STRtree"]
+
+
+class BinaryPredicate(ParamEnum):
+    """The enumeration of GEOS binary predicates types"""
+
+    intersects = 1
+    within = 2
+    contains = 3
+    overlaps = 4
+    crosses = 5
+    touches = 6
+    covers = 7
+    covered_by = 8
+    contains_properly = 9
+
+
+class STRtree:
+    """
+    A query-only R-tree spatial index created using the
+    Sort-Tile-Recursive (STR) [1]_ algorithm.
+
+    The tree indexes the bounding boxes of each geometry.  The tree is
+    constructed directly at initialization and nodes cannot be added or
+    removed after it has been created.
+
+    All operations return indices of the input geometries.  These indices
+    can be used to index into anything associated with the input geometries,
+    including the input geometries themselves, or custom items stored in
+    another object of the same length as the geometries.
+
+    Bounding boxes limited to two dimensions and are axis-aligned (equivalent to
+    the ``bounds`` property of a geometry); any Z values present in geometries
+    are ignored for purposes of indexing within the tree.
+
+    Any mixture of geometry types may be stored in the tree.
+
+    Note: the tree is more efficient for querying when there are fewer
+    geometries that have overlapping bounding boxes and where there is greater
+    similarity between the outer boundary of a geometry and its bounding box.
+    For example, a MultiPolygon composed of widely-spaced individual Polygons
+    will have a large overall bounding box compared to the boundaries of its
+    individual Polygons, and the bounding box may also potentially overlap many
+    other geometries within the tree.  This means that the resulting tree may be
+    less efficient to query than a tree constructed from individual Polygons.
+
+    Parameters
+    ----------
+    geoms : sequence
+        A sequence of geometry objects.
+    node_capacity : int, default 10
+        The maximum number of child nodes per parent node in the tree.
+
+    References
+    ----------
+    .. [1] Leutenegger, Scott T.; Edgington, Jeffrey M.; Lopez, Mario A.
+       (February 1997). "STR: A Simple and Efficient Algorithm for
+       R-Tree Packing".
+       https://ia600900.us.archive.org/27/items/nasa_techdoc_19970016975/19970016975.pdf
+    """
+
+    def __init__(
+        self,
+        geoms: Iterable[BaseGeometry],
+        node_capacity: int = 10,
+    ):
+        # Keep references to geoms in a copied array so that this array is not
+        # modified while the tree depends on it remaining the same
+        self._geometries = np.array(geoms, dtype=np.object_, copy=True)
+
+        # initialize GEOS STRtree
+        self._tree = lib.STRtree(self.geometries, node_capacity)
+
+    def __len__(self):
+        return self._tree.count
+
+    def __reduce__(self):
+        return (STRtree, (self.geometries,))
+
+    @property
+    def geometries(self):
+        """
+        Geometries stored in the tree in the order used to construct the tree.
+
+        The order of this array corresponds to the tree indices returned by
+        other STRtree methods.
+
+        Do not attempt to modify items in the returned array.
+
+        Returns
+        -------
+        ndarray of Geometry objects
+        """
+        return self._geometries
+
+    def query(self, geometry, predicate=None, distance=None):
+        """
+        Return the integer indices of all combinations of each input geometry
+        and tree geometries where the bounding box of each input geometry
+        intersects the bounding box of a tree geometry.
+
+        If the input geometry is a scalar, this returns an array of shape (n, ) with
+        the indices of the matching tree geometries.  If the input geometry is an
+        array_like, this returns an array with shape (2,n) where the subarrays
+        correspond to the indices of the input geometries and indices of the
+        tree geometries associated with each.  To generate an array of pairs of
+        input geometry index and tree geometry index, simply transpose the
+        result.
+
+        If a predicate is provided, the tree geometries are first queried based
+        on the bounding box of the input geometry and then are further filtered
+        to those that meet the predicate when comparing the input geometry to
+        the tree geometry:
+        predicate(geometry, tree_geometry)
+
+        The 'dwithin' predicate requires GEOS >= 3.10.
+
+        Bounding boxes are limited to two dimensions and are axis-aligned
+        (equivalent to the ``bounds`` property of a geometry); any Z values
+        present in input geometries are ignored when querying the tree.
+
+        Any input geometry that is None or empty will never match geometries in
+        the tree.
+
+        Parameters
+        ----------
+        geometry : Geometry or array_like
+            Input geometries to query the tree and filter results using the
+            optional predicate.
+        predicate : {None, 'intersects', 'within', 'contains', 'overlaps', 'crosses',\
+'touches', 'covers', 'covered_by', 'contains_properly', 'dwithin'}, optional
+            The predicate to use for testing geometries from the tree
+            that are within the input geometry's bounding box.
+        distance : number or array_like, optional
+            Distances around each input geometry within which to query the tree
+            for the 'dwithin' predicate.  If array_like, shape must be
+            broadcastable to shape of geometry.  Required if predicate='dwithin'.
+
+        Returns
+        -------
+        ndarray with shape (n,) if geometry is a scalar
+            Contains tree geometry indices.
+
+        OR
+
+        ndarray with shape (2, n) if geometry is an array_like
+            The first subarray contains input geometry indices.
+            The second subarray contains tree geometry indices.
+
+        Examples
+        --------
+        >>> from shapely import box, Point
+        >>> import numpy as np
+        >>> points = [Point(0, 0), Point(1, 1), Point(2,2), Point(3, 3)]
+        >>> tree = STRtree(points)
+
+        Query the tree using a scalar geometry:
+
+        >>> indices = tree.query(box(0, 0, 1, 1))
+        >>> indices.tolist()
+        [0, 1]
+
+        Query using an array of geometries:
+
+        >>> boxes = np.array([box(0, 0, 1, 1), box(2, 2, 3, 3)])
+        >>> arr_indices = tree.query(boxes)
+        >>> arr_indices.tolist()
+        [[0, 0, 1, 1], [0, 1, 2, 3]]
+
+        Or transpose to get all pairs of input and tree indices:
+
+        >>> arr_indices.T.tolist()
+        [[0, 0], [0, 1], [1, 2], [1, 3]]
+
+        Retrieve the tree geometries by results of query:
+
+        >>> tree.geometries.take(indices).tolist()
+        [<POINT (0 0)>, <POINT (1 1)>]
+
+        Retrieve all pairs of input and tree geometries:
+
+        >>> np.array([boxes.take(arr_indices[0]),\
+tree.geometries.take(arr_indices[1])]).T.tolist()
+        [[<POLYGON ((1 0, 1 1, 0 1, 0 0, 1 0))>, <POINT (0 0)>],
+         [<POLYGON ((1 0, 1 1, 0 1, 0 0, 1 0))>, <POINT (1 1)>],
+         [<POLYGON ((3 2, 3 3, 2 3, 2 2, 3 2))>, <POINT (2 2)>],
+         [<POLYGON ((3 2, 3 3, 2 3, 2 2, 3 2))>, <POINT (3 3)>]]
+
+        Query using a predicate:
+
+        >>> tree = STRtree([box(0, 0, 0.5, 0.5), box(0.5, 0.5, 1, 1), box(1, 1, 2, 2)])
+        >>> tree.query(box(0, 0, 1, 1), predicate="contains").tolist()
+        [0, 1]
+        >>> tree.query(Point(0.75, 0.75), predicate="dwithin", distance=0.5).tolist()
+        [0, 1, 2]
+
+        >>> tree.query(boxes, predicate="contains").tolist()
+        [[0, 0], [0, 1]]
+        >>> tree.query(boxes, predicate="dwithin", distance=0.5).tolist()
+        [[0, 0, 0, 1], [0, 1, 2, 2]]
+
+        Retrieve custom items associated with tree geometries (records can
+        be in whatever data structure so long as geometries and custom data
+        can be extracted into arrays of the same length and order):
+
+        >>> records = [
+        ...     {"geometry": Point(0, 0), "value": "A"},
+        ...     {"geometry": Point(2, 2), "value": "B"}
+        ... ]
+        >>> tree = STRtree([record["geometry"] for record in records])
+        >>> items = np.array([record["value"] for record in records])
+        >>> items.take(tree.query(box(0, 0, 1, 1))).tolist()
+        ['A']
+
+
+        Notes
+        -----
+        In the context of a spatial join, input geometries are the "left"
+        geometries that determine the order of the results, and tree geometries
+        are "right" geometries that are joined against the left geometries. This
+        effectively performs an inner join, where only those combinations of
+        geometries that can be joined based on overlapping bounding boxes or
+        optional predicate are returned.
+        """
+
+        geometry = np.asarray(geometry)
+        is_scalar = False
+        if geometry.ndim == 0:
+            geometry = np.expand_dims(geometry, 0)
+            is_scalar = True
+
+        if predicate is None:
+            indices = self._tree.query(geometry, 0)
+            return indices[1] if is_scalar else indices
+
+        # Requires GEOS >= 3.10
+        elif predicate == "dwithin":
+            if lib.geos_version < (3, 10, 0):
+                raise UnsupportedGEOSVersionError(
+                    "dwithin predicate requires GEOS >= 3.10"
+                )
+            if distance is None:
+                raise ValueError(
+                    "distance parameter must be provided for dwithin predicate"
+                )
+            distance = np.asarray(distance, dtype="float64")
+            if distance.ndim > 1:
+                raise ValueError("Distance array should be one dimensional")
+
+            try:
+                distance = np.broadcast_to(distance, geometry.shape)
+            except ValueError:
+                raise ValueError("Could not broadcast distance to match geometry")
+
+            indices = self._tree.dwithin(geometry, distance)
+            return indices[1] if is_scalar else indices
+
+        predicate = BinaryPredicate.get_value(predicate)
+        indices = self._tree.query(geometry, predicate)
+        return indices[1] if is_scalar else indices
+
+    @requires_geos("3.6.0")
+    def nearest(self, geometry) -> Union[Any, None]:
+        """
+        Return the index of the nearest geometry in the tree for each input
+        geometry based on distance within two-dimensional Cartesian space.
+
+        This distance will be 0 when input geometries intersect tree geometries.
+
+        If there are multiple equidistant or intersected geometries in the tree,
+        only a single result is returned for each input geometry, based on the
+        order that tree geometries are visited; this order may be
+        nondeterministic.
+
+        If any input geometry is None or empty, an error is raised.  Any Z
+        values present in input geometries are ignored when finding nearest
+        tree geometries.
+
+        Parameters
+        ----------
+        geometry : Geometry or array_like
+            Input geometries to query the tree.
+
+        Returns
+        -------
+        scalar or ndarray
+            Indices of geometries in tree. Return value will have the same shape
+            as the input.
+
+            None is returned if this index is empty. This may change in
+            version 2.0.
+
+        See also
+        --------
+        query_nearest: returns all equidistant geometries, exclusive geometries, \
+and optional distances
+
+        Examples
+        --------
+        >>> from shapely.geometry import Point
+        >>> tree = STRtree([Point(i, i) for i in range(10)])
+
+        Query the tree for nearest using a scalar geometry:
+
+        >>> index = tree.nearest(Point(2.2, 2.2))
+        >>> index
+        2
+        >>> tree.geometries.take(index)
+        <POINT (2 2)>
+
+        Query the tree for nearest using an array of geometries:
+
+        >>> indices = tree.nearest([Point(2.2, 2.2), Point(4.4, 4.4)])
+        >>> indices.tolist()
+        [2, 4]
+        >>> tree.geometries.take(indices).tolist()
+        [<POINT (2 2)>, <POINT (4 4)>]
+
+        Nearest only return one object if there are multiple equidistant results:
+
+        >>> tree = STRtree ([Point(0, 0), Point(0, 0)])
+        >>> tree.nearest(Point(0, 0))
+        0
+        """
+        if self._tree.count == 0:
+            return None
+
+        geometry_arr = np.asarray(geometry, dtype=object)
+        if is_missing(geometry_arr).any() or is_empty(geometry_arr).any():
+            raise ValueError(
+                "Cannot determine nearest geometry for empty geometry or "
+                "missing value (None)."
+            )
+        # _tree.nearest returns ndarray with shape (2, 1) -> index in input
+        # geometries and index into tree geometries
+        indices = self._tree.nearest(np.atleast_1d(geometry_arr))[1]
+
+        if geometry_arr.ndim == 0:
+            return indices[0]
+        else:
+            return indices
+
+    @requires_geos("3.6.0")
+    def query_nearest(
+        self,
+        geometry,
+        max_distance=None,
+        return_distance=False,
+        exclusive=False,
+        all_matches=True,
+    ):
+        """Return the index of the nearest geometries in the tree for each input
+        geometry based on distance within two-dimensional Cartesian space.
+
+        This distance will be 0 when input geometries intersect tree geometries.
+
+        If there are multiple equidistant or intersected geometries in tree and
+        `all_matches` is True (the default), all matching tree geometries are
+        returned; otherwise only the first matching tree geometry is returned.
+        Tree indices are returned in the order they are visited for each input
+        geometry and may not be in ascending index order; no meaningful order is
+        implied.
+
+        The max_distance used to search for nearest items in the tree may have a
+        significant impact on performance by reducing the number of input
+        geometries that are evaluated for nearest items in the tree.  Only those
+        input geometries with at least one tree geometry within +/- max_distance
+        beyond their envelope will be evaluated.  However, using a large
+        max_distance may have a negative performance impact because many tree
+        geometries will be queried for each input geometry.
+
+        The distance, if returned, will be 0 for any intersected geometries in
+        the tree.
+
+        Any geometry that is None or empty in the input geometries is omitted
+        from the output.  Any Z values present in input geometries are ignored
+        when finding nearest tree geometries.
+
+        Parameters
+        ----------
+        geometry : Geometry or array_like
+            Input geometries to query the tree.
+        max_distance : float, optional
+            Maximum distance within which to query for nearest items in tree.
+            Must be greater than 0.
+        return_distance : bool, default False
+            If True, will return distances in addition to indices.
+        exclusive : bool, default False
+            If True, the nearest tree geometries that are equal to the input
+            geometry will not be returned.
+        all_matches : bool, default True
+            If True, all equidistant and intersected geometries will be returned
+            for each input geometry.
+            If False, only the first nearest geometry will be returned.
+
+        Returns
+        -------
+        tree indices or tuple of (tree indices, distances) if geometry is a scalar
+            indices is an ndarray of shape (n, ) and distances (if present) an
+            ndarray of shape (n, )
+
+        OR
+
+        indices or tuple of (indices, distances)
+            indices is an ndarray of shape (2,n) and distances (if present) an
+            ndarray of shape (n).
+            The first subarray of indices contains input geometry indices.
+            The second subarray of indices contains tree geometry indices.
+
+        See also
+        --------
+        nearest: returns singular nearest geometry for each input
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> from shapely import box, Point
+        >>> points = [Point(0, 0), Point(1, 1), Point(2,2), Point(3, 3)]
+        >>> tree = STRtree(points)
+
+        Find the nearest tree geometries to a scalar geometry:
+
+        >>> indices = tree.query_nearest(Point(0.25, 0.25))
+        >>> indices.tolist()
+        [0]
+
+        Retrieve the tree geometries by results of query:
+
+        >>> tree.geometries.take(indices).tolist()
+        [<POINT (0 0)>]
+
+        Find the nearest tree geometries to an array of geometries:
+
+        >>> query_points = np.array([Point(2.25, 2.25), Point(1, 1)])
+        >>> arr_indices = tree.query_nearest(query_points)
+        >>> arr_indices.tolist()
+        [[0, 1], [2, 1]]
+
+        Or transpose to get all pairs of input and tree indices:
+
+        >>> arr_indices.T.tolist()
+        [[0, 2], [1, 1]]
+
+        Retrieve all pairs of input and tree geometries:
+
+        >>> list(zip(query_points.take(arr_indices[0]), tree.geometries.take(arr_indices[1])))
+        [(<POINT (2.25 2.25)>, <POINT (2 2)>), (<POINT (1 1)>, <POINT (1 1)>)]
+
+        All intersecting geometries in the tree are returned by default:
+
+        >>> tree.query_nearest(box(1,1,3,3)).tolist()
+        [1, 2, 3]
+
+        Set all_matches to False to to return a single match per input geometry:
+
+        >>> tree.query_nearest(box(1,1,3,3), all_matches=False).tolist()
+        [1]
+
+        Return the distance to each nearest tree geometry:
+
+        >>> index, distance = tree.query_nearest(Point(0.5, 0.5), return_distance=True)
+        >>> index.tolist()
+        [0, 1]
+        >>> distance.round(4).tolist()
+        [0.7071, 0.7071]
+
+        Return the distance for each input and nearest tree geometry for an array
+        of geometries:
+
+        >>> indices, distance = tree.query_nearest([Point(0.5, 0.5), Point(1, 1)], return_distance=True)
+        >>> indices.tolist()
+        [[0, 0, 1], [0, 1, 1]]
+        >>> distance.round(4).tolist()
+        [0.7071, 0.7071, 0.0]
+
+        Retrieve custom items associated with tree geometries (records can
+        be in whatever data structure so long as geometries and custom data
+        can be extracted into arrays of the same length and order):
+
+        >>> records = [
+        ...     {"geometry": Point(0, 0), "value": "A"},
+        ...     {"geometry": Point(2, 2), "value": "B"}
+        ... ]
+        >>> tree = STRtree([record["geometry"] for record in records])
+        >>> items = np.array([record["value"] for record in records])
+        >>> items.take(tree.query_nearest(Point(0.5, 0.5))).tolist()
+        ['A']
+        """
+
+        geometry = np.asarray(geometry, dtype=object)
+        is_scalar = False
+        if geometry.ndim == 0:
+            geometry = np.expand_dims(geometry, 0)
+            is_scalar = True
+
+        if max_distance is not None:
+            if not np.isscalar(max_distance):
+                raise ValueError("max_distance parameter only accepts scalar values")
+
+            if max_distance <= 0:
+                raise ValueError("max_distance must be greater than 0")
+
+        # a distance of 0 means no max_distance is used
+        max_distance = max_distance or 0
+
+        if not np.isscalar(exclusive):
+            raise ValueError("exclusive parameter only accepts scalar values")
+
+        if exclusive not in {True, False}:
+            raise ValueError("exclusive parameter must be boolean")
+
+        if not np.isscalar(all_matches):
+            raise ValueError("all_matches parameter only accepts scalar values")
+
+        if all_matches not in {True, False}:
+            raise ValueError("all_matches parameter must be boolean")
+
+        results = self._tree.query_nearest(
+            geometry, max_distance, exclusive, all_matches
+        )
+
+        # output indices are shape (n, )
+        if is_scalar:
+            if not return_distance:
+                return results[0][1]
+
+            else:
+                return (results[0][1], results[1])
+
+        # output indices are shape (2, n)
+        if not return_distance:
+            return results[0]
+
+        return results

billinglayer/python/shapely/testing.py

@@ -0,0 +1,204 @@
+from functools import partial
+
+import numpy as np
+
+import shapely
+
+__all__ = ["assert_geometries_equal"]
+
+
+def _equals_exact_with_ndim(x, y, tolerance):
+    dimension_equals = shapely.get_coordinate_dimension(
+        x
+    ) == shapely.get_coordinate_dimension(y)
+    with np.errstate(invalid="ignore"):
+        # Suppress 'invalid value encountered in equals_exact' with nan coordinates
+        geometry_equals = shapely.equals_exact(x, y, tolerance=tolerance)
+    return dimension_equals & geometry_equals
+
+
+def _replace_nan(arr):
+    return np.where(np.isnan(arr), 0.0, arr)
+
+
+def _assert_nan_coords_same(x, y, tolerance, err_msg, verbose):
+    x, y = np.broadcast_arrays(x, y)
+    x_coords = shapely.get_coordinates(x, include_z=True)
+    y_coords = shapely.get_coordinates(y, include_z=True)
+
+    # Check the shapes (condition is copied from numpy test_array_equal)
+    if x_coords.shape != y_coords.shape:
+        return False
+
+    # Check NaN positional equality
+    x_id = np.isnan(x_coords)
+    y_id = np.isnan(y_coords)
+    if not (x_id == y_id).all():
+        msg = build_err_msg(
+            [x, y],
+            err_msg + "\nx and y nan coordinate location mismatch:",
+            verbose=verbose,
+        )
+        raise AssertionError(msg)
+
+    # If this passed, replace NaN with a number to be able to use equals_exact
+    x_no_nan = shapely.transform(x, _replace_nan, include_z=True)
+    y_no_nan = shapely.transform(y, _replace_nan, include_z=True)
+
+    return _equals_exact_with_ndim(x_no_nan, y_no_nan, tolerance=tolerance)
+
+
+def _assert_none_same(x, y, err_msg, verbose):
+    x_id = shapely.is_missing(x)
+    y_id = shapely.is_missing(y)
+
+    if not (x_id == y_id).all():
+        msg = build_err_msg(
+            [x, y],
+            err_msg + "\nx and y None location mismatch:",
+            verbose=verbose,
+        )
+        raise AssertionError(msg)
+
+    # If there is a scalar, then here we know the array has the same
+    # flag as it everywhere, so we should return the scalar flag.
+    if x.ndim == 0:
+        return bool(x_id)
+    elif y.ndim == 0:
+        return bool(y_id)
+    else:
+        return y_id
+
+
+def assert_geometries_equal(
+    x,
+    y,
+    tolerance=1e-7,
+    equal_none=True,
+    equal_nan=True,
+    normalize=False,
+    err_msg="",
+    verbose=True,
+):
+    """Raises an AssertionError if two geometry array_like objects are not equal.
+
+    Given two array_like objects, check that the shape is equal and all elements of
+    these objects are equal. An exception is raised at shape mismatch or conflicting
+    values. In contrast to the standard usage in shapely, no assertion is raised if
+    both objects have NaNs/Nones in the same positions.
+
+    Parameters
+    ----------
+    x : Geometry or array_like
+    y : Geometry or array_like
+    equal_none : bool, default True
+        Whether to consider None elements equal to other None elements.
+    equal_nan : bool, default True
+        Whether to consider nan coordinates as equal to other nan coordinates.
+    normalize : bool, default False
+        Whether to normalize geometries prior to comparison.
+    err_msg : str, optional
+        The error message to be printed in case of failure.
+    verbose : bool, optional
+        If True, the conflicting values are appended to the error message.
+    """
+    __tracebackhide__ = True  # Hide traceback for py.test
+    if normalize:
+        x = shapely.normalize(x)
+        y = shapely.normalize(y)
+    x = np.array(x, copy=False)
+    y = np.array(y, copy=False)
+
+    is_scalar = x.ndim == 0 or y.ndim == 0
+
+    # Check the shapes (condition is copied from numpy test_array_equal)
+    if not (is_scalar or x.shape == y.shape):
+        msg = build_err_msg(
+            [x, y],
+            err_msg + f"\n(shapes {x.shape}, {y.shape} mismatch)",
+            verbose=verbose,
+        )
+        raise AssertionError(msg)
+
+    flagged = False
+    if equal_none:
+        flagged = _assert_none_same(x, y, err_msg, verbose)
+
+    if not np.isscalar(flagged):
+        x, y = x[~flagged], y[~flagged]
+        # Only do the comparison if actual values are left
+        if x.size == 0:
+            return
+    elif flagged:
+        # no sense doing comparison if everything is flagged.
+        return
+
+    is_equal = _equals_exact_with_ndim(x, y, tolerance=tolerance)
+    if is_scalar and not np.isscalar(is_equal):
+        is_equal = bool(is_equal[0])
+
+    if np.all(is_equal):
+        return
+    elif not equal_nan:
+        msg = build_err_msg(
+            [x, y],
+            err_msg + f"\nNot equal to tolerance {tolerance:g}",
+            verbose=verbose,
+        )
+        raise AssertionError(msg)
+
+    # Optionally refine failing elements if NaN should be considered equal
+    if not np.isscalar(is_equal):
+        x, y = x[~is_equal], y[~is_equal]
+        # Only do the NaN check if actual values are left
+        if x.size == 0:
+            return
+    elif is_equal:
+        # no sense in checking for NaN if everything is equal.
+        return
+
+    is_equal = _assert_nan_coords_same(x, y, tolerance, err_msg, verbose)
+    if not np.all(is_equal):
+        msg = build_err_msg(
+            [x, y],
+            err_msg + f"\nNot equal to tolerance {tolerance:g}",
+            verbose=verbose,
+        )
+        raise AssertionError(msg)
+
+
+## BELOW A COPY FROM numpy.testing._private.utils (numpy version 1.20.2)
+
+
+def build_err_msg(
+    arrays,
+    err_msg,
+    header="Geometries are not equal:",
+    verbose=True,
+    names=("x", "y"),
+    precision=8,
+):
+    msg = ["\n" + header]
+    if err_msg:
+        if err_msg.find("\n") == -1 and len(err_msg) < 79 - len(header):
+            msg = [msg[0] + " " + err_msg]
+        else:
+            msg.append(err_msg)
+    if verbose:
+        for i, a in enumerate(arrays):
+
+            if isinstance(a, np.ndarray):
+                # precision argument is only needed if the objects are ndarrays
+                r_func = partial(np.array_repr, precision=precision)
+            else:
+                r_func = repr
+
+            try:
+                r = r_func(a)
+            except Exception as exc:
+                r = f"[repr failed for <{type(a).__name__}>: {exc}]"
+            if r.count("\n") > 3:
+                r = "\n".join(r.splitlines()[:3])
+                r += "..."
+            msg.append(f" {names[i]}: {r}")
+    return "\n".join(msg)