spatialanalysis/kfunction.py at master · gisalgs/spatialanalysis · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
import sys
sys.path.append('../geom')
sys.path.append('../indexing')

from geom.point import *
from indexing.extent import *
from indexing.kdtree1 import *
from indexing.kdtree2b import *

import random
from math import sqrt, pi
import numpy as np

def kfunc(tree, p, d, density):
    """
    Input
      tree:    a k-D tree
      p:       a point where the K-function is computed
      d:       radius of the circle around p
      density: density of points in the area

    Return
      n: count of points in the circle
      ld: L(d) value
    """
    neighbors = []
    range_query_circular(tree, p, d, neighbors)
    n = len(neighbors)
    kd = float(n)/density
    ld = sqrt(kd/pi)
    return n, ld

def kfunc_monte_carlo(n, area, radii, density, rounds=100):
    """
    Input
      n:            number of points
      area:         Extent object defining the area
      radii:        list containing a set of radii of circles
      density:      density of point events in the area
      rounds:       number of simulations
    Return
      percentiles:  a list of 2.5th and 97.5th percentiles
                    for each d in radii
    """
    alllds = []
    for test in range(rounds):
        points = [ Point(random.uniform(area.xmin, area.xmax),
                         random.uniform(area.ymin, area.ymax))
                   for i in range(n) ]
        t = kdtree2(points)
        lds = [0 for d in radii]
        for i, d in enumerate(radii):
            for p in points:
                ld = kfunc(t, p, d, density)[1]
                lds[i] += ld
        lds = [ld/n for ld in lds]
        alllds.append(lds)
    alllds = np.array(alllds)
    percentiles = []
    for i in range(len(radii)):
        percentiles.append([np.percentile(alllds[:,i], 2.5),
                            np.percentile(alllds[:,i], 97.5)])
    return percentiles

def test(points, area):
    """
    Input
      points: a list of Point objects
      area: an Extent object

    Return
      ds: list of radii
      percentiles: Monte Carlo result of using radii in ds
      lds: L(d) values for each radius in ds
    """
    n = len(points)
    density = n/area.area()
    t = kdtree2(points)
    d = min([area.xmax-area.xmin,area.ymax-area.ymin])*2/3/10
    ds = [ d*(i+1) for i in range(10)]
    lds = [0 for d in ds]
    for i, d in enumerate(ds):
        for p in points:
            ld = kfunc(t, p, d, density)[1]
            lds[i] += ld
    lds = [ld/n for ld in lds]
    percentiles = kfunc_monte_carlo(n, area, ds, density)
    return ds, percentiles, lds

if __name__ == "__main__":
    n = 100
    # random patter
    points = [ Point(random.uniform(20, 30),
                     random.uniform(30, 40))
               for i in range(n) ]
    xmin = min([p.x for p in points])
    ymin = min([p.y for p in points])
    xmax = max([p.x for p in points])
    ymax = max([p.y for p in points])
    area = Extent(xmin=xmin, xmax=xmax,
                  ymin=ymin, ymax=ymax)
    ds, percentiles, lds = test(points, area)
    print("Random")
    for i, v in enumerate(percentiles):
        print(ds[i], v[0], lds[i], v[1])

    # three blocks of points
    points1 = [ Point(random.uniform(10, 20),
                      random.uniform(10, 20))
                for i in range(n/3) ]
    points2 = [ Point(random.uniform(30, 40),
                      random.uniform(10, 20))
                for i in range(n/3) ]
    points3 = [ Point(random.uniform(20, 30),
                      random.uniform(30, 40))
                for i in range(n/3) ]
    points = points1 + points2 + points3
    xmin = min([p.x for p in points])
    ymin = min([p.y for p in points])
    xmax = max([p.x for p in points])
    ymax = max([p.y for p in points])
    area = Extent(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)
    ds, percentiles, lds = test(points, area)
    print("Three blocks")
    for i, v in enumerate(percentiles):
        print(ds[i], v[0], lds[i], v[1])