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