I have written a script based on the standard NNA tool which takes point and polygon features as input and performs NNA for each of the polygons using the points within. It also uses polygon area as study area size for the NNA.
The script was for a one-off use, so there’s no real error-checking etc.
import os
import math
import codecs
from PyQt5.QtCore import QVariant
from qgis.PyQt.QtGui import QIcon
from qgis.core import (QgsFeatureRequest,
QgsDistanceArea,
QgsProject,
QgsField,
QgsFields,
QgsFeatureSink,
QgsFeature,
QgsProcessing,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterFeatureSink,
QgsProcessingParameterNumber,
QgsProcessingParameterFileDestination,
QgsProcessingOutputHtml,
QgsProcessingOutputNumber,
QgsSpatialIndex,
QgsWkbTypes)
from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0]
class NearestNeighbourAnalysisWithPoly(QgisAlgorithm):
INPUT = 'INPUT'
AREA = 'AREA'
OUTPUT = 'OUTPUT'
def icon(self):
return QIcon(os.path.join(pluginPath, 'images', 'ftools', 'neighbour.png'))
def group(self):
return self.tr('Vector analysis')
def groupId(self):
return 'vectoranalysis'
def __init__(self):
super().__init__()
def initAlgorithm(self, config=None):
self.addParameter(QgsProcessingParameterFeatureSource(self.INPUT,
self.tr('Points'), [QgsProcessing.TypeVectorPoint]))
self.addParameter(QgsProcessingParameterFeatureSource(self.AREA,
self.tr('Analysis areas'), [QgsProcessing.TypeVectorPolygon]))
self.addParameter(QgsProcessingParameterFeatureSink(self.OUTPUT,
self.tr('Nearest Neighbour Analysis')))
def name(self):
return 'nearestneighbouranalysispolygons'
def displayName(self):
return self.tr('Nearest neighbour analysis within polygons')
def processAlgorithm(self, parameters, context, feedback):
source = self.parameterAsSource(parameters, self.INPUT, context)
area_poly = self.parameterAsSource(parameters, self.AREA, context)
fields = QgsFields()
fields.append(QgsField('fid', QVariant.LongLong))
fields.append(QgsField('Area', QVariant.Double))
fields.append(QgsField('D_obs', QVariant.Double))
fields.append(QgsField('D_exp', QVariant.Double))
fields.append(QgsField('N_points', QVariant.Int))
fields.append(QgsField('NNI', QVariant.Double))
fields.append(QgsField('Z_score', QVariant.Double))
(sink, dest_id) = self.parameterAsSink(
parameters, self.OUTPUT, context,
fields, QgsWkbTypes.MultiPolygon, area_poly.sourceCrs())
distance = QgsDistanceArea()
distance.setSourceCrs(source.sourceCrs(), context.transformContext())
distance.setEllipsoid(context.project().ellipsoid())
measure = QgsDistanceArea()
measure.setSourceCrs(source.sourceCrs(), context.transformContext())
measure.setEllipsoid(context.project().ellipsoid())
all_features = list(source.getFeatures())
n_poly = area_poly.featureCount()
total = 100.0 / n_poly if n_poly else 1
for current, poly in enumerate(area_poly.getFeatures()):
if feedback.isCanceled():
break
fid = poly['fid']
feedback.pushInfo(f'{current=}: {fid=}')
A = measure.measureArea(poly.geometry())
features = []
for f in all_features:
if f.geometry().intersects(poly.geometry()):
features.append(f)
count = len(features)
if count > 1:
spatialIndex = QgsSpatialIndex()
spatialIndex.addFeatures(features)
sumDist = 0.00
for feat in features:
neighbourID = spatialIndex.nearestNeighbor(
feat.geometry().asPoint(), 2)[1]
request = QgsFeatureRequest().setFilterFid(neighbourID).setSubsetOfAttributes([])
neighbour = next(source.getFeatures(request))
sumDist += distance.measureLine(neighbour.geometry().asPoint(),
feat.geometry().asPoint())
do = float(sumDist) / count
de = float(0.5 / math.sqrt(count / A))
nni = float(do / de)
SE = float(0.26136 / math.sqrt(count ** 2 / A))
zscore = float((do - de) / SE)
else:
do = None
de = None
nni = None
zscore = None
new_feature = QgsFeature()
new_feature.setGeometry(poly.geometry())
new_feature.setAttributes([fid, A, do, de, count, nni, zscore])
sink.addFeature(new_feature, QgsFeatureSink.FastInsert)
feedback.setProgress(int(current * total))
return {self.OUTPUT: dest_id}