Looking for Python equivalent of QGIS "Nearest Neighbour analysis" tool

I've been using QGIS "Nearest Neighbour analysis" tool to find out the following information about my point shapefile

• Observed mean distance
• Expected mean distance
• Nearest neighbour index
• Number of points
• Z-Score

To do this, I've used shapefile as the input layer and used the advanced settings to limit the features processed to 50 (to save on running time).

This tool gives me exactly what I need but I'm wondering is there an equivalent for Python?

I tried having a nosey at GDAL/OGR in the hope that it had something similar but no luck, unfortunately.

edit An example of my attribute data

This attempt is not finished yet, because it requires Ellipsoidal distance instead of Cartesian.

As you can see in the source code

spatialIndex = QgsSpatialIndex(source, feedback)

distance = QgsDistanceArea()
distance.setSourceCrs(source.sourceCrs(), context.transformContext())
distance.setEllipsoid(context.project().ellipsoid())

There is also no spatial index implemented in my solution.

I hope with my answer I do not collide with the author copyrights of this particular tool.

First of all, I will credit here: Victor Olaya and his code QGIS/python/plugins/processing/algs/qgis/NearestNeighbourAnalysis.py that can be partially reproduced.

And secondly partial credits to @rdmolony with this answer.

Let's assume there is a shapefile called 'points' with 10 point features in it, see the image below.

After applying the "Nearest neighbour analysis" I could get the following output:

Using the following code in Python:

import math
import shapely.geometry
import geopandas as gpd
from shapely.ops import nearest_points

absolute_path_to_shapefile = 'P:/Test/qgis_test/points.shp'

count = len(gdf)
total = 100.0 / count if count else 1

bbox = gdf.total_bounds
polygon = shapely.geometry.box(*bbox, ccw=True)
area = polygon.area

def get_nearest_values(row, other_gdf, point_column='geometry', value_column="geometry"):
"""
Find the nearest point and return the corresponding value from specified value column.
"""

# Create an union of the other GeoDataFrame's geometries:
other_points = other_gdf["geometry"].unary_union
other_points = other_points.difference(row[point_column])

# Find the nearest points
nearest_geoms = nearest_points(row[point_column], other_points)

# Get corresponding values from the other df
nearest_data = other_gdf.loc[other_gdf["geometry"] == nearest_geoms[1]]

nearest_value = nearest_data[value_column].values[0]

return nearest_value

gdf['Nearest'] = gdf.apply(lambda row: get_nearest_values(row, gdf), axis=1)
gdf['Distance'] = gdf.apply(lambda row: row.geometry.distance(row['Nearest']), axis=1)

sumDist = gdf['Distance'].sum()

do = float(sumDist) / count
de = float(0.5 / math.sqrt(count / area))
d = float(do / de)
SE = float(0.26136 / math.sqrt(count ** 2 / area))
zscore = float((do - de) / SE)

print(f'Observed mean distance: {do}')
print(f'Expected mean distance: {de}')
print(f'Nearest neighbour index: {d}')
print(f'Number of points: {count}')
print(f'Z-Score: {zscore}')

I could get the result like this :

Observed mean distance: 496517.1068282208
Expected mean distance: 302718.07444028446
Nearest neighbour index: 1.6401964360610586
Number of points: 10
Z-Score: 3.8729700181269147

References:

• This is brilliant - thank you!! My only question is, is there a way to limit the number of points read in from the shapefile? Commented Apr 29, 2022 at 11:36
• One other issue I've noticed is that my shapefiles all return this error IndexError: index 0 is out of bounds for axis 0 with size 0 Their columns are X, Y , Z so perhaps this is the issue Commented Apr 29, 2022 at 11:58
• I've edited my original question to include a screenshot of my data structure Commented Apr 29, 2022 at 12:07
• Yes its a shapefile. Of course I can send you data - how would you like me to send it? Commented Apr 29, 2022 at 12:59

You can call the same tool directly from python

arguments = {
"INPUT": layer,
"OUTPUT_HTML_FILE": "path/to/html.html"
}
output = processing.run("qgis:nearestneighbouranalysis", arguments) # Qgis 3.16
output = processing.run("native:nearestneighbouranalysis", arguments) # Qgis 3.22

#  output = {
#      "OUTPUT_HTML_FILE" -> [html] # HTML file with the computed statistics
#      "OBSERVED_MD" -> [number] # Observed mean distance
#      "EXPECTED_MD" -> [number] # Expected mean distance
#      "NN_INDEX" -> [number] # Nearest neighbour index
#      "POINT_COUNT" -> [number] #Number of points
#      "Z_SCORE" -> [number] # Z-Score
#  }

nearestneighbouranalysis documentation

see this to know how to use the processing algorithms to run outside qgis

btw the algorithm is really fast even with a 100 000 points it only takes like 5 seconds, i don't know why you would need to limit your analysis to 50 points

• unfortunately, due to me having many versions of python I can't get qgis.core to work in my python script. This is why I was looking an equivalent function in another module. I also limit the analysis because I am sometimes running 35 million points which takes over an hour Commented Apr 29, 2022 at 9:21