Coverage for src/sensai/geoanalytics/geopandas/geometry.py: 31%

1import logging 

2import numpy as np 

3from itertools import combinations 

4from scipy.spatial.distance import euclidean 

5from scipy.spatial.qhull import Delaunay 

6from shapely.geometry import MultiLineString, Polygon 

7from shapely.ops import polygonize, unary_union 

8 

9from .coordinates import extract_coordinates_array, TCoordinates 

10from .graph import delaunay_graph 

11 

12log = logging.getLogger(__name__) 

13 

14 

15# after already having implemented this, I found the following package: https://github.com/bellockk/alphashape 

16# It will compute the same polygons (I have verified it). It also contains an optimizer for alpha, which is, however, 

17# extremely slow and therefore unusable in most practical applications. 

18def alpha_shape(coordinates: TCoordinates, alpha=0.5): 

19 """ 

20 Compute the `alpha shape`_ of a set of points. Based on `this implementation`_. In contrast to the standard 

21 definition of the parameter alpha here we normalize it by the mean edge size of the cluster. This results in 

22 similar "concavity properties" of the resulting shapes for different coordinate sets and a fixed alpha. 

23 

24 .. _this implementation: https://sgillies.net/2012/10/13/the-fading-shape-of-alpha.html 

25 .. _alpha shape: https://en.wikipedia.org/wiki/Alpha_shape 

26 

27 :param coordinates: a suitable iterable of 2-dimensional coordinates 

28 :param alpha: alpha value to influence the gooeyness of the border. Larger numbers 

29 don't fall inward as much as smaller numbers. 

30 :return: a shapely Polygon 

31 """ 

32 coordinates = extract_coordinates_array(coordinates) 

33 

34 edge_index_pairs = set() 

35 edge_vertex_pairs = [] 

36 graph = delaunay_graph(coordinates) 

37 mean_edge_size = graph.size(weight="weight") / graph.number_of_edges() 

38 

39 def add_edge(edge_index_pair, edge_vertex_pair): 

40 """ 

41 Add a line between the i-th and j-th points, 

42 if not in the list already 

43 """ 

44 edge_index_pair = tuple(sorted(edge_index_pair)) 

45 if edge_index_pair in edge_index_pairs: 

46 # already added 

47 return 

48 edge_index_pairs.add(edge_index_pair) 

49 edge_vertex_pairs.append(edge_vertex_pair) 

50 

51 tri = Delaunay(coordinates) 

52 for simplex in tri.simplices: 

53 vertices = tri.points[simplex] 

54 area = Polygon(vertices).area 

55 edges = combinations(vertices, 2) 

56 product_edges_lengths = 1 

57 for vertex_1, vertex_2 in edges: 

58 product_edges_lengths *= euclidean(vertex_1, vertex_2) 

59 # this is the radius of the circumscribed circle of the triangle 

60 # see https://en.wikipedia.org/wiki/Circumscribed_circle#Triangles 

61 circum_r = product_edges_lengths / (4.0 * area) 

62 

63 if circum_r < mean_edge_size/alpha: 

64 for index_pair in combinations(simplex, 2): 

65 add_edge(index_pair, tri.points[np.array(index_pair)]) 

66 

67 remaining_edges = MultiLineString(edge_vertex_pairs) 

68 

69 return unary_union(list(polygonize(remaining_edges)))