pythonGISToolbox/centerline.py at master · gvalenciam/pythonGISToolbox · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
# -*- encoding: utf-8 -*-
from shapely.geometry import LineString
from shapely.ops import unary_union
from scipy.spatial import Voronoi
import numpy as np


class Centerline(object):

	def __init__(self, inputGEOM, dist=0.5):
		self.inputGEOM = inputGEOM
		self.dist = abs(dist)

	def createCenterline(self):
		"""
		Calculates the centerline of a polygon.

		Densifies the border of a polygon which is then represented by a Numpy
		array of points necessary for creating the Voronoi diagram. Once the
		diagram is created, the ridges located within the polygon are
		joined and returned.

		Returns:
			a union of lines that are located within the polygon.
		"""

		minx = int(min(self.inputGEOM.envelope.exterior.xy[0]))
		miny = int(min(self.inputGEOM.envelope.exterior.xy[1]))

		border = np.array(self.densifyBorder(self.inputGEOM, minx, miny))

		# print("BORDER")
		# print(border)

		vor = Voronoi(border)
		vertex = vor.vertices

		lst_lines = []
		for j, ridge in enumerate(vor.ridge_vertices):
			if -1 not in ridge:
				line = LineString([ \
					(vertex[ridge[0]][0] + minx, vertex[ridge[0]][1] + miny), \
					(vertex[ridge[1]][0] + minx, vertex[ridge[1]][1] + miny)])

				if line.within(self.inputGEOM) and len(line.coords[0]) > 1:
					lst_lines.append(line)

		return unary_union(lst_lines)

		# if len(border) > 3:
		# 	vor = Voronoi(border)
		# 	vertex = vor.vertices
		#
		# 	lst_lines = []
		# 	for j, ridge in enumerate(vor.ridge_vertices):
		# 		if -1 not in ridge:
		# 			line = LineString([ \
		# 				(vertex[ridge[0]][0] + minx, vertex[ridge[0]][1] + miny), \
		# 				(vertex[ridge[1]][0] + minx, vertex[ridge[1]][1] + miny)])
		#
		# 			if line.within(self.inputGEOM) and len(line.coords[0]) > 1:
		# 				lst_lines.append(line)
		#
		# 	return unary_union(lst_lines)


	def densifyBorder(self, polygon, minx, miny):
		"""
		Densifies the border of a polygon by a given factor (by default: 0.5).

		The function tests the complexity of the polygons geometry, i.e. does
		the polygon have holes or not. If the polygon doesn't have any holes,
		its exterior is extracted and densified by a given factor. If the
		polygon has holes, the boundary of each hole as well as its exterior is
		extracted and densified by a given factor.

		Returns:
			a list of points where each point is represented by a list of its
			reduced coordinates.

		Example:
			[[X1, Y1], [X2, Y2], ..., [Xn, Yn]
		"""

		if len(polygon.interiors) == 0:
			exterIN = LineString(polygon.exterior)
			points = self.fixedInterpolation(exterIN, minx, miny)

		else:
			exterIN = LineString(polygon.exterior)
			points = self.fixedInterpolation(exterIN, minx, miny)

			for j in range(len(polygon.interiors)):
				interIN = LineString(polygon.interiors[j])
				points += self.fixedInterpolation(interIN, minx, miny)

		return points


	def fixedInterpolation(self, line, minx, miny):
		"""
		A helping function which is used in densifying the border of a polygon.

		It places points on the border at the specified distance. By default the
		distance is 0.5 (meters) which means that the first point will be placed
		0.5 m from the starting point, the second point will be placed at the
		distance of 1.0 m from the first point, etc. Naturally, the loop breaks
		when the summarized distance exceeds the length of the line.

		Returns:
			a list of points where each point is represented by a list of its
			reduced coordinates.

		Example:
			[[X1, Y1], [X2, Y2], ..., [Xn, Yn]
		"""

		count = self.dist
		newline = []

		startpoint = [line.xy[0][0] - minx, line.xy[1][0] - miny]
		endpoint = [line.xy[0][-1] - minx, line.xy[1][-1] - miny]
		newline.append(startpoint)

		while count < line.length:
			point = line.interpolate(count)
			newline.append([point.x - minx, point.y - miny])
			count += self.dist

		newline.append(endpoint)

		return newline