How to measure mean distance in search radius by PyQGIS?

Question

I have a point layer with 250 features. I want to find mean distance for each point feature base on 500 search radius. In QGIS by "Distance Matrix" I can measure mean distance for each feature based on all points not in specific search radius. Also with below code (Measure distances between specific features in QGIS for Desktop) I can measure mean distance like "Distance Matrix":

from math import sqrt
layer = iface.activeLayer()
features = layer.getFeatures()
points = []
for feature in features:
    geom = feature.geometry().asPoint()
    points.append(geom)

n = len(points)

for i in range(n-1):
    for j in range(n):
        if i < j:
            print i, j, sqrt(points[i].sqrDist(points[j]))

I tried to improved this code like below to measure mean distance for each point in 500 meter search radius, but the results is not correct.

from math import sqrt
layer = iface.activeLayer()
features = layer.getFeatures()

points = []

for feature in features:
    geom = feature.geometry().asPoint()
    points.append(geom)

n = len(points)
i_count = 0
sum = 0

for i in range(n-1):
    for j in range(n):
        if i < j and (sqrt(points[i].sqrDist(points[j]))) <= 500:
                print i, j, sqrt(points[i].sqrDist(points[j]))
                sum += (sqrt(points[i].sqrDist(points[j])))
                i_count += 1
                print sum / i_count

what is the correct code for this purpose?

Answer 1

The more suitable code for this purpose uses itertools python module for avoiding unnecessary repetitions in distance determinations. Distances that match the requirement (< 500) are grouped in a list for posterior mean determinations. Another list of indexes was used for verification purposes. Next code:

import itertools
import numpy as np

layer = iface.activeLayer()

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

n = len(feats)

comb = range(n)

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

for i, j in itertools.combinations(comb, 2):

    dist = feats[i].geometry().distance(feats[j].geometry())

    if dist < 500:
        distances[i].append(dist)
        indexes[i].append([i,j])

for i, group in enumerate(distances):
    print i, np.mean(group)

print "Done!"

was ran with a point layer with 250 features (see next image).

Results (feature ids and means) were printed at Python Console of QGIS.

0 247.465031352
1 214.562246297
2 nan
3 189.574942919
4 nan
5 217.824344254
6 349.317512486
7 488.048322459
8 361.248613788
9 335.670536721
10 334.301813118
11 331.55702782
12 330.198383101
13 271.797320613
14 351.092830155
15 332.478092381
16 325.4604826
17 323.834616308
18 325.408633423
19 380.609141447
20 327.041047149
21 337.774583784
22 315.274001931
23 305.319935776
24 404.008395136
25 257.484342471
26 326.91975402
27 343.252715434
28 360.400255971
29 322.254034852
30 364.237642274
31 336.930010805
32 356.044460947
33 290.990703848
34 336.36440097
35 476.761971931
36 339.348599874
37 324.07930633
38 370.656980684
39 453.636320183
40 351.254056004
41 474.157655599
42 375.314164631
43 306.364583952
44 299.918081058
45 253.009511802
46 422.482973393
47 469.133917394
48 143.879722927
49 261.979034615
50 181.044111762
51 287.663365184
52 364.508896568
53 231.152521284
54 342.356241239
55 346.586923474
56 343.369152984
57 360.318508131
58 321.78302739
59 nan
60 239.5655864
61 403.71209171
62 364.982715484
63 269.104504035
64 358.271674983
65 371.089942145
66 290.933381774
67 375.100768156
68 441.809464951
69 155.61391372
70 341.070215657
71 440.657626204
72 405.936685729
73 367.491993624
74 287.263625555
75 375.537272425
76 414.110678988
77 267.274594508
78 329.779316558
79 374.234301314
80 nan
81 381.912946129
82 321.845082221
83 404.920663696
84 453.031248814
85 414.331177132
86 293.171318961
87 359.987866123
88 257.159529418
89 343.74939263
90 341.614176993
91 nan
92 313.158111625
93 266.907919749
94 495.216593186
95 403.881726084
96 279.465232144
97 nan
98 339.471903516
99 318.686307553
100 305.156507879
101 450.078167027
102 382.9545129
103 210.968025246
104 307.609312287
105 184.857476899
106 309.366078641
107 307.215074131
108 313.555574773
109 268.748215233
110 347.02004555
111 393.311819568
112 449.550350115
113 334.66447438
114 319.19916216
115 267.182460682
116 378.806724629
117 459.364915029
118 410.790492576
119 447.533507176
120 353.025878937
121 403.301592415
122 387.367121968
123 446.284169266
124 263.045276441
125 316.414394214
126 353.205451115
127 440.107369574
128 470.2322466
129 381.64970501
130 nan
131 370.698426183
132 200.262534215
133 303.279301551
134 315.492559589
135 415.244794115
136 287.087218913
137 375.986761985
138 320.549174839
139 362.305570542
140 98.6390928652
141 216.106459279
142 397.804472302
143 306.132513161
144 282.667737399
145 380.839647486
146 260.326774002
147 340.852588565
148 219.667837193
149 345.326797612
150 423.291171307
151 491.200500139
152 299.155007831
153 236.904296491
154 360.576559435
155 413.515994347
156 256.589776211
157 322.131569393
158 nan
159 nan
160 307.030229451
161 nan
162 308.600730648
163 nan
164 175.918363135
165 nan
166 406.647761353
167 456.157174422
168 356.621227508
169 nan
170 480.764456018
171 440.409967461
172 279.265597101
173 nan
174 202.810210209
175 343.736204062
176 229.055033735
177 nan
178 nan
179 nan
180 379.365568897
181 283.857290617
182 211.310850558
183 nan
184 314.685866836
185 230.257416602
186 nan
187 450.198703463
188 236.528643485
189 334.491501456
190 nan
191 368.026809326
192 nan
193 356.526369035
194 104.567172266
195 399.149296168
196 nan
197 nan
198 313.222681772
199 410.57608185
200 nan
201 453.074278044
202 363.40960241
203 304.891753463
204 183.635649004
205 nan
206 490.466979717
207 nan
208 280.437793607
209 319.068302497
210 nan
211 137.489027251
212 nan
213 341.566332757
214 433.546874239
215 nan
216 nan
217 466.526203964
218 282.28621432
219 174.092213979
220 342.816797481
221 472.338783981
222 nan
223 nan
224 nan
225 412.331021077
226 nan
227 nan
228 nan
229 397.254280768
230 nan
231 nan
232 nan
233 nan
234 nan
235 nan
236 nan
237 331.414378912
238 nan
239 nan
240 nan
241 nan
242 nan
243 nan
244 nan
245 nan
246 nan
247 nan
248 nan
249 nan

For example, above list indicates that for point with id equal 2 there is not mean. To corroborate this, at next image, it can be observed that there is not exist any point at a distance less than 500 m.

Another corroboration was realized with first point and distance tool and it was as expected.

>>>distances[0]
[130.83078719752868, 286.9391760435126, 181.1325815019598, 486.55852003895126, 116.50750219654094, 350.03091547175427, 152.08651552382634, 180.24248462812838, 342.85679956243325]
>>>indexes[0]
[[0, 16], [0, 45], [0, 63], [0, 78], [0, 100], [0, 132], [0, 146], [0, 206], [0, 216]] 

Answer 2

Above answer of xunilk done correctly but assume you have 2 points (A and B) that are in 500 meter search radius together, so based on this code only for first point (i.e. A) mean distance measured and printed but for second point B mean distance reported as 'nan', while, point B have same mean distance that reported for point A. So I improved this code like below to report for two (or more) points.

import itertools
import numpy as np

layer = iface.activeLayer()

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

n = len(feats)

comb = range(n)

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

for i, j in itertools.combinations(comb, 2):

    dist = feats[i].geometry().distance(feats[j].geometry())

    if dist < 500:
        i_dist= distances[i].append(dist)
        i_index= indexes[i].append([i,j])
        j_dist= distances[j].append(dist)
        j_index= indexes[j].append([i,j])

for i, group in enumerate(distances):
    print i, np.mean(group)