aisutils.grid

105 ''' 106 Provide a grid data structure. Will resize up to include the full 107 range on the max size if not an interger size of stepSize. Units are 108 whatever you desire. It is up to you to project them. 109 110 0,0 is at the lower left and (xNumCells-1,yNumCells-1) is the upper right cell 111 ''' 112 epsilon = .000001

113 - def __init__(self,minx,miny,maxx,maxy,stepSize,gridType='occurrence',verbose=False):

114 ''' Prepare a grid. 115 Readjust the grid such that the stepSize divides evenly into the ranges. 116 Compute and cache the number of cells. 117 ''' 118 self.minx=minx 119 self.miny=miny 120 self.maxx=maxx 121 self.maxy=maxy 122 self.verbose = verbose 123 stepSize=float(stepSize) 124 125 assert gridType in gridTypes 126 self.gridType = gridType 127 128 if verbose: 129 print ' == IN Params ==' 130 print ' x:',minx,maxx 131 print ' y:',miny,maxy 132 print ' step:',stepSize 133 134 assert(minx<maxx) 135 assert(miny<maxy) 136 assert(stepSize>0) 137 self.stepSize = stepSize 138 if (maxx-minx)/stepSize > floor((maxx-minx)/stepSize)+self.epsilon: 139 self.maxx = minx + stepSize*floor((maxx-minx)/stepSize) + stepSize 140 #self.maxx=maxx 141 142 if (maxy-miny)/stepSize > floor((maxy-miny)/stepSize)+self.epsilon: 143 self.maxy = miny + stepSize*floor((maxy-miny)/stepSize) + stepSize 144 #self.maxy=maxy 145 146 if verbose: 147 print ' Range changes:' 148 print ' maxx:',maxx,'->',self.maxx 149 print ' maxy:',maxy,'->',self.maxy 150 151 self.xNumCells = int(ceil((self.maxx-self.minx)/self.stepSize)) 152 self.yNumCells = int(ceil((self.maxy-self.miny)/self.stepSize)) 153 154 if verbose: 155 print ' step size num cells deltaMeters' 156 print ' x: ',self.stepSize,' ->',self.xNumCells, ' ', self.maxx - self.minx 157 print ' y: ',self.stepSize,' ->',self.yNumCells, ' ', self.maxy - self.miny 158 159 # FIX: why should I have to do add +1? Rounding/edge error? 160 # Will this cause errors down the road in other functions? 161 162 if gridType=='occurrence': 163 self.grid=numpy.zeros((self.xNumCells+1,self.yNumCells+1),dtype=int) 164 else: 165 self.grid=numpy.zeros((self.xNumCells+1,self.yNumCells+1),dtype=float)

166

167 - def describe(self):

168 print ' === GRID === ' 169 print ' Bounds ... (%.2f,%.2f) -> (%.2f,%.2f)'%(self.minx,self.miny,self.maxx,self.maxy) 170 print ' Width ...',self.xNumCells 171 print ' Height ...',self.yNumCells 172 print ' CellSize ...',self.stepSize

173

174 - def getCell(self,x,y,verbose=False):

175 '@return: the i,j of the cell containing this coordinate' 176 i = int(floor( (x-self.minx)/self.stepSize )) 177 j = int(floor( (y-self.miny)/self.stepSize )) 178 if verbose: print x,y,'->',i,j, (x-self.miny), self.stepSize 179 return i,j

180

181 - def getCellCenter(self,i,j):

182 '''@param i: the horizonal cell count from the left starting at 0 183 @return the x,y of the center of specified cell''' 184 x = self.minx + self.stepSize * i + 0.5*self.stepSize 185 y = self.miny + self.stepSize * j + 0.5*self.stepSize 186 return x,y

187

188 - def getLineCells2pt(self, p1, p2, verbose=False):

189 return self.getLineCells(p1[0],p1[1],p2[0],p2[1],verbose)

190 191

192 - def getLineCellsWithCrossingsPts(self, p1, p2, verbose=False):

193 return self.getLineCellsWithCrossings(p1[0],p1[1],p2[0],p2[1])

194

195 - def getLineCellsWithCrossings(self, x0,y0, x1,y1,verbose=False):

196 ''' 197 Scan convert a single line segment with two vertices. 198 199 Used for when we need to calculated paraeters off the line to 200 add a value to the cell. e.g. distance or speed wieghted by 201 distance. 202 203 A line inside a cell returns the cell, 0., 1., and the distance 204 205 @return: a list of (cell(an i,j), startFrac, endFrac, distance) for a line 206 @todo: make this actually be fast 207 @todo: switch to delta Y for N-S lines 208 ''' 209 if verbose: print '\n\n===============\n getLineCells',x0,y0, x1,y1 210 #from IPython.Shell import IPShellEmbed 211 #ipshell = IPShellEmbed(argv=[]) 212 #ipshell() 213 214 # only scan left to right direction 215 flippedX = False 216 if x0>x1: 217 if verbose: print 'flipping X' 218 x0,y0,x1,y1=x1,y1,x0,y0 219 flippedX = True 220 221 startCell = self.getCell(x0,y0,verbose) 222 endCell = self.getCell(x1,y1,verbose) 223 if verbose: print 'cellrange',startCell,endCell,x0,y0,'->',x1,y1 224 225 # just one cell 226 if startCell==endCell: 227 dist = distance(x0,y0, x1,y1) 228 return ( (startCell,0.,1.,dist),) 229 230 dx = x1-x0 231 dy = y1-y0 232 m = slope = (dy/float(dx)) 233 b = y0 - slope * x0 234 # Else, we have have to cross diagonally 235 236 stepSize = self.stepSize 237 minx=self.minx; maxx=self.maxx 238 miny=self.miny; maxy=self.maxy 239 240 flippedY = False 241 if y0>y1: 242 if verbose: print 'flipping Y' 243 flippedY = True 244 y0,y1=y1,y0 245 y_first_ycrossing = miny + ceil ((y0 - miny)/stepSize) * stepSize 246 y_last_ycrossing = miny + floor((y1 - miny)/stepSize) * stepSize 247 steps = int( ceil((y_last_ycrossing - y_first_ycrossing)/stepSize + 1)) 248 249 if verbose: 250 print 'steps y',y_first_ycrossing,y_last_ycrossing ,'->',steps 251 252 x_ycrossings = [ ] 253 for step in range(steps): 254 y = step*stepSize + y_first_ycrossing 255 x = (y-b)/m 256 x_ycrossings.append(x) 257 del y_first_ycrossing 258 del y_last_ycrossing 259 260 if verbose: print 'x_ycrossings',x_ycrossings 261 262 x_first_xcrossing = minx + ceil((x0 - minx)/stepSize) * stepSize 263 if verbose: 264 print ' x_xcross params:', x0,minx, stepSize, ceil((x0 - minx)/stepSize), 265 print '--->', x_first_xcrossing 266 x_last_xcrossing = minx + floor((x1 - minx)/stepSize) * stepSize 267 268 steps = int( ceil((x_last_xcrossing - x_first_xcrossing)/float(stepSize) +1) ) 269 if verbose: 270 print 'steps2',x_first_xcrossing,x_last_xcrossing ,'->',steps 271 print 'steps2',x_last_xcrossing - x_first_xcrossing, (x_last_xcrossing - x_first_xcrossing)/stepSize 272 x_xcrossings = [ ] 273 for step in range(steps): 274 x = step*stepSize + x_first_xcrossing 275 if verbose: print ' X:', step,stepSize,x_first_xcrossing 276 x_xcrossings.append(x) 277 if verbose: print 'x_xcrossings',x_xcrossings 278 279 crossings = x_ycrossings+x_xcrossings 280 if flippedX: 281 if verbose: print 'reverse sort crossings' 282 crossings.sort(reverse=True) 283 else: 284 crossings.sort() 285 if verbose: print '\ncrossings',crossings 286 287 # remove duplicates 288 cnew=[crossings[0],] 289 for i in range(1,len(crossings)): 290 #if crossings[i]==crossings[i-1]: continue # Faster... first? 291 if almostEqual(crossings[i],crossings[i-1]): continue 292 cnew.append(crossings[i]) 293 #print crossings,'->',cnew 294 if verbose: 295 print 'orig crossings:',crossings 296 print ' new crossings:',cnew 297 crossings=cnew 298 299 # Be careful which side to add that start cell 300 cells=[] 301 if not flippedX: 302 cells = [startCell,] 303 304 fractions = [0.,] 305 306 totXRange = x1 - x0 307 print 'xrange',totXRange , x1,x0 308 for x in crossings: 309 fractions.append( (x-x0) / totXRange ) 310 x = x+0.0001 311 y = m * x + b 312 cells.append(self.getCell(x,y)) 313 if flippedX: 314 cells.append(startCell) 315 316 fractions.append(1.) 317 print 'fractions',fractions 318 distances = [] 319 for i in range(0,len(fractions)-1): 320 x = x0 + totXRange * fractions[i] 321 y = m * x + b 322 p1 = (x,y) 323 x = x0 + totXRange * fractions[i+1] 324 y = m * x + b 325 p2 = (x,y) 326 327 distances.append(distancePt(p1,p2)) 328 print 'dists',distances 329 330 331 assert(len(cells)>1) 332 # Make sure the line does not extend beyond the cell range 333 minCellX = min(startCell[0],endCell[0]) 334 maxCellX = max(startCell[0],endCell[0]) 335 minCellY = min(startCell[1],endCell[1]) 336 maxCellY = max(startCell[1],endCell[1]) 337 338 print 'minmaxX',minCellX,maxCellX 339 print 'minmaxY',minCellY,maxCellY 340 341 342 # if rounding errors happen to push us out, then toss a cell 343 344 if not inBoundingBox(cells[0][0],cells[0][1],minCellX,minCellY,maxCellX,maxCellY): 345 if verbose: print 'chomp front' 346 assert False # need to also remove from the fractions 347 cells=cells[1:] 348 if not inBoundingBox(cells[-1][0],cells[-1][1],minCellX,minCellY,maxCellX,maxCellY): 349 if verbose: print 'chomp tail' 350 assert False # need to also remove from the fractions 351 cells=cells[:-1] 352 353 #c = cells[0] 354 #if not inRange(c[0],minCellX,maxCellX) or not inRange(c[1],minCellY,maxCellY): 355 # if verbose: print 'chomp front' 356 # cells=cells[1:] 357 #c = cells[-1] 358 #if not inRange(c[0],minCellX,maxCellX) or not inRange(c[1],minCellY,maxCellY): 359 # 360 # cells=cells[:-1] 361 362 assert(len(cells)>1) 363 if verbose: print cells 364 #return cells 365 r = [] 366 print 'lens', 367 assert len(fractions)==len(cells)+1 368 assert len(cells)==len(distances) 369 for i in range(len(cells)): 370 r.append((cells[i],fractions[i],fractions[i+1],distances[i])) 371 return r

372 373 ######################################################################

374 - def getLineCells(self, x0,y0, x1,y1,verbose=False):

375 ''' 376 Scan convert a single line segment with two vertices. 377 @return: a list of cells for a line 378 @todo: make this actually be fast 379 ''' 380 if verbose: print '\n\n===============\n getLineCells',x0,y0, x1,y1 381 #from IPython.Shell import IPShellEmbed 382 #ipshell = IPShellEmbed(argv=[]) 383 #ipshell() 384 385 # only scan left to right direction 386 flippedX = False 387 if x0>x1: 388 if verbose: print 'flipping X' 389 x0,y0,x1,y1=x1,y1,x0,y0 390 flippedX = True 391 392 startCell = self.getCell(x0,y0,verbose) 393 endCell = self.getCell(x1,y1,verbose) 394 if verbose: print 'cellrange',startCell,endCell,x0,y0,'->',x1,y1 395 396 # Vertical Line or just one cell? 397 if startCell[0]==endCell[0]: 398 yFlipped=False 399 if startCell[1]>endCell[1]: 400 startCell,endCell=endCell,startCell 401 yFlipped=True 402 if verbose: 403 print 'CASE vert or single cell', startCell[1],endCell[1]+1 404 i = startCell[0] 405 cells = [ (i,j) for j in range(startCell[1],endCell[1]+1) ] 406 if yFlipped: 407 cells.reverse() 408 return cells 409 410 # Horizontal Line? Do this fast 411 if startCell[1]==endCell[1]: 412 if verbose: print 'CASE horiz' 413 j = startCell[1] 414 #print range(startCell[0],endCell[0]+1) 415 cells = [ (i,j) for i in range(startCell[0],endCell[0]+1) ] 416 if flippedX: 417 # Line goes from right to left 418 if verbose: print 'reversing based on x (horiz line)' 419 cells.reverse() 420 #print 'going to return',cells 421 return cells 422 423 if verbose: print 'CASE not easy... must calculate' 424 425 dx = x1-x0 426 dy = y1-y0 427 m = slope = (dy/float(dx)) 428 b = y0 - slope * x0 429 # Else, we have have to cross diagonally 430 431 stepSize = self.stepSize 432 minx=self.minx; maxx=self.maxx 433 miny=self.miny; maxy=self.maxy 434 435 flippedY = False 436 if y0>y1: 437 if verbose: print 'flipping Y' 438 flippedY = True 439 y0,y1=y1,y0 440 y_first_ycrossing = miny + ceil ((y0 - miny)/stepSize) * stepSize 441 y_last_ycrossing = miny + floor((y1 - miny)/stepSize) * stepSize 442 steps = int( ceil((y_last_ycrossing - y_first_ycrossing)/stepSize + 1)) 443 444 if verbose: 445 print 'steps y',y_first_ycrossing,y_last_ycrossing ,'->',steps 446 447 x_ycrossings = [ ] 448 for step in range(steps): 449 y = step*stepSize + y_first_ycrossing 450 x = (y-b)/m 451 x_ycrossings.append(x) 452 del y_first_ycrossing 453 del y_last_ycrossing 454 455 if verbose: print 'x_ycrossings',x_ycrossings 456 457 x_first_xcrossing = minx + ceil((x0 - minx)/stepSize) * stepSize 458 if verbose: 459 print ' x_xcross params:', x0,minx, stepSize, ceil((x0 - minx)/stepSize), 460 print '--->', x_first_xcrossing 461 x_last_xcrossing = minx + floor((x1 - minx)/stepSize) * stepSize 462 463 steps = int( ceil((x_last_xcrossing - x_first_xcrossing)/float(stepSize) +1) ) 464 if verbose: 465 print 'steps2',x_first_xcrossing,x_last_xcrossing ,'->',steps 466 print 'steps2',x_last_xcrossing - x_first_xcrossing, (x_last_xcrossing - x_first_xcrossing)/stepSize 467 x_xcrossings = [ ] 468 for step in range(steps): 469 x = step*stepSize + x_first_xcrossing 470 if verbose: print ' X:', step,stepSize,x_first_xcrossing 471 x_xcrossings.append(x) 472 if verbose: print 'x_xcrossings',x_xcrossings 473 474 crossings = x_ycrossings+x_xcrossings 475 if flippedX: 476 if verbose: print 'reverse sort crossings' 477 crossings.sort(reverse=True) 478 else: 479 crossings.sort() 480 if verbose: print '\ncrossings',crossings 481 482 # remove duplicates 483 cnew=[crossings[0],] 484 for i in range(1,len(crossings)): 485 #if crossings[i]==crossings[i-1]: continue # Faster... first? 486 if almostEqual(crossings[i],crossings[i-1]): continue 487 cnew.append(crossings[i]) 488 #print crossings,'->',cnew 489 if verbose: 490 print 'orig crossings:',crossings 491 print ' new crossings:',cnew 492 crossings=cnew 493 494 495 # Be careful which side to add that start cell 496 cells=[] 497 if not flippedX: 498 cells = [startCell,] 499 for x in crossings: 500 x = x+0.0001 501 y = m * x + b 502 cells.append(self.getCell(x,y)) 503 if flippedX: 504 cells.append(startCell) 505 506 assert(len(cells)>1) 507 # Make sure the line does not extend beyond the cell range 508 minCellX = min(startCell[0],endCell[0]) 509 maxCellX = max(startCell[0],endCell[0]) 510 minCellY = min(startCell[1],endCell[1]) 511 maxCellY = max(startCell[1],endCell[1]) 512 513 c = cells[0] 514 if not inRange(c[0],minCellX,maxCellX) or not inRange(c[1],minCellY,maxCellY): 515 if verbose: print 'chomp front' 516 cells=cells[1:] 517 c = cells[-1] 518 if not inRange(c[0],minCellX,maxCellX) or not inRange(c[1],minCellY,maxCellY): 519 if verbose: print 'chomp tail' 520 cells=cells[:-1] 521 522 assert(len(cells)>1) 523 if verbose: print cells 524 return cells

525

526 - def getMultiSegLineCells(self,multiSegLine,verbose=False):

527 ''' 528 Return the cells for a multi vertex line. Handles the doubling that will happen at each endpoint 529 @param multiSegLine: ((x1,y1),(x2,y2),(x3,y3),(x4,y4)...) 530 ''' 531 assert len(multiSegLine)>1 532 533 cells = self.getLineCells2pt(multiSegLine[0],multiSegLine[1]) 534 if len(multiSegLine)==2: 535 return cells 536 for i in range(1,len(multiSegLine)-1): 537 # Make sure to skip the first cell on each additional segment 538 newCells = self.getLineCells2pt(multiSegLine[i],multiSegLine[i+1]) 539 if len(newCells)<1: 540 print 'ACK... must have at least one cell!!! ERROR fail death' 541 print ' ',i,multiSegLine[i],multiSegLine[i+1] 542 print ' ',newCells 543 assert False 544 if newCells[0] == cells[-1]: 545 cells+=newCells[1:] 546 else: 547 cells+=newCells[:-1] 548 return cells

549 550

551 - def getMultiSegLineCellsWithCrossings(self,multiSegLine,verbose=False):

552 ''' 553 Return the cells for a multi vertex line. Returns distances 554 within cells, so that doubling is not an issue 555 556 @param multiSegLine: ((x1,y1),(x2,y2),(x3,y3),(x4,y4)...) 557 @return: list of (cell, startFrac, endFrac, distance) 558 ''' 559 assert len(multiSegLine)>1 560 561 #cells = getLineCellsWithCrossingsPts(multiSegLine[0],multiSegLine[1]) 562 cells=[] 563 for i in range(len(multiSegLine)-1): 564 newCells = self.getLineCellsWithCrossingsPts(multiSegLine[i],multiSegLine[i+1]) 565 cells+=newCells 566 return cells

567 568

569 - def writeLayoutGnuplot(self,filename):

570 'Write out the grid lines as gnuplot dat file' 571 o = file(filename,'w') 572 # Horizonal lines 573 minxStr = str(self.minx)+' ' 574 maxxStr = str(self.maxx)+' ' 575 for yCell in range(self.yNumCells+1): 576 yStr = str(self.miny + yCell*self.stepSize)+' 0\n' 577 o.write(minxStr+yStr) 578 o.write(maxxStr+yStr+'\n') 579 # Vertical Lines 580 minyStr = ' '+str(self.miny)+' 0\n' 581 maxyStr = ' '+str(self.maxy)+' 0\n' 582 for xCell in range(self.xNumCells+1): 583 xStr = str(self.minx + xCell*self.stepSize) 584 o.write(xStr+minyStr) 585 o.write(xStr+maxyStr+'\n')

586 587

588 - def addMultiSegLine(self,multiSegLine,verbose=False):

589 grid=self.grid 590 #print '\naddMultiSegLine line ',multiSegLine 591 if verbose: sys.stderr.write('addMultiSegLine cell inserts: (using type '+self.gridType+')\n') 592 593 if 'occurrence' == self.gridType: 594 cells = self.getMultiSegLineCells(multiSegLine) 595 for cell in cells: 596 try: 597 grid[cell[0],cell[1]]+=1 598 if verbose: 599 print ' ',cell[0],cell[1],grid[cell[0],cell[1]] 600 except: 601 print multiSegLine 602 print 'CRAP... grid failure' 603 print ' ',cell, self.xNumCells, self.yNumCells 604 assert False 605 elif 'distance' == self.gridType: 606 result = self.getMultiSegLineCellsWithCrossings(multiSegLine) 607 if verbose: sys.stderr.write(str(result)+'\n') 608 for seg in result: 609 # (0, 0), 0.0, 0.24874371859296465, 0.99124918032832443) 610 cell = seg[0] 611 dist = seg[3] 612 print 'adding:', cell, dist 613 grid[cell[0],cell[1]] += dist 614 print grid 615 elif 'distanceWeightedSpeed' == self.gridType: 616 assert False 617 else: 618 assert False

619

620 - def writeCellsGnuplot(self,filename,useSquares=False):

621 ''' 622 @param useSquares: if true then write out the height of each cell as a square. False then it writes a point 623 ''' 624 assert useSquares==False # FIX: implement this feature 625 grid = self.grid 626 o = file(filename,'w') 627 for i in range(grid.shape[0]): 628 for j in range(grid.shape[1]): 629 x,y = self.getCellCenter(i,j) 630 o.write('%f %f %f\n' % (x,y,grid[i,j]))

631 632

633 - def writeArcAsciiGrid(self,filename):

634 g = self.grid 635 o = file(filename,'w') 636 o.write('ncols '+str(self.xNumCells)+'\n') 637 o.write('nrows '+str(self.yNumCells)+'\n') 638 o.write('xllcorner '+str(self.minx)+'\n') 639 o.write('yllcorner '+str(self.miny)+'\n') 640 o.write('cellsize '+str(self.stepSize)+'\n') 641 for j in range(self.yNumCells-1,-1,-1): 642 #print j 643 zPoints=[] 644 for i in range(self.xNumCells): 645 zPoints.append('%3d' % (g[i,j])) 646 o.write(' '.join(zPoints)) 647 o.write('\n')

Source Code for Module aisutils.grid