Distance matrix returns points in wrong order QGIS 3.4

Question

I'm using QGIS 3.4 on Windows 10 and scripting in pyqgis. I have a layer representing building centroids that is in a projected CRS. I want to find the azimuth between each building and its nearest neighbour in the same layer, so the input and target points being in the correct order is important. However, when I run the following script the closest pairs are output each time with the points in the same order rather than reversed. This renders the output useless. I have tried changing the parameters for MATRIX_TYPE but have kept the other parameters as shown. The data is a vector point layer in projected format.

import os
import ogr

from qgis.core import (
    QgsVectorLayer,
    QgsProject,
    QgsProperty
)
import processing

QgsProject.instance().removeAllMapLayers() # remove all of the layers from the active map

Building_centroids = QgsVectorLayer('C:/Users/ADMIN/Desktop/Mod04AssessedEx/data/bldg_centroids.shp', "Building_centroids", "ogr") # create vector layer, ogr (or something else) needed

if not Building_centroids.isValid():
    print("Layer failed to load!")

QgsProject.instance().addMapLayer(Building_centroids) # add layer to map

distance_matrix = processing.run("qgis:distancematrix",

{
'INPUT': Building_centroids,
'INPUT_FIELD': 'OBJECTID',
'TARGET': Building_centroids,
'TARGET_FIELD': 'OBJECTID',
'MATRIX_TYPE': 0,
'NEAREST_POINTS': 1,
'OUTPUT': 'memory:'
}
)['OUTPUT']

QgsProject.instance().addMapLayer(distance_matrix) # adds the results as a temporary layer to the active map

# this calculates the azimuth between each structure and its nearest neighbour

for structure in distance_matrix.getFeatures():

    structure_geometry = structure.geometry()
    print (structure_geometry)
    structure_geometry = structure.geometry().asMultiPoint()# this line unpacks the multipoint geometry into an iterable or indexable list

    print(structure_geometry[0])
    print(structure_geometry[1])


    azimuth_2_neighbour = structure_geometry[0].azimuth(structure_geometry[1]) # is this working?? no buildings are due west of their nearest neighbour, or is that coincidence


    print(azimuth_2_neighbour)

Answer 1

Following script does the job:

import numpy as np

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]

n = len(feats)

distances = [ [] for i in range(n) ]
indices = [ [] for i in range(n) ]

##determining distances and indices for each pair of points
##itertools python module (for avoiding repeated distances values) produces wrong results
for i in range(n):
    for j in range(n):
        if i != j:
            distances[i].append(feats[i].geometry().distance(feats[j].geometry()))
            indices[i].append([i,j])

min_distance = []

#determining min distance for each pair of points
for i, item in enumerate(distances):
    min_distance.append((np.min(item)))

index_dist = []

#determining index for min distance in distances list
for i, item in enumerate(distances):
    for j, element in enumerate(item):
        if min_distance[i] == element:
            index_dist.append(j)

index_closest_pairs = []

#determining indices for min distance (closest pairs) in indices list
for i, index in enumerate(indices):
    index_closest_pairs.append(index[index_dist[i]])

#determining azimuth for closest_pairs
for item in index_closest_pairs:
    print(feats[item[0]].geometry().asPoint().azimuth(feats[item[1]].geometry().asPoint()))

It determines:

• distances, indices and min distance for each pair of points

• index for min distance in distances list

• indices for min distance (closest pairs) in indices list

• azimuth for closest_pairs

I tested above script with points layer (only 20 points) of following image:

Each closest pair was manually corroborated with my index_closest_pairs list:

[[0, 18], [1, 18], [2, 5], [3, 11], [4, 13], [5, 14], [6, 9], [7, 11], [8, 17], [9, 6], [10, 16], [11, 7], [12, 8], [13, 4], [14, 5], [15, 13], [16, 10], [17, 8], [18, 0], [19, 9]]

and obtained pairs were as expected.

After running above script, azimuths values printed in Python console, in the correct order, were:

-15.17811103645111
151.15349295781988
52.82637775713254
158.15282481491224
14.954150539186188
-12.955243712564535
17.40354742643329
74.83217815029472
98.89967179370501
-162.59645257356672
106.2360708587406
-105.16782184970529
172.15558560326318
-165.04584946081383
167.04475628743546
-155.89332967752063
-73.76392914125941
-81.100328206295
164.8218889635489
-178.84549335216983