"""
This module implements some of the spatial analysis techniques and processes used to
understand the patterns and relationships of geographic features.
This is mainly inspired by turf.js.
link: http://turfjs.org/
"""
import copy
import itertools
import math
from math import floor, sqrt
from typing import List, Optional, Union
import numpy as np
from geojson import Feature, FeatureCollection, LineString, MultiLineString
from geojson import Point as GeoPoint
from geojson import Polygon
from scipy.spatial import Delaunay, Voronoi
from shapely import geometry as geometry
from shapely.geometry import LineString as ShapelyLineString
from shapely.geometry import MultiPoint, MultiPolygon, Point, mapping, shape
from shapely.ops import cascaded_union, clip_by_rect, polygonize, unary_union
from turfpy.helper import get_coord, get_coords, get_geom, get_type, length_to_degrees
from turfpy.measurement import (
bbox,
bbox_polygon,
center,
centroid,
destination,
rhumb_bearing,
rhumb_destination,
rhumb_distance,
)
from turfpy.meta import coord_each, feature_each, flatten_each
from .dev_lib.earcut import earcut
from .dev_lib.spline import Spline
[docs]def circle(
center: Feature, radius: int, steps: int = 64, units: str = "km", **kwargs
) -> Polygon:
"""
Takes a Point and calculates the circle polygon given a radius in degrees,
radians, miles, or kilometers; and steps for precision.
:param center: A `Point` object representing center point of circle.
:param radius: An int representing radius of the circle.
:param steps: An int representing number of steps.
:param units: A string representing units of distance e.g. 'mi', 'km',
'deg' and 'rad'.
:param kwargs: A dict representing additional properties.
:return: A polygon feature object.
Example:
>>> from turfpy.transformation import circle
>>> from geojson import Feature, Point
>>> circle(center=Feature(geometry=Point((-75.343, 39.984))), radius=5, steps=10)
"""
coordinates = []
options = dict(steps=steps, units=units)
options.update(kwargs)
for i in range(steps):
bearing = i * -360 / steps
pt = destination(center, radius, bearing, options=options)
cords = pt.geometry.coordinates
coordinates.append(cords)
coordinates.append(coordinates[0])
return Feature(geometry=Polygon([coordinates], **kwargs))
[docs]def bbox_clip(geojson: Feature, bbox: list) -> Feature:
"""
Takes a Feature or geometry and a bbox and clips the feature to the bbox
:param geojson: Geojson data
:param bbox: Bounding Box which is used to clip the geojson
:return: Clipped geojson
Example:
>>> from turfpy.transformation import bbox_clip
>>> from geojson import Feature
>>> f = Feature(geometry={"coordinates": [[[2, 2], [8, 4],
>>> [12, 8], [3, 7], [2, 2]]], "type": "Polygon"})
>>> bbox = [0, 0, 10, 10]
>>> clip = bbox_clip(f, bbox)
"""
bb_polygon = bbox_polygon(bbox)
bb_clip = intersect([geojson, bb_polygon])
if not bb_clip:
return bb_clip
if "properties" in geojson:
bb_clip.properties = geojson["properties"]
return bb_clip
[docs]def intersect(features: Union[List[Feature], FeatureCollection]) -> Feature:
"""
Takes polygons and finds their intersection
:param features: List of features of Feature Collection
:return: Intersection Geojson Feature
Example:
>>> from turfpy.transformation import intersect
>>> from geojson import Feature
>>> f = Feature(geometry={"coordinates": [
>>> [[-122.801742, 45.48565], [-122.801742, 45.60491],
>>> [-122.584762, 45.60491], [-122.584762, 45.48565],
>>> [-122.801742, 45.48565]]], "type": "Polygon"})
>>> b = Feature(geometry={"coordinates": [
>>> [[-122.520217, 45.535693], [-122.64038, 45.553967],
>>> [-122.720031, 45.526554], [-122.669906, 45.507309],
>>> [-122.723464, 45.446643], [-122.532577, 45.408574],
>>> [-122.487258, 45.477466], [-122.520217, 45.535693]
>>> ]], "type": "Polygon"})
>>> inter = intersect([f, b])
"""
properties_list = []
if isinstance(features, list):
shapes = []
for f in features:
poly = get_geom(f)
s = shape(poly)
shapes.append(s)
if "properties" in f.keys():
properties_list.append(f["properties"])
else:
if "features" not in features.keys():
raise Exception("Invalid FeatureCollection")
if "properties" in features.keys():
properties_list.append(features["properties"])
shapes = []
for f in features["features"]:
poly = get_geom(f)
s = shape(poly)
shapes.append(s)
if "properties" in f.keys():
properties_list.append(f["properties"])
intersection = shapes[0]
for shape_value in shapes:
intersection = shape_value.intersection(intersection)
intersection = mapping(intersection)
if (
len(intersection.get("coordinates", [])) == 0
and len(intersection.get("geometries", [])) == 0
):
return None
properties = merge_dict(properties_list)
intersection_feature = Feature(geometry=intersection, properties=properties)
return intersection_feature
[docs]def bezier_spline(line: Feature, resolution=10000, sharpness=0.85) -> Feature:
"""
Takes a line and returns a curved version by applying a Bezier spline algorithm
:param line: LineString Feature which is used to draw the curve
:param resolution: time in milliseconds between points
:param sharpness: a measure of how curvy the path should be between splines
:return: Curve as LineString Feature
Example:
>>> from geojson import LineString, Feature
>>> from turfpy.transformation import bezier_spline
>>> ls = LineString([(-76.091308, 18.427501),
>>> (-76.695556, 18.729501),
>>> (-76.552734, 19.40443),
>>> (-74.61914, 19.134789),
>>> (-73.652343, 20.07657),
>>> (-73.157958, 20.210656)])
>>> f = Feature(geometry=ls)
>>> bezier_spline(f)
"""
coords = []
points = []
geom = get_geom(line)
for c in geom["coordinates"]:
points.append({"x": c[0], "y": c[1]})
spline = Spline(points_data=points, resolution=resolution, sharpness=sharpness)
i = 0
while i < spline.duration:
pos = spline.pos(i)
if floor(i / 100) % 2 == 0:
coords.append((pos["x"], pos["y"]))
i = i + 10
return Feature(geometry=LineString(coords))
[docs]def merge_dict(dicts: list):
super_dict: dict = {}
for d in dicts:
for k, v in d.items():
if k not in super_dict.keys():
super_dict[k] = v
else:
if isinstance(super_dict[k], list):
if v not in super_dict[k]:
super_dict[k].append(v)
else:
if super_dict[k] != v:
super_dict[k] = [super_dict[k], v]
return super_dict
[docs]def union(
features: Union[List[Feature], FeatureCollection]
) -> Union[Feature, FeatureCollection]:
"""
Given list of features or ``FeatureCollection`` return union of those.
:param features: A list of GeoJSON features or FeatureCollection.
:return: A GeoJSON Feature or FeatureCollection.
Example:
>>> from turfpy.transformation import union
>>> from geojson import Feature, Polygon, FeatureCollection
>>> f1 = Feature(geometry=Polygon([[
... [-82.574787, 35.594087],
... [-82.574787, 35.615581],
... [-82.545261, 35.615581],
... [-82.545261, 35.594087],
... [-82.574787, 35.594087]
... ]]), properties={"fill": "#00f"})
>>> f2 = Feature(geometry=Polygon([[
... [-82.560024, 35.585153],
... [-82.560024, 35.602602],
... [-82.52964, 35.602602],
... [-82.52964, 35.585153],
... [-82.560024, 35.585153]]]), properties={"fill": "#00f"})
>>> union(FeatureCollection([f1, f2], properties={"combine": "yes"}))
"""
shapes = []
properties_list = []
if isinstance(features, list):
for f in features:
if f.type != "Feature":
raise Exception("Not a valid feature")
geom = get_geom(f)
s = shape(geom)
shapes.append(s)
if "properties" in f.keys():
properties_list.append(f["properties"])
else:
if "features" not in features.keys():
raise Exception("Invalid FeatureCollection")
if "properties" in features.keys():
properties_list.append(features["properties"])
for f in features["features"]:
geom = get_geom(f)
s = shape(geom)
shapes.append(s)
if "properties" in f.keys():
properties_list.append(f["properties"])
result = cascaded_union(shapes)
result = mapping(result)
properties = merge_dict(properties_list)
if result["type"] == "GeometryCollection":
features = []
for geom in result["geometries"]:
features.append(Feature(geometry=geom))
return FeatureCollection(features, properties=properties)
return Feature(geometry=result, properties=properties)
[docs]def _alpha_shape(points, alpha):
"""
Compute the alpha shape (concave hull) of a set of points.
:param points: Iterable container of points.
:param alpha: alpha value to influence the gooeyness of the border. Smaller
numbers don't fall inward as much as larger numbers. Too large,
and you lose everything!
"""
if len(points) < 4:
# When you have a triangle, there is no sense in computing an alpha
# shape.
return geometry.MultiPoint(list(points)).convex_hull
def add_edge(edges, edge_points, coords, i, j):
"""Add a line between the i-th and j-th points, if not in the list already"""
if (i, j) in edges or (j, i) in edges:
# already added
return
edges.add((i, j))
edge_points.append(coords[[i, j]])
coords = np.array([point.coords[0] for point in points])
tri = Delaunay(coords)
edges = set()
edge_points = []
# loop over triangles:
# ia, ib, ic = indices of corner points of the triangle
for ia, ib, ic in tri.vertices:
pa = coords[ia]
pb = coords[ib]
pc = coords[ic]
# Lengths of sides of triangle
a = math.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
b = math.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
c = math.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
# Semiperimeter of triangle
s = (a + b + c) / 2.0
# Area of triangle by Heron's formula
area = math.sqrt(s * (s - a) * (s - b) * (s - c))
if area > 0:
circum_r = a * b * c / (4.0 * area)
else:
circum_r = 0
# Here's the radius filter.
# print circum_r
if circum_r < 1.0 / alpha:
add_edge(edges, edge_points, coords, ia, ib)
add_edge(edges, edge_points, coords, ib, ic)
add_edge(edges, edge_points, coords, ic, ia)
m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))
return cascaded_union(triangles), edge_points
[docs]def get_points(features):
points = []
if "type" not in features.keys():
raise Exception("Invalid Feature")
if features["type"] == "Feature":
get_ext_points(geometry.shape(features["geometry"]), points)
else:
if "features" not in features.keys():
raise Exception("Invalid FeatureCollection")
for feature in features["features"]:
get_ext_points(geometry.shape(feature["geometry"]), points)
return points
[docs]def get_ext_points(geom, points):
if geom.type == "Point":
for p in geom.coords:
points.append(Point(p))
elif geom.type == "MultiPoint":
for p in geom.geoms:
points.append(p)
elif geom.type == "LineString":
for p in geom.coords:
points.append(Point(p))
elif geom.type == "MultiLineString":
for g in geom.geoms:
for p in g.coords:
points.append(Point(p))
elif geom.type == "Polygon":
for p in geom.exterior.coords:
points.append(Point(p))
elif geom.type == "MultiPolygon":
for g in geom.geoms:
for p in g.exterior.coords:
points.append(Point(p))
else:
raise Exception("Invalid Geometry")
[docs]def concave(features: Union[Feature, FeatureCollection], alpha=2):
"""Generate concave hull for the given feature or Feature Collection.
:param features: It can be a feature or Feature Collection
:param alpha: Alpha determines the shape of concave hull,
greater values will make shape more tighten
:return: Feature of concave hull polygon
Example:
>>> from turfpy.transformation import concave
>>> from geojson import FeatureCollection, Feature, Point
>>> f1 = Feature(geometry=Point((-63.601226, 44.642643)))
>>> f2 = Feature(geometry=Point((-63.591442, 44.651436)))
>>> f3 = Feature(geometry=Point((-63.580799, 44.648749)))
>>> f4 = Feature(geometry=Point((-63.573589, 44.641788)))
>>> f5 = Feature(geometry=Point((-63.587665, 44.64533)))
>>> f6 = Feature(geometry=Point((-63.595218, 44.64765)))
>>> fc = [f1, f2, f3, f4, f5, f6]
>>> concave(FeatureCollection(fc), alpha=100)
"""
points = get_points(features)
concave_hull, edges = _alpha_shape(points, alpha)
return Feature(geometry=mapping(concave_hull))
[docs]def convex(features: Union[Feature, FeatureCollection]):
"""Generate convex hull for the given feature or Feature Collection
:param features: It can be a feature or Feature Collection
:return: Feature of convex hull polygon
Example:
>>> from turfpy.transformation import convex
>>> from geojson import FeatureCollection, Feature, Point
>>> f1 = Feature(geometry=Point((10.195312, 43.755225)))
>>> f2 = Feature(geometry=Point((10.404052, 43.8424511)))
>>> f3 = Feature(geometry=Point((10.579833, 43.659924)))
>>> f4 = Feature(geometry=Point((10.360107, 43.516688)))
>>> f5 = Feature(geometry=Point((10.14038, 43.588348)))
>>> f6 = Feature(geometry=Point((10.195312, 43.755225)))
>>> fc = [f1, f2, f3, f4, f5, f6]
>>> convex(FeatureCollection(fc))
"""
points = get_points(features)
point_collection = geometry.MultiPoint(list(points))
return Feature(geometry=mapping(point_collection.convex_hull))
[docs]def dissolve(
features: Union[List[Feature], FeatureCollection], property_name: str = None
) -> FeatureCollection:
"""
Take FeatureCollection or list of features to dissolve based on
property_name provided.
:param features: A list of GeoJSON features or FeatureCollection.
:param property_name: Name of property based on which to dissolve.
:return: A GeoJSON Feature or FeatureCollection.
Example:
>>> from geojson import Polygon, Feature, FeatureCollection
>>> from turfpy.transformation import dissolve
>>> f1 = Feature(geometry=Polygon([[
>>> [0, 0],
>>> [0, 1],
>>> [1, 1],
>>> [1, 0],
>>> [0, 0]]]), properties={"combine": "yes", "fill": "#00f"})
>>> f2 = Feature(geometry=Polygon([[
>>> [0, -1],
>>> [0, 0],
>>> [1, 0],
>>> [1, -1],
>>> [0,-1]]]), properties={"combine": "yes"})
>>> f3 = Feature(geometry=Polygon([[
>>> [1,-1],
>>> [1, 0],
>>> [2, 0],
>>> [2, -1],
>>> [1, -1]]]), properties={"combine": "no"})
>>> dissolve(FeatureCollection([f1, f2, f3]), property_name='combine')
"""
if isinstance(features, list):
features = FeatureCollection(features)
if "features" not in features.keys():
raise Exception("Invalid FeatureCollection")
dissolve_feature_list = []
if property_name:
for k, g in itertools.groupby(
features["features"], key=lambda x: x["properties"].get(property_name)
):
fc = FeatureCollection(list(g))
# if "properties" in features.keys():
# fc['properties'] = features['properties']
result = union(fc)
if result["type"] == "FeatureCollection":
for f in result["features"]:
dissolve_feature_list.append(f)
else:
dissolve_feature_list.append(result)
else:
return union(features)
if "properties" in features.keys():
return FeatureCollection(dissolve_feature_list, properties=features["properties"])
else:
return FeatureCollection(dissolve_feature_list)
[docs]def difference(feature_1: Feature, feature_2: Feature) -> Feature:
"""
Find the difference between given two features.
:param feature_1: A GeoJSON feature
:param feature_2: A GeoJSON feature
:return: A GeoJSON feature
Example:
>>> from geojson import Polygon, Feature
>>> from turfpy.transformation import difference
>>> f1 = Feature(geometry=Polygon([[
>>> [128, -26],
>>> [141, -26],
>>> [141, -21],
>>> [128, -21],
>>> [128, -26]]]), properties={"combine": "yes", "fill": "#00f"})
>>> f2 = Feature(geometry=Polygon([[
>>> [126, -28],
>>> [140, -28],
>>> [140, -20],
>>> [126, -20],
>>> [126, -28]]]), properties={"combine": "yes"})
>>> difference(f1, f2)
"""
properties_list = []
if "properties" in feature_1.keys():
properties_list.append(feature_1["properties"])
if "properties" in feature_2.keys():
properties_list.append(feature_2["properties"])
shape_1 = shape(get_geom(feature_1))
shape_2 = shape(get_geom(feature_2))
difference_result = shape_1.difference(shape_2)
difference_result = mapping(difference_result)
if len(difference_result["coordinates"]) == 0:
return None
properties = merge_dict(properties_list)
difference_feature = Feature(geometry=difference_result, properties=properties)
return difference_feature
[docs]def scale(feature, factor, origin):
is_point = get_type(feature) == "Point"
origin = define_origin(feature, origin)
if factor == 1 or is_point:
return feature
def _callback_coord_each(
coord, coord_index, feature_index, multi_feature_index, geometry_index
):
nonlocal factor, origin
original_distance = rhumb_distance(GeoPoint(origin), GeoPoint(coord))
bearing = rhumb_bearing(GeoPoint(origin), GeoPoint(coord))
new_distance = original_distance * factor
new_coord = get_coord(rhumb_destination(GeoPoint(origin), new_distance, bearing))
coord[0] = new_coord[0]
coord[1] = new_coord[1]
if len(coord) == 3:
coord[2] = coord[2] * factor
coord_each(feature, _callback_coord_each)
return feature
[docs]def define_origin(geojson, origin):
if not origin:
origin = "centroid"
if isinstance(origin, list):
return get_coord(origin)
bb = bbox(geojson)
west = bb[0]
south = bb[1]
east = bb[2]
north = bb[3]
if (
origin == "sw"
or origin == "southwest"
or origin == "westsouth"
or origin == "bottomleft"
):
return [west, south]
elif (
origin == "se"
or origin == "southeast"
or origin == "eastsouth"
or origin == "bottomright"
):
return [east, south]
elif (
origin == "nw"
or origin == "northwest"
or origin == "westnorth"
or origin == "topleft"
):
return [west, north]
elif (
origin == "ne"
or origin == "northeast"
or origin == "eastnorth"
or origin == "topright"
):
return [east, north]
elif origin == "center":
return center(geojson)["geometry"]["coordinates"]
elif origin is None or origin == "centroid":
return centroid(geojson)["geometry"]["coordinates"]
else:
raise Exception("invalid origin")
[docs]def tesselate(poly: Feature) -> FeatureCollection:
"""Tesselates a Feature into a FeatureCollection of triangles using earcut.
:param poly: A GeoJSON feature ``class:geojson.Polygon``.
:return: A GeoJSON FeatureCollection of triangular polygons.
Example:
>>> from geojson import Feature
>>> from turfpy.transformation import tesselate
>>> polygon = Feature(geometry={"coordinates": [[[11, 0], [22, 4], [31, 0], [31, 11],
... [21, 15], [11, 11], [11, 0]]], "type": "Polygon"})
>>> tesselate(polygon)
"""
if (
poly["geometry"]["type"] != "Polygon"
and poly["geometry"]["type"] != "MultiPolygon"
):
raise ValueError("Geometry must be Polygon or MultiPolygon")
fc = FeatureCollection([])
if poly.geometry.type == "Polygon":
fc["features"] = __process_polygon(poly.geometry.coordinates)
else:
for co in poly.geometry.coordinates:
fc["features"].extend(__process_polygon(co))
return fc
[docs]def __process_polygon(coordinates):
data = __flatten_coords(coordinates)
dim = 2
result = earcut(data["vertices"], data["holes"], dim)
features = []
vertices = []
for i, val in enumerate(result):
index = val
vertices.append(
[data["vertices"][index * dim], data["vertices"][index * dim + 1]]
)
i = 0
while i < len(vertices):
coords = vertices[i : i + 3]
coords.append(vertices[i])
features.append(Feature(geometry={"coordinates": [coords], "type": "Polygon"}))
i += 3
return features
[docs]def __flatten_coords(data):
dim = len(data[0][0])
result = {"vertices": [], "holes": [], "dimensions": dim}
hole_index = 0
for i, val in enumerate(data):
for j, _ in enumerate(val):
for d in range(dim):
result["vertices"].append(data[i][j][d])
if i > 0:
hole_index += len(data[i - 1])
result["holes"].append(hole_index)
return result
[docs]def line_offset(geojson: Feature, distance: float, unit: str = "km") -> Feature:
"""
Takes a linestring or multilinestring and returns
a line at offset by the specified distance.
:param geojson: input GeoJSON
:param distance: distance to offset the line (can be of negative value)
:param unit: Units in which distance to be calculated, values can be 'deg', 'rad',
'mi', 'km', default is 'km'
:return: Line feature offset from the input line
Example:
>>> from geojson import MultiLineString, Feature
>>> from turfpy.transformation import line_offset
>>> ls = Feature(geometry=MultiLineString([
... [(3.75, 9.25), (-130.95, 1.52)],
... [(23.15, -34.25), (-1.35, -4.65), (3.45, 77.95)]
... ]))
>>> line_offset(ls, 2, unit='mi')
"""
if not geojson:
raise Exception("geojson is required")
if not distance:
raise Exception("distance is required")
type = get_type(geojson)
properties = geojson.get("properties", {})
if type == "LineString":
return line_offset_feature(geojson, distance, unit)
elif type == "MultiLineString":
coords = []
def callback_flatten_each(feature, feature_index, multi_feature_index):
nonlocal coords
coords.append(
line_offset_feature(feature, distance, unit).geometry.coordinates
)
return True
flatten_each(geojson, callback_flatten_each)
return Feature(geometry=MultiLineString(coords), properties=properties)
[docs]def line_offset_feature(line, distance, units):
segments = []
offset_degrees = length_to_degrees(distance, units)
coords = get_coords(line)
final_coords = []
for index, current_coords in enumerate(coords):
if index != len(coords) - 1:
segment = _process_segment(current_coords, coords[index + 1], offset_degrees)
segments.append(segment)
if index > 0:
seg2_coords = segments[index - 1]
intersects = _intersection(segment, seg2_coords)
if intersects:
seg2_coords[1] = intersects
segment[0] = intersects
final_coords.append(seg2_coords[0])
if index == len(coords) - 2:
final_coords.append(segment[0])
final_coords.append(segment[1])
if len(coords) == 2:
final_coords.append(segment[0])
final_coords.append(segment[1])
return Feature(
geometry=LineString(final_coords), properties=line.get("properties", {})
)
[docs]def _process_segment(point1, point2, offset):
L = sqrt(
(point1[0] - point2[0]) * (point1[0] - point2[0])
+ (point1[1] - point2[1]) * (point1[1] - point2[1])
)
out1x = point1[0] + offset * (point2[1] - point1[1]) / L
out2x = point2[0] + offset * (point2[1] - point1[1]) / L
out1y = point1[1] + offset * (point1[0] - point2[0]) / L
out2y = point2[1] + offset * (point1[0] - point2[0]) / L
return [[out1x, out1y], [out2x, out2y]]
[docs]def _intersection(a, b):
if _is_parallel(a, b):
return False
return _intersect_segments(a, b)
[docs]def _is_parallel(a, b):
r = _ab(a)
s = _ab(b)
return _cross_product(r, s) == 0
[docs]def _ab(segment):
start = segment[0]
end = segment[1]
return [end[0] - start[0], end[1] - start[1]]
[docs]def _cross_product(v1, v2):
return (v1[0] * v2[1]) - (v2[0] * v1[1])
[docs]def _intersect_segments(a, b):
p = a[0]
r = _ab(a)
q = b[0]
s = _ab(b)
cross = _cross_product(r, s)
qmp = _sub(q, p)
numerator = _cross_product(qmp, s)
t = numerator / cross
intersection = _add(p, _scalar_mult(t, r))
return intersection
[docs]def _add(v1, v2):
return [v1[0] + v2[0], v1[1] + v2[1]]
[docs]def _sub(v1, v2):
return [v1[0] - v2[0], v1[1] - v2[1]]
[docs]def _scalar_mult(s, v):
return [s * v[0], s * v[1]]
[docs]def voronoi(
points: Union[FeatureCollection, List], bbox: Optional[list] = None
) -> Feature:
"""Takes a FeatureCollection of points, and a bounding box,
and returns a FeatureCollection of Voronoi polygons.
:param points: To find the Voronoi polygons around. Points should be either
FeatureCollection of points or list of points.
:param bbox: A bounding box to clip.
:return: A GeoJSON Feature.
Example:
>>> from turfpy.transformation import voronoi
>>> points = [
... [-66.9703, 40.3183],
... [-63.7763, 40.4500],
... [-65.4196, 42.13985310302137],
... [-69.5813, 43.95405461286195],
... [-65.66337553550034, 55.97088945355232],
... [-60.280418548905, 56.240669185466146],
... [-68.5129561347689, 50.12984589640148],
... [-64.2393519226657, 59.66235385923687],
... ]
>>> bbox = [-70, 40, -60, 60]
>>> voronoi(points, bbox)
"""
if isinstance(points, FeatureCollection):
coords = []
for feature in points["features"]:
coords.append(feature["features"][0]["geometry"]["coordinates"])
points = np.array(coords)
elif isinstance(points, list):
points = np.array(points)
else:
raise ValueError(
"points should be either FeatureCollection of points of List of Points"
)
vor = Voronoi(points)
lines = [
ShapelyLineString(vor.vertices[line])
for line in vor.ridge_vertices
if -1 not in line
]
convex_hull = MultiPoint([Point(i) for i in points]).convex_hull.buffer(2)
result = MultiPolygon([poly.intersection(convex_hull) for poly in polygonize(lines)])
result = MultiPolygon(
[p for p in result] + [p for p in convex_hull.difference(unary_union(result))]
)
if bbox is not None:
w, s, e, n = bbox
cliped_result = clip_by_rect(result, w, s, e, n)
return Feature(geometry=cliped_result)
return Feature(geometry=result)