Connecting points with parallel lines

Question

Background

I'm trying to create parallel lines that go through each parallel points. In this points arrangement, it's North-South lines:

I use Python and the libraries: pygeos, geopandas, and shapely

Difficulties

My current approach is to create Delaunay lines, then filter it based on its azimuth. I keep the deviation around +- 10 degree from 0 and 180 (0 & 180 based on minimum_rotated_rectagle orientation)

Currently here's my result:

Here I see there're two things that should be done.

1. Connect the lines
2. Connect points that are not in any lines

To solve both of this problem, I'm thinking to create a connecting line for each line/point to nearby lines (filtered by threshold distance (x and y); ex: line 8 has nearby lines -> 6, 7, 9). When connecting the lines the connecting line created must be either at the start or end of the line.

Using the created connecting lines, then I filter it again using the azimuth.

Questions

1. What line should I use, is it hausdorff lines? (if it's, is there library that create the line? not only calculate the distance).
2. Is there any better approach for this problem?

I've tried to use voronoi and buffer to simulate neighboring points, however it's not possible because there could be deviation in points position, therefore not creating a straight line.

Desired Output

Here's the desired output. The lines connected to points in its parallel.

• For point in green circle, it will be combined to neighboring line based on the shortest distance.

Answer 1

Also you can try the following:

1. Create vertical polygons to indicate the boundaries of the future line.
2. Number each point within these boundaries (for example, all red points are 1, all green points are 2, etc.)
3. Number each point sequentially from bottom to top.
4. Connect the points sequentially.

A Voronoi diagram can help you identify a unique line number:

Answer 2

Assuming that units in your comment are in meters, I digitized your image (approximate proportion of 6:4) to a dimension of 165x110 square meters (taking in account 7.5 meters as separation). It looks as follows:

Axis is an arbitrary line put in that location for testing my algorithm based in distances with PyQGIS3. You can adapt it for using with pygeos, geopandas, and shapely python modules. Delta variable is determined for searching where distance differences show greater variation (about 7.5-9.0 meters). It is used for slicing points features based in this criterion. Complete code looks as follows:

from math import sqrt, fabs

registry = QgsProject.instance()

line = registry.mapLayersByName('axis')
points = registry.mapLayersByName('points_to_line')

feat_line = [ feat for feat in line[0].getFeatures() ]

feat_points = [ feat for feat in points[0].getFeatures() ]

points = [ feat.geometry().asPoint() for feat in feat_points ]

distances = [ sqrt(feat_line[0].geometry().closestSegmentWithContext(feat.geometry().asPoint())[0])
              for feat in feat_points ]

delta = [ fabs(distances[i] - distances[i+1]) for i in range(len(distances)-1) ]

cut_index = [ i for i, item in enumerate(delta) if item > 5 ]

cut_index.insert(len(cut_index), len(distances)-1)

cut_index2 = [ item + 1 for item in cut_index ]
cut_index2.insert(0, 0)

del(cut_index2[-1])

slices = [ [cut_index2[i], cut_index[i]] for i in range(len(cut_index)) ]

epsg = line[0].crs().postgisSrid()

uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'line',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(slices)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolylineXY(points[slices[i][0]:slices[i][1]+1]))

prov.addFeatures(feats)

registry.addMapLayer(mem_layer)

After running above code in Python Console of QGIS, I got quickly the result of following image.

Answer 3

So got one approach that could do the job.

Few conditions:

1. The points are equally spaced. Although there might be deviations.
2. The points are in rectangle bounding area.

Solution:

1. Using shapely, create minimum_rotated_rectangle
2. Calculate the rotation of the bounding rectangle
3. Affine transform based on the rotation angle
4. Create delaunay triangles
5. Filter the delaunay lines/edges based on angle & distance
6. Bin the lines based on xbar of the lines (Because points position is rotated, it's posible to bin the points/lines)
7. Join lines in the same bin
8. Create fitted line from points in (7) -> then rotate it to be at 0 deg angle
9. Join unconnected points by smallest deviation to the rotated fitted line

Result:

But I think it's only can be implemented if the points are in straight line arrangement..