Snapping a point to a unique hexagon in QGIS

Question

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:

Answer 1

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:

Before:

After:

Answer 2

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)