Recursion error in catchment polygons

[rtsaboya]

Hi!

Main issue:

I am having trouble generating catchment polygons for points over an OSM network. Here is the error message, and below I provide more information about my analysis and setup:

2021-10-05T13:49:25 CRITICAL Catchment Analyser raised an exception: Traceback (most recent call last): File "C:\Users/renato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\esstoolkit\catchment_analyser\catchment_analysis.py", line 111, in analysis output_polygon_features = self.polygon_writer( File "C:\Users/renato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\esstoolkit\catchment_analyser\catchment_analysis.py", line 417, in polygon_writer hull = self.concave_hull.concave_hull(points, polygon_tolerance) File "C:\Users/renato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\esstoolkit\catchment_analyser\analysis_tools.py", line 457, in concave_hull return self.concave_hull(points_list, kk + 1) File "C:\Users/renato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\esstoolkit\catchment_analyser\analysis_tools.py", line 457, in concave_hull return self.concave_hull(points_list, kk + 1) File "C:\Users/renato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\esstoolkit\catchment_analyser\analysis_tools.py", line 457, in concave_hull return self.concave_hull(points_list, kk + 1) [Previous line repeated 989 more times] File "C:\Users/renato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\esstoolkit\catchment_analyser\analysis_tools.py", line 406, in concave_hull c_points = self.sort_by_angle(k_nearest_points, current_point, previous_angle) File "C:\Users/renato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\esstoolkit\catchment_analyser\analysis_tools.py", line 340, in sort_by_angle vertex_list = sorted(list_of_points, key=getkey, reverse=True) File "C:\Users/renato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\esstoolkit\catchment_analyser\analysis_tools.py", line 338, in getkey return self.angle_difference(last_angle, self.angle(last_point, item)) File "C:\Users/renato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\esstoolkit\catchment_analyser\analysis_tools.py", line 121, in angle return math.atan2(to_point[1] - from_point[1], to_point[0] - from_point[0]) RecursionError: maximum recursion depth exceeded while calling a Python object

My point layer has 1120 features; my network was downloaded and cleaned in QGIS and further processed using the "Road Network Cleaner" from SST. It has 34,390 segments.

My QGIS version is 3.16.11, the latest LTR available at the time of writing. I'm on Windows 10.

Something very weird is that this problem only occurs when I select a 500m radius. When I try it with 1200m, it works fine.

Possibly related issue

Testing with only one point in the 500m radius, I found that the algorithm went all the way through, but the result is not very accurate:

The distance in the segments look alright, but the resulting polygon does not.

Any help is appreciated. Best, Renato.

[jorgegil]

Hi @rtsaboya

We need to look at this specific instance, but the catchment polygons are a known problem… we know they work most of the time, but can sometimes fail. We’ve adopted a simplified concave hull algorithm that is not bullet proof. I’ll see if I can reproduce and patch the issue.

In essence there is not one real solution for drawing such catchment polygons, and there are many approaches. Here we went for the concave hull solution, that potentially looks the nicest and is how people expect a catchment polygon to look. BUT it’s the most computationally intensive and the algorithm is recursive. In some cases (not related to size but geometry) it hits configurations where it struggles.

In my own work I often just use dissolved buffers on the segments, then remove holes. Doesn’t look as nice but is more correct. For a super fast and simple solution, I might also use convex hulls, but is way less correct. Others do voronoi tesselations and dissolve them to make the polygons…

As you can see there are many approaches, and no one is 100% correct. The only correct solution are the segments themselves, and those the tool produces. From that, any of the above methods can be applied. And QGIS has the tools for it.

I know this doesn’t give you a solution to your problem, but rather points to ways to move forward from the results that you get.

Best, Jorge

[rtsaboya]

Dear @jorgegil

Thank you for the detailed and thoughtful explanation.

I think maybe your suggestion about using buffers may solve my problem. However, I cannot figure out how to produce the buffers while preserving the IDs of the original points. I examined the segments table and it seems that the point ID is not carried through to the segments as it is for the polygons.

Ideally, then, I would have one dissolved polygon for each original point, with the ID of the point in the buffer's table. Is it possible?

[jorgegil]

Dear @rtsaboya

The ID of the original point is the field name of the distances on the segments. That is, if when you calculate catchment you check the box 'Origin Names'. This way you can create a buffer for the segments under a certain value in each distance column.

The min distance you can use if you want dissolved buffers, not per origin.

Jorge