8

I created a grid of hexagons in QGIS using MMQGIS. I have a set of points and I want to match each with a unique hexagon. No hexagon should have more than 1 point assigned. Is there a way that I can systematically move each point to the nearest center of the hexagon (or 2nd or 3rd nearest if the center of that hexagon has already been assigned to a closer point)?

Here is a snapshot of what I am looking at:

enter image description here

2
  • Is this option suitable - first create a point for each hex, and then enter the necessary atiributive information? Dec 22, 2021 at 10:05
  • Is there a unique id number for points in the point layer? Jan 16 at 20:14

2 Answers 2

10

You can try the script below.

Note: Before you do anything else, please make a backup copy of your point layer file in case of unexpected or undesired results (don't just duplicate the layer in the TOC) because this script will permanently modify the point layer's geometry.

Additionally, both layers must be in the same CRS.

As per @bugmenot123's comment you may also like to save the original geometries as x& y columns in your point layer to allow for analysing the displacement later.

Once you have made a backup and checked that both layers are using the same crs, paste the script below into a blank editor in the Python console. Change the layer names in lines 3 & 4 to match your own grid and point layer names. Click Run.

Disclaimer- I wrote this script fairly hastily so it could probably be cleaned up a bit, but I have tested in QGIS 3.20 and it works for me on a couple of test layers.

project = QgsProject.instance()
#Change to match your layer names below
grid = project.mapLayersByName('Hex_grid')[0]
points = project.mapLayersByName('Random_points')[0]

point_feats = [f for f in points.getFeatures()]
grid_feats = [f for f in grid.getFeatures()]

def find_nearest_empty_hexagon(feat):
    index = QgsSpatialIndex(grid.getFeatures())
    n = 0
    found = False
    while not found:
        n+=1
        nn_id = index.nearestNeighbor(feat.geometry(), n)[n-1]
        geom = grid.getFeature(nn_id).geometry()
        if not geom.intersects(QgsGeometry().collectGeometry([f.geometry() for f in points.getFeatures()])):
            found = True
    return nn_id

points_moved = []
for hex_ft in grid_feats:
    points_within = [f for f in point_feats if f.geometry().within(hex_ft.geometry())]
    if len(points_within) > 1:
        for p in points_within:
            current_centroid = hex_ft.geometry().centroid()
            min_dist = min([p.geometry().distance(current_centroid) for p in points_within])
            if p.geometry().distance(current_centroid) == min_dist:
                continue
            nearest_empty_id = find_nearest_empty_hexagon(p)
            nearest_empty_feat = [f for f in grid_feats if f.id() == nearest_empty_id][0]
            points.dataProvider().changeGeometryValues({p.id(): nearest_empty_feat.geometry().centroid()})
            points_moved.append(p.id())

points_outside_grid = []

remaining_points = [p for p in point_feats if p.id() not in points_moved]
for p in remaining_points:
    enclosing_hex = [f for f in grid_feats if f.geometry().contains(p.geometry())]
    if not enclosing_hex:
        points_outside_grid.append(p)
    elif len(enclosing_hex)>0:
        enclosing_hex_centroid = enclosing_hex[0].geometry().centroid()
    points.dataProvider().changeGeometryValues({p.id(): enclosing_hex_centroid})
    
if points_outside_grid:
    for p in points_outside_grid:
        nearest_empty_id = find_nearest_empty_hexagon(p)
        nearest_empty_feat = [f for f in grid_feats if f.id() == nearest_empty_id][0]
        points.dataProvider().changeGeometryValues({p.id(): nearest_empty_feat.geometry().centroid()})
    
points.triggerRepaint()

Screencast demo below:

enter image description here

Before:

enter image description here

After:

enter image description here

6
  • 2
    I'd recommend storing the original locations as WKT or x/y columns so one can analyse the displacement later. Dec 22, 2021 at 8:27
  • @bugmenot123, good idea :-)
    – Ben W
    Dec 22, 2021 at 9:55
  • Hi Ben, thank you for this, it's exactly what I need. However, I am getting this error in the Python Console: "Traceback (most recent call last): File "C:\OSGeo4W\apps\Python39\lib\code.py", line 90, in runcode exec(code, self.locals) File "<input>", line 1, in <module> File "<string>", line 37, in <module> IndexError: list index out of range" Any suggestions? Dec 29, 2021 at 23:13
  • @Jess Siegel, checking line 37 in the code block, my guess is that you have one or more points which do not actually lie within any hexagon feature (maybe they are just outside the hexagon grid). I have edited the code in my answer which should deal with the error (but any points outside the hex grid will not be moved). Can you let me know if my guess is correct- if so, I will try to add some logic to the algorithm to deal with that case.
    – Ben W
    Dec 30, 2021 at 0:59
  • 1
    @Jess Siegel, I have edited my answer again. This version of the script should work also if some points lie outside of hexagon grid features. I hope it works for you.
    – Ben W
    Dec 30, 2021 at 1:52
8

You can use the following script. It creates a new layer for the relocated points.

# h: hexagon, p: point
h_lyr = QgsProject.instance().mapLayersByName("Hexagons")[0]
p_lyr = QgsProject.instance().mapLayersByName("Points")[0]

distances = {p.id(): 
                {h.id(): h.geometry().distance(p.geometry().centroid())
                    for h in h_lyr.getFeatures()}
             for p in p_lyr.getFeatures()}

sort = lambda x: sorted(x.items(), key=lambda item: item[1])
sorted_distances = {k: {m: n for m, n in sort(v)} for k, v in distances.items()}

new_positions = {}
assigned =  [] # assigned a point to the hexagon id?
for pid, dists in sorted_distances.items(): # pid: point id
    iter_ = iter(dists.keys())
    hid = next(iter_) # hid: hexagon id
    
    while hid in assigned:
        hid = next(iter_)
        
    hex_cent = h_lyr.getFeature(hid).geometry().centroid()
    new_positions[pid] = hex_cent
    assigned.append(hid)

# Crate a new layer and add the relocated points
crs = p_lyr.crs().authid()
lyr = QgsVectorLayer("Point?crs=" + crs, "New Locations", "memory")
for new_position in new_positions.values():
    feature = QgsFeature()
    feature.setGeometry(new_position)
    lyr.dataProvider().addFeatures([feature])
    lyr.updateExtents()
    
QgsProject.instance().addMapLayer(lyr)

enter image description here

3
  • Since my English is not so good, explaining how the script works takes me longer than writing the code. Sorry. :( Jan 2 at 15:57
  • Thank you! This worked! The only issue is, I need a way to connect the new locations with the old locations so I know which points were assigned where. Jan 5 at 16:46
  • @JessSiegel I am not sure if it works but add a new field for both layers using Field Calculator, use $id for expression. And use "join by lines (hub lines)" to connect points (use newly created fields as ID fields in the tool window) Jan 16 at 20:12

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.