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

Question

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

Answer 1

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'

gdf = gpd.read_file(absolute_path_to_shapefile)

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:

• pandas python how to count the number of records or rows in a dataframe
• Calculate Distance to Nearest Feature with Geopandas
• Finding nearest point in other GeoDataFrame using GeoPandas
• Getting polygon breadth in Shapely
• Calculating bounding box size?
• Filter by bounding box in geopandas?
• GeoPandas Set CRS on Points

Answer 2

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