1import math
2import scipy
3import numpy
4from .shapes import Rectangle, Point, LineSegment
5from .standalone import get_segment_point_dist, get_bounding_box
6import random
7import time
8import warnings
9
10
11dep_msg = "is deprecated and will be reoved in libpysal 4.4.0."
12
13__all__ = ["SegmentGrid", "SegmentLocator", "Polyline_Shapefile_SegmentLocator"]
14DEBUG = False
15
16
17class BruteSegmentLocator(object):
18    def __init__(self, segments):
19        self.data = segments
20        self.n = len(segments)
21
22    def nearest(self, pt):
23        d = self.data
24        distances = [get_segment_point_dist(d[i], pt)[0] for i in range(self.n)]
25        return numpy.argmin(distances)
26
27
28class SegmentLocator(object):
29    def __init__(self, segments, nbins=500):
30        warnings.warn("SegmentLocator " + dep_msg, DeprecationWarning)
31        self.data = segments
32        if hasattr(segments, "bounding_box"):
33            bbox = segment.bounding_box
34        else:
35            bbox = get_bounding_box(segments)
36        self.bbox = bbox
37        res = max((bbox.right - bbox.left), (bbox.upper - bbox.lower)) / float(nbins)
38        self.grid = SegmentGrid(bbox, res)
39        for i, seg in enumerate(segments):
40            self.grid.add(seg, i)
41
42    def nearest(self, pt):
43        d = self.data
44        possibles = self.grid.nearest(pt)
45        distances = [get_segment_point_dist(d[i], pt)[0] for i in possibles]
46        # print "possibles",possibles
47        # print "distances",distances
48        # print "argmin", numpy.argmin(distances)
49        return possibles[numpy.argmin(distances)]
50
51
52class Polyline_Shapefile_SegmentLocator(object):
53    def __init__(self, shpfile, nbins=500):
54        warnings.warn(
55            "Polyline_Shapefile_SegmentLocator " + dep_msg, DeprecationWarning
56        )
57        self.data = shpfile
58        bbox = Rectangle(*shpfile.bbox)
59        res = max((bbox.right - bbox.left), (bbox.upper - bbox.lower)) / float(nbins)
60        self.grid = SegmentGrid(bbox, res)
61        for i, polyline in enumerate(shpfile):
62            for p, part in enumerate(polyline.segments):
63                for j, seg in enumerate(part):
64                    self.grid.add(seg, (i, p, j))
65
66    def nearest(self, pt):
67        d = self.data
68        possibles = self.grid.nearest(pt)
69        distances = [
70            get_segment_point_dist(d[i].segments[p][j], pt)[0]
71            for (i, p, j) in possibles
72        ]
73        # print "possibles",possibles
74        # print "distances",distances
75        # print "argmin", numpy.argmin(distances)
76        return possibles[numpy.argmin(distances)]
77
78
79class SegmentGrid(object):
80    """
81    Notes:
82        SegmentGrid is a low level Grid class.
83        This class does not maintain a copy of the geometry in the grid.
84        It returns only approx. Solutions.
85        This Grid should be wrapped by a locator.
86    """
87
88    def __init__(self, bounds, resolution):
89        """
90        Returns a grid with specified properties.
91
92        __init__(Rectangle, number) -> SegmentGrid
93
94        Parameters
95        ----------
96        bounds      : the area for the grid to encompass
97        resolution  : the diameter of each bin
98
99        Examples
100        --------
101        TODO: complete this doctest
102        >>> g = SegmentGrid(Rectangle(0, 0, 10, 10), 1)
103        """
104        warnings.warn("SegmentGrid " + dep_msg, DeprecationWarning)
105        if resolution == 0:
106            raise Exception("Cannot create grid with resolution 0")
107        self.res = resolution
108        self.hash = {}
109        self._kd = None
110        self._kd2 = None
111        self._hashKeys = None
112        self.x_range = (bounds.left, bounds.right)
113        self.y_range = (bounds.lower, bounds.upper)
114        try:
115            self.i_range = (
116                int(math.ceil((self.x_range[1] - self.x_range[0]) / self.res)) + 1
117            )
118            self.j_range = (
119                int(math.ceil((self.y_range[1] - self.y_range[0]) / self.res)) + 1
120            )
121            self.mask = numpy.zeros((self.i_range, self.j_range), bool)
122            self.endMask = numpy.zeros((self.i_range, self.j_range), bool)
123        except Exception:
124            raise Exception(
125                "Invalid arguments for SegmentGrid(): ("
126                + str(self.x_range)
127                + ", "
128                + str(self.y_range)
129                + ", "
130                + str(self.res)
131                + ")"
132            )
133
134    @property
135    def hashKeys(self):
136        if self._hashKeys is None:
137            self._hashKeys = numpy.array(list(self.hash.keys()), dtype=float)
138        return self._hashKeys
139
140    @property
141    def kd(self):
142        if self._kd is None:
143            self._kd = scipy.spatial.cKDTree(self.hashKeys)
144        return self._kd
145
146    @property
147    def kd2(self):
148        if self._kd2 is None:
149            self._kd2 = scipy.spatial.KDTree(self.hashKeys)
150        return self._kd2
151
152    def in_grid(self, loc):
153        """
154        Returns whether a 2-tuple location _loc_ lies inside the grid bounds.
155        """
156        return (
157            self.x_range[0] <= loc[0] <= self.x_range[1]
158            and self.y_range[0] <= loc[1] <= self.y_range[1]
159        )
160
161    def _grid_loc(self, loc):
162        i = int((loc[0] - self.x_range[0]) / self.res)  # floored
163        j = int((loc[1] - self.y_range[0]) / self.res)  # floored
164        # i = min(self.i_range-1, max(int((loc[0] - self.x_range[0])/self.res), 0))
165        # j = min(self.j_range-1, max(int((loc[1] - self.y_range[0])/self.res), 0))
166        # print "bin:", loc, " -> ", (i,j)
167        return (i, j)
168
169    def _real_loc(self, grid_loc):
170        x = (grid_loc[0] * self.res) + self.x_range[0]
171        y = (grid_loc[1] * self.res) + self.y_range[0]
172        return x, y
173
174    def bin_loc(self, loc, id):
175        grid_loc = self._grid_loc(loc)
176        if grid_loc not in self.hash:
177            self.hash[grid_loc] = set()
178            self.mask[grid_loc] = True
179        self.hash[grid_loc].add(id)
180        return grid_loc
181
182    def add(self, segment, id):
183        """
184        Adds segment to the grid.
185
186        add(segment, id) -> bool
187
188        Parameters
189        ----------
190        id -- id to be stored int he grid.
191        segment -- the segment which identifies where to store 'id' in the grid.
192
193        Examples
194        --------
195        >>> g = SegmentGrid(Rectangle(0, 0, 10, 10), 1)
196        >>> g.add(LineSegment(Point((0.2, 0.7)), Point((4.2, 8.7))), 0)
197        True
198        """
199        if not (self.in_grid(segment.p1) and self.in_grid(segment.p2)):
200            raise Exception(
201                "Attempt to insert item at location outside grid bounds: "
202                + str(segment)
203            )
204        i, j = self.bin_loc(segment.p1, id)
205        I, J = self.bin_loc(segment.p2, id)
206        self.endMask[i, j] = True
207        self.endMask[I, J] = True
208
209        bbox = segment.bounding_box
210        left = bbox.left
211        lower = bbox.lower
212        res = self.res
213        line = segment.line
214        tiny = res / 1000.0
215        for i in range(1 + min(i, I), max(i, I)):
216            # print 'i',i
217            x = self.x_range[0] + (i * res)
218            y = line.y(x)
219            self.bin_loc((x - tiny, y), id)
220            self.bin_loc((x + tiny, y), id)
221        for j in range(1 + min(j, J), max(j, J)):
222            # print 'j',j
223            y = self.y_range[0] + (j * res)
224            x = line.x(y)
225            self.bin_loc((x, y - tiny), id)
226            self.bin_loc((x, y + tiny), id)
227        self._kd = None
228        self._kd2 = None
229        return True
230
231    def remove(self, segment):
232        self._kd = None
233        self._kd2 = None
234        pass
235
236    def nearest(self, pt):
237        """
238        Return a set of ids.
239
240        The ids identify line segments within a radius of the query point.
241        The true nearest segment is guaranteed to be within the set.
242
243        Filtering possibles is the responsibility of the locator not the grid.
244        This means the Grid doesn't need to keep a reference to the underlying segments,
245        which in turn means the Locator can keep the segments on disk.
246
247        Locators can be customized to different data stores (shape files, SQL, etc.)
248        """
249        grid_loc = numpy.array(self._grid_loc(pt))
250        possibles = set()
251
252        if DEBUG:
253            print("in_grid:", self.in_grid(pt))
254            i = pylab.matshow(
255                self.mask, origin="lower", extent=self.x_range + self.y_range, fignum=1
256            )
257        # Use KD tree to search out the nearest filled bin.
258        # it may be faster to not use kdtree, or at least check grid_loc first
259        # The KD tree is build on the keys of self.hash, a dictionary of stored bins.
260        dist, i = self.kd.query(grid_loc, 1)
261
262        ### Find non-empty bins within a radius of the query point.
263        # Location of Q point
264        row, col = grid_loc
265        # distance to nearest filled cell +2.
266        # +1 returns inconsistent results (compared to BruteSegmentLocator)
267        # +2 seems to do the trick.
268        radius = int(math.ceil(dist)) + 2
269        if radius < 30:
270            a, b = numpy.ogrid[
271                -radius : radius + 1, -radius : radius + 1
272            ]  # build square index arrays centered at 0,0
273            index = (
274                a ** 2 + b ** 2 <= radius ** 2
275            )  # create a boolean mask to filter indicies outside radius
276            a, b = index.nonzero()
277            # grad the (i,j)'s of the elements within radius.
278            rows, cols = (
279                row + a - radius,
280                col + b - radius,
281            )  # recenter the (i,j)'s over the Q point
282            #### Filter indicies by bounds of the grid.
283            ### filters must be applied one at a time
284            ### I havn't figure out a way to group these
285            filter = rows >= 0
286            rows = rows[filter]
287            cols = cols[filter]  # i >= 0
288            filter = rows < self.i_range
289            rows = rows[filter]
290            cols = cols[filter]  # i < i_range
291            filter = cols >= 0
292            rows = rows[filter]
293            cols = cols[filter]  # j >= 0
294            filter = cols < self.j_range
295            rows = rows[filter]
296            cols = cols[filter]  # j < j_range
297            if DEBUG:
298                maskCopy = self.mask.copy().astype(float)
299                maskCopy += self.endMask.astype(float)
300                maskCopy[rows, cols] += 1
301                maskCopy[row, col] += 3
302                i = pylab.matshow(
303                    maskCopy,
304                    origin="lower",
305                    extent=self.x_range + self.y_range,
306                    fignum=1,
307                )
308                # raw_input('pause')
309            ### All that was just setup for this one line...
310            idx = self.mask[rows, cols].nonzero()[0]  # Filter out empty bins.
311            rows, cols = (
312                rows[idx],
313                cols[idx],
314            )  # (i,j)'s of the filled grid cells within radius.
315
316            for t in zip(rows, cols):
317                possibles.update(self.hash[t])
318
319            if DEBUG:
320                print("possibles", possibles)
321        else:
322            ### The old way...
323            ### previously I was using kd.query_ball_point on, but the performance was terrible.
324            I = self.kd2.query_ball_point(grid_loc, radius)
325            for i in I:
326                t = tuple(self.kd.data[i])
327                possibles.update(self.hash[t])
328        return list(possibles)
329
330
331def random_segments(n):
332    segs = []
333    for i in range(n):
334        a, b, c, d = [random.random() for x in [1, 2, 3, 4]]
335        seg = LineSegment(Point((a, b)), Point((c, d)))
336        segs.append(seg)
337    return segs
338
339
340def random_points(n):
341    return [Point((random.random(), random.random())) for x in range(n)]
342
343
344def combo_check(bins, segments, qpoints):
345    G = SegmentLocator(segments, bins)
346    G2 = BruteSegmentLocator(segs)
347    for pt in qpoints:
348        a = G.nearest(pt)
349        b = G2.nearest(pt)
350        if a != b:
351            print(a, b, a == b)
352            global DEBUG
353            DEBUG = True
354            a = G.nearest(pt)
355            print(a)
356            a = segments[a]
357            b = segments[b]
358            print("pt to a (grid)", get_segment_point_dist(a, pt))
359            print("pt to b (brut)", get_segment_point_dist(b, pt))
360            input()
361            pylab.clf()
362            DEBUG = False
363
364
365def brute_check(segments, qpoints):
366    t0 = time.time()
367    G2 = BruteSegmentLocator(segs)
368    t1 = time.time()
369    print("Created Brute in %0.4f seconds" % (t1 - t0))
370    t2 = time.time()
371    q = list(map(G2.nearest, qpoints))
372    t3 = time.time()
373    print("Brute Found %d matches in %0.4f seconds" % (len(qpoints), t3 - t2))
374    print("Total Brute Time:", t3 - t0)
375    print()
376    return q
377
378
379def grid_check(bins, segments, qpoints, visualize=False):
380    t0 = time.time()
381    G = SegmentLocator(segments, bins)
382    t1 = time.time()
383    G.grid.kd
384    t2 = time.time()
385    print("Created Grid in %0.4f seconds" % (t1 - t0))
386    print("Created KDTree in %0.4f seconds" % (t2 - t1))
387    if visualize:
388        i = pylab.matshow(
389            G.grid.mask, origin="lower", extent=G.grid.x_range + G.grid.y_range
390        )
391
392    t2 = time.time()
393    q = list(map(G.nearest, qpoints))
394    t3 = time.time()
395    print("Grid Found %d matches in %0.4f seconds" % (len(qpoints), t3 - t2))
396    print("Total Grid Time:", t3 - t0)
397    qps = len(qpoints) / (t3 - t2)
398    print("q/s:", qps)
399    # print
400    return qps
401
402
403def binSizeTest():
404    q = 100
405    minN = 1000
406    maxN = 10000
407    stepN = 1000
408    minB = 250
409    maxB = 2000
410    stepB = 250
411    sizes = list(range(minN, maxN, stepN))
412    binSizes = list(range(minB, maxB, stepB))
413    results = numpy.zeros((len(sizes), len(binSizes)))
414    for row, n in enumerate(sizes):
415        segs = random_segments(n)
416        qpts = random_points(q)
417        for col, bins in enumerate(binSizes):
418            print("N, Bins:", n, bins)
419            qps = test_grid(bins, segs, qpts)
420            results[row, col] = qps
421    return results
422
423
424if __name__ == "__main__":
425    import pylab
426
427    pylab.ion()
428
429    n = 100
430    q = 1000
431
432    t0 = time.time()
433    segs = random_segments(n)
434    t1 = time.time()
435    qpts = random_points(q)
436    t2 = time.time()
437    print("segments:", t1 - t0)
438    print("points:", t2 - t1)
439    # test_brute(segs,qpts)
440    # test_grid(50, segs, qpts)
441
442    SG = SegmentLocator(segs)
443    grid = SG.grid
444