Setting the study area boundaries in QGIS for nearest neighbor analysis?

Question

I've made a cursory (i.e., two hour) search of the Stack Exchange and the Internet to find out how we can set the study area boundary for the nearest neighbor analysis command in QGIS 3.16. Sadly, I have found no definitive answers.

So here are my questions:

1. Is it possible in QGIS to set the study area boundaries? I have already tried running the NNA operation with my study area boundary on, my study area boundary off, then reprojected to WGS 84 (to make sure it met the lat longitude of the point data layer instead of relying on projections on the fly) and then by setting the scale of the window so that the study area just fit inside the window. In all four cases, I always got the same results. Sadly, there appears to be no option in the Nearest Neighbor Analysis window that let's us set the study area boundary.

2. If it is possible to set the study area boundaries for the NNA, how do we do it?

Answer 1

Using Nearest neighbor analysis, the study area is automatically limited by the extent of the layer you use as input. Thus the x_min/x_max and y_min/y_max value define the limit: northernmost point for North, southernmost point for South, westmost point for West, eastmost point for East.

Help says:

The output describes how the data are distributed (clustered, randomly or distributed).

Output is a list of values: Observed mean distance, Expected mean distance, Nearest neighbour index, Number of points and Z-Score. No need to define any study area.

If you want to do an analysis for just some, but not all points of a layer, make a selection and check the box next to Selected features only (available in many QGIS tools, by the way), see screeshot (yellow points are selected):

Answer 2

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}