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

Initial Condition

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


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: enter image description here

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.


  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.

enter image description here

  • Which is the average distance between expected lines? Are points representing plants? – xunilk Mar 15 at 12:40
  • Yes. It's around 7.5 however that're some deviation between lines (as green circle/ line 13, 15, etc). And there might be overlapping in x/y if we try to bin them. – dan Mar 15 at 13:04
  • 7.5 meters or feet? – xunilk Mar 15 at 19:00

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: enter image description here enter image description here

  • Thanks. However, there might be deviation in points location. So then thresholding to create polygon to indicate the boundaries might be difficult.. Also there might be lost points, so if we solely depends on voronoi, the polygons might be unconnected.. – dan Mar 15 at 6:52

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:

enter image description here

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)


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,

prov = mem_layer.dataProvider()

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

for i, feat in enumerate(feats):



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

enter image description here

  • Could you please describe the algorithm? – Comrade Che Mar 16 at 5:57
  • Thanks a lot for your answer. So you use 5 as the distance threshold? I think this also could do the job. But the problem I found that the there might be overlapping in x values between columns. The column rotates and eventually will have same start x values(line) with the end x values at the neighboring lines. Therefore its hard to only depends on delta distance.. cmiiw – dan Mar 16 at 6:20
  • I used 5 as threshold to guarantee an expected result (I analyzed delta list previous to select this value) but it could also be chosen another value (i.e. 6). For small areas algorithm is robust, it was a good election and works well. The problem could be the variability of point separation for bigger areas (for selecting a good threshold). However, in my first image it can be observed a clear hexagonal pattern (except for areas without points ). – xunilk Mar 16 at 12:37
  • If the column rotates and eventually will have same start x values(line) with the end x values at the neighboring lines (for bigger areas), you can compensate this issue rotating complete area before running the algorithm and restoring it at the ending. – xunilk Mar 16 at 12:38
  • @Comrade Che Algorithm is based in distribution of delta values for selecting a threshold and slicing point features based in this criterion. – xunilk Mar 16 at 13:23

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.


  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


enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.