I have address points connected to a network of roads that I have built using GRASS v.net.connect. I would like to cluster addresses into fixed size groups along the network.
So far, I have written an algorithm that:
- Uses the southern most address point as the "start point" to begin.
- Finds the shortest path between the start point and all other address points using the
native:shortestpathpointtopointalgorithm - Associates these two points (the start point and nearest end point) to a cluster
- Set the previous nearest end point as the next start point
- Find the nearest address point using the shortest path from this start point to get the next address point to add to the cluster
- repeat this process until all address points have been clustered incrementing the cluster number every N points
I could not tell you if this method works because the algorithm I have written is estimated to take about 55 hours using 416 address points. I have included the code below and would like any suggestions you might have on how I can improve its speed. Or if there is a better solution to doing what I am trying to achieve.
def get_layer_by_name_quick(namepattern):
layers = []
for lyr in QgsProject.instance().mapLayers().values():
if namepattern in str(lyr.name()):
layers.append(lyr)
if len(layers) == 0:
print(f"No Layer named : {namepattern}")
return None
elif len(layers) > 1:
print(f"multiple layers returned: \nnamepattern: {namepattern}, \nlayers: {layers}\n")
return layers[0]
else:
return layers[0]
class ClusterAddressesV2(object):
def __init__(self, addresses, network, clustersize = 30):
self.addresses = addresses
self.network = network
self.clustersize = clustersize
self.cluster = 0
self.clusters = {}
self.collected_points = []
self.startPoint, self.start_point = self.get_min_y_point()
self.all_points = {feat.id() : feat.geometry().asPoint() for feat in self.addresses.getFeatures()}
def get_min_y_point(self):
pnts = [(addr.id(), addr.geometry().asPoint()) for addr in self.addresses.getFeatures() if addr.id() not in self.collected_points]
xys = [(pnt[1].x(), pnt[1].y()) for pnt in pnts]
xysorted = sorted(xys, key = lambda item: item[1]) # sort by y
first_point = xysorted[0]
pnt0 = [(p[0], p[1]) for p in pnts if p[1].y() == first_point[1]] [0]
pid, point = pnt0
return pid, point
def get_all_points(self):
the_points = {feat.id() : feat.geometry().asPoint() for feat in self.addresses.getFeatures() if feat.id() not in self.collected_points}
return the_points
def get_shortest_paths(self):
the_points = self.get_all_points() # get all points
shortest_paths = []
for k, v in the_points.items():
if v != self.start_point:
end_point = v
endPoint = k
else:
# k is startPoint
self.startPoint = k
endPoint = k
end_point = v
if self.startPoint not in self.collected_points:
self.collected_points.append(self.startPoint)
feedback = QgsProcessingFeedback()
parameters = {'INPUT': self.network,
'STRATEGY': 0,
'DIRECTION_FIELD': '',
'VALUE_FORWARD': '',
'VALUE_BACKWARD': '',
'VALUE_BOTH': '',
'DEFAULT_DIRECTION': 2,
'SPEED_FIELD': '',
'DEFAULT_SPEED': 1,
'TOLERANCE': 15, #tolerance
'START_POINT': self.start_point,
'END_POINT': end_point,
'OUTPUT': 'memory:shortestpath'}
try:
branch = processing.run('native:shortestpathpointtopoint', parameters, feedback = feedback)
shortest_paths.append((self.startPoint, endPoint, branch['TRAVEL_COST']))
except Exception as e:
print(f'Exception: {e} startPoint : {self.startPoint}, endPoint: {endPoint}')
# after running through all the points
s = sorted( shortest_paths, key = lambda item: item[2])
s1 = s[0][1] #endPoint
s0 = s[0][0] #startPoint
if s1 not in list(self.clusters.keys()):
self.clusters[s1] = self.cluster
if s0 not in list(self.clusters.keys()):
self.clusters[s0] = self.cluster
if s1 not in self.collected_points:
self.collected_points.append(s1)
if s0 not in self.collected_points:
self.collected_points.append(s0)
self.startPoint = s1
self.start_point = self.all_points[self.startPoint]
if len(list(self.clusters.keys())) % self.clustersize == 0:
self.cluster += 1
return
def write_attributes(self):
all_features = [f for f in self.addresses.getFeatures()]
layer_provider = self.addresses.dataProvider()
self.addresses.startEditing()
fields = self.addresses.fields()
clusterIdx = fields.indexFromName('Cluster')
for f in all_features:
id = f.id()
if id in list(self.clusters.keys()):
cluster = self.clusters[id]
#print(cluster)
attr_value={clusterIdx : cluster}
layer_provider.changeAttributeValues({id:attr_value})
self.addresses.commitChanges()
return
def run_algo(self):
all_features = [f for f in self.addresses.getFeatures()]
while len(list(self.clusters.keys())) < len(all_features):
self.get_shortest_paths()
if len(list(self.clusters.keys())) >= len(all_features):
break
self.write_attributes()
addresses = get_layer_by_name_quick('repr_addresses')
network = get_layer_by_name_quick('Network')
cav = ClusterAddressesV2(addresses, network, 25)
cav.run_algo()
