Changeset 1289
- Timestamp:
- May 6, 2005, 10:16:02 AM (20 years ago)
- File:
-
- 1 edited
-
inundation/ga/storm_surge/pyvolution/util.py (modified) (69 diffs)
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/util.py
r1280 r1289 2 2 3 3 It is also a clearing house for functions that may later earn a module 4 of their own. 4 of their own. 5 5 """ 6 6 … … 17 17 v2 = v[1]/l 18 18 19 try: 19 try: 20 20 theta = acos(v1) 21 21 except ValueError, e: … … 31 31 theta = 3.1415926535897931 32 32 print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\ 33 %(v1, theta) 34 33 %(v1, theta) 34 35 35 if v2 < 0: 36 36 #Quadrant 3 or 4 … … 43 43 """Compute difference between angle of vector x0, y0 and x1, y1. 44 44 This is used for determining the ordering of vertices, 45 e.g. for checking if they are counter clockwise. 46 45 e.g. for checking if they are counter clockwise. 46 47 47 Always return a positive value 48 48 """ 49 49 50 50 from math import pi 51 51 52 52 a0 = angle(v0) 53 53 a1 = angle(v1) … … 56 56 if a0 < a1: 57 57 a0 += 2*pi 58 58 59 59 return a0-a1 60 60 … … 65 65 66 66 67 def point_on_line(x, y, x0, y0, x1, y1): 68 """Determine whether a point is on a line segment 69 67 def point_on_line(x, y, x0, y0, x1, y1): 68 """Determine whether a point is on a line segment 69 70 70 Input: x, y, x0, x0, x1, y1: where 71 71 point is given by x, y 72 72 line is given by (x0, y0) and (x1, y1) 73 74 """ 75 73 74 """ 75 76 76 from Numeric import array, dot, allclose 77 77 from math import sqrt 78 78 79 a = array([x - x0, y - y0]) 79 a = array([x - x0, y - y0]) 80 80 a_normal = array([a[1], -a[0]]) 81 81 82 82 b = array([x1 - x0, y1 - y0]) 83 83 84 84 if dot(a_normal, b) == 0: 85 85 #Point is somewhere on the infinite extension of the line 86 86 87 len_a = sqrt(sum(a**2)) 88 len_b = sqrt(sum(b**2)) 87 len_a = sqrt(sum(a**2)) 88 len_b = sqrt(sum(b**2)) 89 89 if dot(a, b) >= 0 and len_a <= len_b: 90 90 return True 91 else: 91 else: 92 92 return False 93 93 else: 94 return False 95 94 return False 95 96 96 def ensure_numeric(A, typecode = None): 97 97 """Ensure that sequence is a Numeric array. … … 117 117 if type(A) == ArrayType: 118 118 if A.typecode == typecode: 119 return array(A) 119 return array(A) 120 120 else: 121 121 return A.astype(typecode) 122 122 else: 123 123 return array(A).astype(typecode) 124 125 124 125 126 126 127 127 def file_function(filename, domain=None, quantities = None, interpolation_points = None, verbose = False): … … 133 133 #In fact, this is where origin should be converted to that of domain 134 134 #Also, check that file covers domain fully. 135 135 136 136 assert type(filename) == type(''),\ 137 137 'First argument to File_function must be a string' … … 148 148 149 149 if domain is not None: 150 quantities = domain.conserved_quantities 150 quantities = domain.conserved_quantities 151 151 else: 152 152 quantities = None 153 154 #Choose format 153 154 #Choose format 155 155 #FIXME: Maybe these can be merged later on 156 156 if line[:3] == 'CDF': 157 157 return File_function_NetCDF(filename, domain, quantities, interpolation_points, verbose = verbose) 158 158 else: 159 return File_function_ASCII(filename, domain, quantities, interpolation_points) 160 161 159 return File_function_ASCII(filename, domain, quantities, interpolation_points) 160 161 162 162 class File_function_NetCDF: 163 """Read time history of spatial data from NetCDF sww file and 164 return a callable object f(t,x,y) 165 which will return interpolated values based on the input file. 166 163 """Read time history of spatial data from NetCDF sww file and 164 return a callable object f(t,x,y) 165 which will return interpolated values based on the input file. 166 167 167 x, y may be either scalars or vectors 168 168 169 169 #FIXME: More about format, interpolation and order of quantities 170 170 171 171 The quantities returned by the callable objects are specified by 172 the list quantities which must contain the names of the 172 the list quantities which must contain the names of the 173 173 quantities to be returned and also reflect the order, e.g. for 174 the shallow water wave equation, on would have 174 the shallow water wave equation, on would have 175 175 quantities = ['stage', 'xmomentum', 'ymomentum'] 176 177 interpolation_points decides at which points interpolated 176 177 interpolation_points decides at which points interpolated 178 178 quantities are to be computed whenever object is called. 179 180 181 179 180 181 182 182 If None, return average value 183 183 """ 184 184 185 185 def __init__(self, filename, domain=None, quantities=None, 186 186 interpolation_points=None, verbose = False): … … 189 189 If domain is specified, model time (domain.starttime) 190 190 will be checked and possibly modified 191 191 192 192 All times are assumed to be in UTC 193 193 """ … … 196 196 #(both in UTM coordinates) 197 197 #If not - modify those from file to match domain 198 198 199 199 import time, calendar 200 200 from config import time_format 201 from Scientific.IO.NetCDF import NetCDFFile 201 from Scientific.IO.NetCDF import NetCDFFile 202 202 203 203 #Open NetCDF file … … 213 213 if not fid.variables.has_key(quantity): 214 214 missing.append(quantity) 215 215 216 216 if len(missing) > 0: 217 217 msg = 'Quantities %s could not be found in file %s'\ 218 218 %(str(missing), filename) 219 219 raise msg 220 220 221 221 222 222 #Get first timestep 223 223 try: 224 224 self.starttime = fid.starttime[0] 225 except ValueError: 225 except ValueError: 226 226 msg = 'Could not read starttime from file %s' %filename 227 227 raise msg … … 229 229 #Get origin 230 230 self.xllcorner = fid.xllcorner[0] 231 self.yllcorner = fid.yllcorner[0] 232 231 self.yllcorner = fid.yllcorner[0] 232 233 233 self.number_of_values = len(quantities) 234 234 self.fid = fid 235 235 self.filename = filename 236 237 238 236 self.quantities = quantities 237 self.domain = domain 238 self.interpolation_points = interpolation_points 239 239 240 240 if domain is not None: … … 249 249 if verbose: print 'Domain starttime is now set to %f' %domain.starttime 250 250 251 #Read all data in and produce values for desired data points at each timestep 252 self.spatial_interpolation(interpolation_points, verbose = verbose) 251 #Read all data in and produce values for desired data points at each timestep 252 self.spatial_interpolation(interpolation_points, verbose = verbose) 253 253 fid.close() 254 254 255 255 def spatial_interpolation(self, interpolation_points, verbose = False): 256 """For each timestep in fid: read surface, interpolate to desired points 257 256 """For each timestep in fid: read surface, interpolate to desired points 257 and store values for use when object is called. 258 258 """ 259 259 260 260 from Numeric import array, zeros, Float, alltrue, concatenate, reshape 261 261 262 262 from config import time_format 263 263 from least_squares import Interpolation 264 264 import time, calendar 265 265 266 266 fid = self.fid 267 267 … … 273 273 triangles = fid.variables['volumes'][:] 274 274 time = fid.variables['time'][:] 275 276 275 276 #Check 277 277 msg = 'File %s must list time as a monotonuosly ' %self.filename 278 278 msg += 'increasing sequence' 279 assert alltrue(time[1:] - time[:-1] > 0 ), msg 280 281 282 if interpolation_points is not None: 283 279 assert alltrue(time[1:] - time[:-1] > 0 ), msg 280 281 282 if interpolation_points is not None: 283 284 284 try: 285 P = ensure_numeric(interpolation_points) 285 P = ensure_numeric(interpolation_points) 286 286 except: 287 287 msg = 'Interpolation points must be an N x 2 Numeric array or a list of points\n' 288 288 msg += 'I got: %s.' %( str(interpolation_points)[:60] + '...') 289 289 raise msg 290 291 290 291 292 292 self.T = time[:] #Time assumed to be relative to starttime 293 self.index = 0 #Initial time index 294 295 self.values = {} 293 self.index = 0 #Initial time index 294 295 self.values = {} 296 296 for name in self.quantities: 297 self.values[name] = zeros( (len(self.T), 298 len(interpolation_points)), 299 Float) 297 self.values[name] = zeros( (len(self.T), 298 len(interpolation_points)), 299 Float) 300 300 301 301 #Build interpolator 302 302 if verbose: print 'Build interpolation matrix' 303 303 x = reshape(x, (len(x),1)) 304 y = reshape(y, (len(y),1)) 304 y = reshape(y, (len(y),1)) 305 305 vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array 306 306 307 interpol = Interpolation(vertex_coordinates, 307 interpol = Interpolation(vertex_coordinates, 308 308 triangles, 309 309 point_coordinates = P, 310 alpha = 0, 310 alpha = 0, 311 311 verbose = verbose) 312 312 313 313 314 314 if verbose: print 'Interpolate' 315 315 for i, t in enumerate(self.T): 316 316 #Interpolate quantities at this timestep 317 317 #print ' time step %d of %d' %(i, len(self.T)) 318 318 for name in self.quantities: 319 320 interpol.interpolate(fid.variables[name][i,:])319 self.values[name][i, :] =\ 320 interpol.interpolate(fid.variables[name][i,:]) 321 321 322 322 #Report … … 327 327 print ' Reference:' 328 328 print ' Lower left corner: [%f, %f]'\ 329 %(self.xllcorner, self.yllcorner) 329 %(self.xllcorner, self.yllcorner) 330 330 print ' Start time (file): %f' %self.starttime 331 331 #print ' Start time (domain): %f' %domain.starttime … … 342 342 print ' %s in [%f, %f]' %(name, min(q), max(q)) 343 343 344 344 345 345 print ' Interpolation points (xi, eta):'\ 346 346 'number of points == %d ' %P.shape[0] … … 349 349 print ' Interpolated quantities (over all timesteps):' 350 350 for name in self.quantities: 351 q = self.values[name][:].flat 351 q = self.values[name][:].flat 352 352 print ' %s at interpolation points in [%f, %f]'\ 353 %(name, min(q), max(q)) 353 %(name, min(q), max(q)) 354 355 356 354 357 355 358 print '------------------------------------------------' … … 364 367 def __call__(self, t, x=None, y=None, point_id = None): 365 368 """Evaluate f(t, point_id) 366 367 368 t: time - Model time (tau = domain.starttime-self.starttime+t) 369 must lie within existing timesteps 370 point_id: index of one of the preprocessed points. 371 372 373 FIXME: point_id could also be a slice 374 FIXME: One could allow arbitrary x, y coordinates 375 376 377 FIXME: What if x and y are vectors? 369 370 Inputs: 371 t: time - Model time (tau = domain.starttime-self.starttime+t) 372 must lie within existing timesteps 373 point_id: index of one of the preprocessed points. 374 If point_id is None all preprocessed points are computed 375 376 FIXME: point_id could also be a slice 377 FIXME: One could allow arbitrary x, y coordinates 378 (requires computation of a new interpolator) 379 Maybe not,.,. 380 FIXME: What if x and y are vectors? 378 381 """ 379 382 … … 388 391 389 392 #Find time tau such that 390 391 392 393 393 # 394 # domain.starttime + t == self.starttime + tau 395 396 if self.domain is not None: 394 397 tau = self.domain.starttime-self.starttime+t 395 398 else: 396 399 #print 'DOMAIN IS NONE!!!!!!!!!' 397 400 tau = t … … 400 403 #print 'S start', self.starttime 401 404 #print 't', t 402 #print 'tau', tau 403 405 #print 'tau', tau 406 404 407 msg = 'Time interval derived from file %s [%s:%s]'\ 405 408 %(self.filename, self.T[0], self.T[1]) … … 408 411 409 412 if tau < self.T[0]: raise msg 410 if tau > self.T[-1]: raise msg 413 if tau > self.T[-1]: raise msg 411 414 412 415 … … 421 424 #Protect against case where tau == T[-1] (last time) 422 425 # - also works in general when tau == T[i] 423 ratio = 0 426 ratio = 0 424 427 else: 425 #t is now between index and index+1 428 #t is now between index and index+1 426 429 ratio = (tau - self.T[self.index])/\ 427 430 (self.T[self.index+1] - self.T[self.index]) … … 429 432 430 433 431 434 432 435 #Compute interpolated values 433 436 q = zeros( len(self.quantities), Float) 434 437 435 438 for i, name in enumerate(self.quantities): 436 Q = self.values[name] 439 Q = self.values[name] 437 440 438 441 if ratio > 0: … … 447 450 #Return vector of interpolated values 448 451 return q 449 450 452 453 451 454 class File_function_ASCII: 452 455 """Read time series from file and return a callable object: 453 454 455 456 457 458 time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...459 460 461 462 time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 #FIXME: This should work with netcdf (e.g. sww) and thus render the482 483 484 #FIXME: Specified quantites not used here -485 486 """ 487 488 456 f(t,x,y) 457 which will return interpolated values based on the input file. 458 459 The file format is assumed to be either two fields separated by a comma: 460 461 time [DD/MM/YY hh:mm:ss], value0 value1 value2 ... 462 463 or four comma separated fields 464 465 time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ... 466 467 In either case, the callable object will return a tuple of 468 interpolated values, one each value stated in the file. 469 470 471 E.g 472 473 31/08/04 00:00:00, 1.328223 0 0 474 31/08/04 00:15:00, 1.292912 0 0 475 476 will provide a time dependent function f(t,x=None,y=None), 477 ignoring x and y, which returns three values per call. 478 479 480 NOTE: At this stage the function is assumed to depend on 481 time only, i.e no spatial dependency!!!!! 482 When that is needed we can use the least_squares interpolation. 483 484 #FIXME: This should work with netcdf (e.g. sww) and thus render the 485 #spatio-temporal boundary condition in shallow water fully general 486 487 #FIXME: Specified quantites not used here - 488 #always return whatever is in the file 489 """ 490 491 489 492 def __init__(self, filename, domain=None, quantities = None, interpolation_points=None): 490 493 """Initialise callable object from a file. … … 504 507 line = fid.readline() 505 508 fid.close() 506 509 507 510 fields = line.split(',') 508 511 msg = 'File %s must have the format date, value0 value1 value2 ...' … … 513 516 try: 514 517 starttime = calendar.timegm(time.strptime(fields[0], time_format)) 515 except ValueError: 518 except ValueError: 516 519 msg = 'First field in file %s must be' %filename 517 520 msg += ' date-time with format %s.\n' %time_format 518 msg += 'I got %s instead.' %fields[0] 521 msg += 'I got %s instead.' %fields[0] 519 522 raise msg 520 523 … … 526 529 527 530 q = ensure_numeric(values) 528 531 529 532 msg = 'ERROR: File must contain at least one independent value' 530 533 assert len(q.shape) == 1, msg 531 534 532 535 self.number_of_values = len(q) 533 536 534 537 self.filename = filename 535 538 self.starttime = starttime 536 539 self.domain = domain 537 540 538 541 if domain is not None: … … 546 549 ##if verbose: print msg 547 550 domain.starttime = self.starttime #Modifying model time 548 551 549 552 #if domain.starttime is None: 550 553 # domain.starttime = self.starttime … … 560 563 # domain.starttime = self.starttime #Modifying model time 561 564 562 if mode == 2: 565 if mode == 2: 563 566 self.read_times() #Now read all times in. 564 567 else: 565 568 raise 'x,y dependency not yet implemented' 566 569 567 570 568 571 def read_times(self): 569 """Read time and values 572 """Read time and values 570 573 """ 571 574 from Numeric import zeros, Float, alltrue 572 575 from config import time_format 573 576 import time, calendar 574 575 fid = open(self.filename) 577 578 fid = open(self.filename) 576 579 lines = fid.readlines() 577 580 fid.close() 578 581 579 582 N = len(lines) 580 583 d = self.number_of_values … … 582 585 T = zeros(N, Float) #Time 583 586 Q = zeros((N, d), Float) #Values 584 587 585 588 for i, line in enumerate(lines): 586 589 fields = line.split(',') … … 600 603 self.index = 0 #Initial index 601 604 602 605 603 606 def __repr__(self): 604 607 return 'File function' … … 610 613 result is a vector of same length as x and y 611 614 FIXME: Naaaa 612 613 615 616 FIXME: Rethink semantics when x,y are vectors. 614 617 """ 615 618 616 619 from math import pi, cos, sin, sqrt 617 620 618 621 619 622 #Find time tau such that 620 621 622 623 623 # 624 # domain.starttime + t == self.starttime + tau 625 626 if self.domain is not None: 624 627 tau = self.domain.starttime-self.starttime+t 625 628 else: 626 629 tau = t 627 628 630 631 629 632 msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\ 630 633 %(self.filename, self.T[0], self.T[1], tau) 631 634 if tau < self.T[0]: raise msg 632 if tau > self.T[-1]: raise msg 635 if tau > self.T[-1]: raise msg 633 636 634 637 oldindex = self.index … … 644 647 #Compute interpolated values 645 648 q = self.Q[self.index,:] +\ 646 ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 649 ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 647 650 648 651 #Return vector of interpolated values … … 659 662 N = len(x) 660 663 assert len(y) == N, 'x and y must have same length' 661 664 662 665 res = [] 663 666 for col in q: 664 667 res.append(col*ones(N, Float)) 665 668 666 669 return res 667 668 669 #FIXME: TEMP 670 671 672 #FIXME: TEMP 670 673 class File_function_copy: 671 674 """Read time series from file and return a callable object: 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 #FIXME: This should work with netcdf (e.g. sww) and thus render the701 702 """ 703 704 675 f(t,x,y) 676 which will return interpolated values based on the input file. 677 678 The file format is assumed to be either two fields separated by a comma: 679 680 time [DD/MM/YY hh:mm:ss], value0 value1 value2 ... 681 682 or four comma separated fields 683 684 time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ... 685 686 In either case, the callable object will return a tuple of 687 interpolated values, one each value stated in the file. 688 689 690 E.g 691 692 31/08/04 00:00:00, 1.328223 0 0 693 31/08/04 00:15:00, 1.292912 0 0 694 695 will provide a time dependent function f(t,x=None,y=None), 696 ignoring x and y, which returns three values per call. 697 698 699 NOTE: At this stage the function is assumed to depend on 700 time only, i.e no spatial dependency!!!!! 701 When that is needed we can use the least_squares interpolation. 702 703 #FIXME: This should work with netcdf (e.g. sww) and thus render the 704 #spatio-temporal boundary condition in shallow water fully general 705 """ 706 707 705 708 def __init__(self, filename, domain=None): 706 709 """Initialise callable object from a file. … … 739 742 try: 740 743 starttime = calendar.timegm(time.strptime(fields[0], time_format)) 741 except ValueError: 744 except ValueError: 742 745 msg = 'First field in file %s must be' %filename 743 746 msg += ' date-time with format %s.\n' %time_format 744 msg += 'I got %s instead.' %fields[0] 747 msg += 'I got %s instead.' %fields[0] 745 748 raise msg 746 749 … … 752 755 753 756 q = ensure_numeric(values) 754 757 755 758 msg = 'ERROR: File must contain at least one independent value' 756 759 assert len(q.shape) == 1, msg 757 760 758 761 self.number_of_values = len(q) 759 762 760 763 self.filename = filename 761 764 self.starttime = starttime 762 765 self.domain = domain 763 766 764 767 if domain is not None: … … 776 779 domain.starttime = self.starttime #Modifying model time 777 780 778 if mode == 2: 781 if mode == 2: 779 782 self.read_times() #Now read all times in. 780 783 else: 781 784 raise 'x,y dependency not yet implemented' 782 785 783 786 784 787 def read_times(self): 785 """Read time and values 788 """Read time and values 786 789 """ 787 790 from Numeric import zeros, Float, alltrue 788 791 from config import time_format 789 792 import time, calendar 790 791 fid = open(self.filename) 793 794 fid = open(self.filename) 792 795 lines = fid.readlines() 793 796 fid.close() 794 797 795 798 N = len(lines) 796 799 d = self.number_of_values … … 798 801 T = zeros(N, Float) #Time 799 802 Q = zeros((N, d), Float) #Values 800 803 801 804 for i, line in enumerate(lines): 802 805 fields = line.split(',') … … 816 819 self.index = 0 #Initial index 817 820 818 821 819 822 def __repr__(self): 820 823 return 'File function' 821 824 822 825 823 826 824 827 def __call__(self, t, x=None, y=None): … … 831 834 832 835 from math import pi, cos, sin, sqrt 833 836 834 837 835 838 #Find time tau such that 836 837 838 839 839 # 840 # domain.starttime + t == self.starttime + tau 841 842 if self.domain is not None: 840 843 tau = self.domain.starttime-self.starttime+t 841 844 else: 842 845 tau = t 843 844 846 847 845 848 msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\ 846 849 %(self.filename, self.T[0], self.T[1], tau) 847 850 if tau < self.T[0]: raise msg 848 if tau > self.T[-1]: raise msg 851 if tau > self.T[-1]: raise msg 849 852 850 853 oldindex = self.index … … 860 863 #Compute interpolated values 861 864 q = self.Q[self.index,:] +\ 862 ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 865 ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 863 866 864 867 #Return vector of interpolated values … … 875 878 N = len(x) 876 879 assert len(y) == N, 'x and y must have same length' 877 880 878 881 res = [] 879 882 for col in q: 880 883 res.append(col*ones(N, Float)) 881 884 882 885 return res 883 886 884 887 885 888 def read_xya(filename, format = 'netcdf'): … … 889 892 890 893 Format can be either ASCII or NetCDF 891 894 892 895 Return list of points, list of attributes and 893 896 attribute names if present otherwise None 894 897 """ 895 898 #FIXME: Probably obsoleted by similar function in load_ASCII 896 899 897 900 from Scientific.IO.NetCDF import NetCDFFile 898 901 899 if format.lower() == 'netcdf': 902 if format.lower() == 'netcdf': 900 903 #Get NetCDF 901 902 fid = NetCDFFile(filename, 'r') 903 904 905 fid = NetCDFFile(filename, 'r') 906 904 907 # Get the variables 905 908 points = fid.variables['points'] … … 910 913 #Don't close - arrays are needed outside this function, 911 914 #alternatively take a copy here 912 else: 915 else: 913 916 #Read as ASCII file assuming that it is separated by whitespace 914 917 fid = open(filename) … … 927 930 attribute_names = ['elevation'] #HACK 928 931 929 attributes = {} 932 attributes = {} 930 933 for key in attribute_names: 931 934 attributes[key] = [] … … 937 940 for i, key in enumerate(attribute_names): 938 941 attributes[key].append( float(fields[2+i]) ) 939 942 940 943 return points, attributes 941 944 942 945 943 946 ##################################### … … 947 950 948 951 949 950 #Redefine these to use separate_by_polygon which will put all inside indices 952 953 #Redefine these to use separate_by_polygon which will put all inside indices 951 954 #in the first part of the indices array and outside indices in the last 952 955 … … 954 957 """See separate_points_by_polygon for documentation 955 958 """ 956 957 from Numeric import array, Float, reshape 958 959 if verbose: print 'Checking input to inside_polygon' 959 960 from Numeric import array, Float, reshape 961 962 if verbose: print 'Checking input to inside_polygon' 960 963 try: 961 points = ensure_numeric(points, Float) 964 points = ensure_numeric(points, Float) 962 965 except: 963 966 msg = 'Points could not be converted to Numeric array' 964 raise msg 965 967 raise msg 968 966 969 try: 967 polygon = ensure_numeric(polygon, Float) 970 polygon = ensure_numeric(polygon, Float) 968 971 except: 969 972 msg = 'Polygon could not be converted to Numeric array' 970 raise msg 971 972 973 973 raise msg 974 975 976 974 977 if len(points.shape) == 1: 975 978 one_point = True 976 979 points = reshape(points, (1,2)) 977 else: 980 else: 978 981 one_point = False 979 982 980 indices, count = separate_points_by_polygon(points, polygon, 981 closed, verbose) 982 983 indices, count = separate_points_by_polygon(points, polygon, 984 closed, verbose) 985 983 986 if one_point: 984 987 return count == 1 985 988 else: 986 return indices[:count] 987 989 return indices[:count] 990 988 991 def outside_polygon(points, polygon, closed = True, verbose = False): 989 992 """See separate_points_by_polygon for documentation 990 993 """ 991 992 from Numeric import array, Float, reshape 993 994 if verbose: print 'Checking input to outside_polygon' 994 995 from Numeric import array, Float, reshape 996 997 if verbose: print 'Checking input to outside_polygon' 995 998 try: 996 points = ensure_numeric(points, Float) 999 points = ensure_numeric(points, Float) 997 1000 except: 998 1001 msg = 'Points could not be converted to Numeric array' 999 raise msg 1000 1002 raise msg 1003 1001 1004 try: 1002 polygon = ensure_numeric(polygon, Float) 1005 polygon = ensure_numeric(polygon, Float) 1003 1006 except: 1004 1007 msg = 'Polygon could not be converted to Numeric array' 1005 raise msg 1006 1007 1008 1008 raise msg 1009 1010 1011 1009 1012 if len(points.shape) == 1: 1010 1013 one_point = True 1011 1014 points = reshape(points, (1,2)) 1012 else: 1015 else: 1013 1016 one_point = False 1014 1017 1015 indices, count = separate_points_by_polygon(points, polygon, 1016 closed, verbose) 1017 1018 indices, count = separate_points_by_polygon(points, polygon, 1019 closed, verbose) 1020 1018 1021 if one_point: 1019 1022 return count != 1 1020 1023 else: 1021 return indices[count:][::-1] #return reversed 1022 1023 1024 def separate_points_by_polygon(points, polygon, 1024 return indices[count:][::-1] #return reversed 1025 1026 1027 def separate_points_by_polygon(points, polygon, 1025 1028 closed = True, verbose = False): 1026 1029 """Determine whether points are inside or outside a polygon 1027 1030 1028 1031 Input: 1029 1032 points - Tuple of (x, y) coordinates, or list of tuples 1030 1033 polygon - list of vertices of polygon 1031 closed - (optional) determine whether points on boundary should be 1032 regarded as belonging to the polygon (closed = True) 1033 or not (closed = False) 1034 1034 closed - (optional) determine whether points on boundary should be 1035 regarded as belonging to the polygon (closed = True) 1036 or not (closed = False) 1037 1035 1038 Outputs: 1036 indices: array of same length as points with indices of points falling 1037 inside the polygon listed from the beginning and indices of points 1039 indices: array of same length as points with indices of points falling 1040 inside the polygon listed from the beginning and indices of points 1038 1041 falling outside listed from the end. 1039 1042 1040 1043 count: count of points falling inside the polygon 1041 1044 1042 1045 The indices of points inside are obtained as indices[:count] 1043 The indices of points outside are obtained as indices[count:] 1044 1045 1046 The indices of points outside are obtained as indices[count:] 1047 1048 1046 1049 Examples: 1047 1050 U = [[0,0], [1,0], [1,1], [0,1]] #Unit square 1048 1051 1049 1052 separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U) 1050 will return the indices [0, 2, 1] and count == 2 as only the first 1053 will return the indices [0, 2, 1] and count == 2 as only the first 1051 1054 and the last point are inside the unit square 1052 1055 1053 1056 Remarks: 1054 1057 The vertices may be listed clockwise or counterclockwise and 1055 the first point may optionally be repeated. 1058 the first point may optionally be repeated. 1056 1059 Polygons do not need to be convex. 1057 Polygons can have holes in them and points inside a hole is 1058 regarded as being outside the polygon. 1059 1060 Algorithm is based on work by Darel Finley, 1061 http://www.alienryderflex.com/polygon/ 1062 1063 """ 1060 Polygons can have holes in them and points inside a hole is 1061 regarded as being outside the polygon. 1062 1063 Algorithm is based on work by Darel Finley, 1064 http://www.alienryderflex.com/polygon/ 1065 1066 """ 1064 1067 1065 1068 from Numeric import array, Float, reshape, Int, zeros 1066 1067 1069 1070 1071 #Input checks 1072 try: 1073 points = ensure_numeric(points, Float) 1074 except: 1075 msg = 'Points could not be converted to Numeric array' 1076 raise msg 1077 1078 try: 1079 polygon = ensure_numeric(polygon, Float) 1080 except: 1081 msg = 'Polygon could not be converted to Numeric array' 1082 raise msg 1083 1084 assert len(polygon.shape) == 2,\ 1085 'Polygon array must be a 2d array of vertices' 1086 1087 assert polygon.shape[1] == 2,\ 1088 'Polygon array must have two columns' 1089 1090 assert len(points.shape) == 2,\ 1091 'Points array must be a 2d array' 1092 1093 assert points.shape[1] == 2,\ 1094 'Points array must have two columns' 1095 1096 N = polygon.shape[0] #Number of vertices in polygon 1097 M = points.shape[0] #Number of points 1098 1099 px = polygon[:,0] 1100 py = polygon[:,1] 1101 1102 #Used for an optimisation when points are far away from polygon 1103 minpx = min(px); maxpx = max(px) 1104 minpy = min(py); maxpy = max(py) 1105 1106 1107 #Begin main loop 1108 indices = zeros(M, Int) 1109 1110 inside_index = 0 #Keep track of points inside 1111 outside_index = M-1 #Keep track of points outside (starting from end) 1112 1113 for k in range(M): 1114 1115 if verbose: 1116 if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M) 1117 1118 #for each point 1119 x = points[k, 0] 1120 y = points[k, 1] 1121 1122 inside = False 1123 1124 if not x > maxpx or x < minpx or y > maxpy or y < minpy: 1125 #Check polygon 1126 for i in range(N): 1127 j = (i+1)%N 1128 1129 #Check for case where point is contained in line segment 1130 if point_on_line(x, y, px[i], py[i], px[j], py[j]): 1131 if closed: 1132 inside = True 1133 else: 1134 inside = False 1135 break 1136 else: 1137 #Check if truly inside polygon 1138 if py[i] < y and py[j] >= y or\ 1139 py[j] < y and py[i] >= y: 1140 if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x: 1141 inside = not inside 1142 1143 if inside: 1144 indices[inside_index] = k 1145 inside_index += 1 1146 else: 1147 indices[outside_index] = k 1148 outside_index -= 1 1149 1150 return indices, inside_index 1151 1152 1153 def separate_points_by_polygon_c(points, polygon, 1154 closed = True, verbose = False): 1155 """Determine whether points are inside or outside a polygon 1156 1157 C-wrapper 1158 """ 1159 1160 from Numeric import array, Float, reshape, zeros, Int 1161 1162 1163 if verbose: print 'Checking input to separate_points_by_polygon' 1068 1164 #Input checks 1069 1165 try: … … 1073 1169 raise msg 1074 1170 1171 #if verbose: print 'Checking input to separate_points_by_polygon 2' 1075 1172 try: 1076 1173 polygon = ensure_numeric(polygon, Float) 1077 1174 except: 1078 1175 msg = 'Polygon could not be converted to Numeric array' 1079 raise msg 1080 1176 raise msg 1177 1178 if verbose: print 'check' 1179 1081 1180 assert len(polygon.shape) == 2,\ 1082 1181 'Polygon array must be a 2d array of vertices' … … 1087 1186 assert len(points.shape) == 2,\ 1088 1187 'Points array must be a 2d array' 1089 1188 1090 1189 assert points.shape[1] == 2,\ 1091 1190 'Points array must have two columns' 1092 1093 N = polygon.shape[0] #Number of vertices in polygon 1191 1192 N = polygon.shape[0] #Number of vertices in polygon 1094 1193 M = points.shape[0] #Number of points 1095 1194 1096 px = polygon[:,0]1097 py = polygon[:,1]1098 1099 #Used for an optimisation when points are far away from polygon1100 minpx = min(px); maxpx = max(px)1101 minpy = min(py); maxpy = max(py)1102 1103 1104 #Begin main loop1105 indices = zeros(M, Int)1106 1107 inside_index = 0 #Keep track of points inside1108 outside_index = M-1 #Keep track of points outside (starting from end)1109 1110 for k in range(M):1111 1112 if verbose:1113 if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)1114 1115 #for each point1116 x = points[k, 0]1117 y = points[k, 1]1118 1119 inside = False1120 1121 if not x > maxpx or x < minpx or y > maxpy or y < minpy:1122 #Check polygon1123 for i in range(N):1124 j = (i+1)%N1125 1126 #Check for case where point is contained in line segment1127 if point_on_line(x, y, px[i], py[i], px[j], py[j]):1128 if closed:1129 inside = True1130 else:1131 inside = False1132 break1133 else:1134 #Check if truly inside polygon1135 if py[i] < y and py[j] >= y or\1136 py[j] < y and py[i] >= y:1137 if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:1138 inside = not inside1139 1140 if inside:1141 indices[inside_index] = k1142 inside_index += 11143 else:1144 indices[outside_index] = k1145 outside_index -= 11146 1147 return indices, inside_index1148 1149 1150 def separate_points_by_polygon_c(points, polygon,1151 closed = True, verbose = False):1152 """Determine whether points are inside or outside a polygon1153 1154 C-wrapper1155 """1156 1157 from Numeric import array, Float, reshape, zeros, Int1158 1159 1160 if verbose: print 'Checking input to separate_points_by_polygon'1161 #Input checks1162 try:1163 points = ensure_numeric(points, Float)1164 except:1165 msg = 'Points could not be converted to Numeric array'1166 raise msg1167 1168 #if verbose: print 'Checking input to separate_points_by_polygon 2'1169 try:1170 polygon = ensure_numeric(polygon, Float)1171 except:1172 msg = 'Polygon could not be converted to Numeric array'1173 raise msg1174 1175 if verbose: print 'check'1176 1177 assert len(polygon.shape) == 2,\1178 'Polygon array must be a 2d array of vertices'1179 1180 assert polygon.shape[1] == 2,\1181 'Polygon array must have two columns'1182 1183 assert len(points.shape) == 2,\1184 'Points array must be a 2d array'1185 1186 assert points.shape[1] == 2,\1187 'Points array must have two columns'1188 1189 N = polygon.shape[0] #Number of vertices in polygon1190 M = points.shape[0] #Number of points1191 1192 1195 from util_ext import separate_points_by_polygon 1193 1196 … … 1196 1199 indices = zeros( M, Int ) 1197 1200 1198 if verbose: print 'Calling C-version of inside poly' 1201 if verbose: print 'Calling C-version of inside poly' 1199 1202 count = separate_points_by_polygon(points, polygon, indices, 1200 1203 int(closed), int(verbose)) 1201 1204 1202 return indices, count 1205 return indices, count 1203 1206 1204 1207 1205 1208 1206 1209 class Polygon_function: 1207 """Create callable object f: x,y -> z, where a,y,z are vectors and 1208 where f will return different values depending on whether x,y belongs 1210 """Create callable object f: x,y -> z, where a,y,z are vectors and 1211 where f will return different values depending on whether x,y belongs 1209 1212 to specified polygons. 1210 1213 1211 1214 To instantiate: 1212 1215 1213 1216 Polygon_function(polygons) 1214 1217 1215 1218 where polygons is a list of tuples of the form 1216 1219 1217 1220 [ (P0, v0), (P1, v1), ...] 1218 1219 with Pi being lists of vertices defining polygons and vi either 1221 1222 with Pi being lists of vertices defining polygons and vi either 1220 1223 constants or functions of x,y to be applied to points with the polygon. 1221 1222 The function takes an optional argument, default which is the value 1224 1225 The function takes an optional argument, default which is the value 1223 1226 (or function) to used for points not belonging to any polygon. 1224 1227 For example: 1225 1226 Polygon_function(polygons, default = 0.03) 1227 1228 If omitted the default value will be 0.0 1229 1228 1229 Polygon_function(polygons, default = 0.03) 1230 1231 If omitted the default value will be 0.0 1232 1230 1233 Note: If two polygons overlap, the one last in the list takes precedence 1231 1234 1232 """ 1233 1235 FIXME? : Currently, coordinates specified here are assumed to be relative to the origin (georeference) used by domain. This function is more general than domain so perhaps it'd be an idea to allow the optional argument georeference??? 1236 1237 """ 1238 1234 1239 def __init__(self, regions, default = 0.0): 1235 1240 1236 1241 try: 1237 1242 len(regions) 1238 1243 except: 1239 msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons 1244 msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons 1240 1245 raise msg 1241 1246 1242 1247 1243 T = regions[0] 1248 T = regions[0] 1244 1249 try: 1245 1250 a = len(T) 1246 except: 1247 msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons 1248 raise msg 1249 1250 assert a == 2, 'Must have two component each: %s' %T 1251 1252 self.regions = regions 1251 except: 1252 msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons 1253 raise msg 1254 1255 assert a == 2, 'Must have two component each: %s' %T 1256 1257 self.regions = regions 1253 1258 self.default = default 1254 1255 1259 1260 1256 1261 def __call__(self, x, y): 1257 1262 from util import inside_polygon 1258 1263 from Numeric import ones, Float, concatenate, array, reshape, choose 1259 1264 1260 1265 x = array(x).astype(Float) 1261 y = array(y).astype(Float) 1262 1266 y = array(y).astype(Float) 1267 1263 1268 N = len(x) 1264 1269 assert len(y) == N 1265 1266 points = concatenate( (reshape(x, (N, 1)), 1270 1271 points = concatenate( (reshape(x, (N, 1)), 1267 1272 reshape(y, (N, 1))), axis=1 ) 1268 1273 1269 1274 if callable(self.default): 1270 1275 z = self.default(x,y) 1271 else: 1276 else: 1272 1277 z = ones(N, Float) * self.default 1273 1278 1274 1279 for polygon, value in self.regions: 1275 1280 indices = inside_polygon(points, polygon) 1276 1277 #FIXME: This needs to be vectorised 1281 1282 #FIXME: This needs to be vectorised 1278 1283 if callable(value): 1279 1284 for i in indices: … … 1281 1286 yy = array([y[i]]) 1282 1287 z[i] = value(xx, yy)[0] 1283 else: 1288 else: 1284 1289 for i in indices: 1285 z[i] = value 1286 1287 return z 1290 z[i] = value 1291 1292 return z 1288 1293 1289 1294 def read_polygon(filename): 1290 1295 """Read points assumed to form a polygon 1291 1296 There must be exactly two numbers in each line 1292 """ 1293 1297 """ 1298 1294 1299 #Get polygon 1295 1300 fid = open(filename) … … 1301 1306 polygon.append( [float(fields[0]), float(fields[1])] ) 1302 1307 1303 return polygon 1304 1308 return polygon 1309 1305 1310 def populate_polygon(polygon, number_of_points, seed = None): 1306 1311 """Populate given polygon with uniformly distributed points. … … 1308 1313 Input: 1309 1314 polygon - list of vertices of polygon 1310 number_of_points - (optional) number of points 1315 number_of_points - (optional) number of points 1311 1316 seed - seed for random number generator (default=None) 1312 1317 1313 1318 Output: 1314 1319 points - list of points inside polygon 1315 1320 1316 1321 Examples: 1317 1322 populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 ) 1318 will return five randomly selected points inside the unit square 1323 will return five randomly selected points inside the unit square 1319 1324 """ 1320 1325 … … 1329 1334 max_y = min_y = polygon[0][1] 1330 1335 for point in polygon[1:]: 1331 x = point[0] 1336 x = point[0] 1332 1337 if x > max_x: max_x = x 1333 if x < min_x: min_x = x 1334 y = point[1] 1338 if x < min_x: min_x = x 1339 y = point[1] 1335 1340 if y > max_y: max_y = y 1336 if y < min_y: min_y = y 1337 1338 1339 while len(points) < number_of_points: 1341 if y < min_y: min_y = y 1342 1343 1344 while len(points) < number_of_points: 1340 1345 x = uniform(min_x, max_x) 1341 y = uniform(min_y, max_y) 1346 y = uniform(min_y, max_y) 1342 1347 1343 1348 if inside_polygon( [x,y], polygon ): … … 1348 1353 def tsh2sww(filename, verbose=True): 1349 1354 """ 1350 to check if a tsh/msh file 'looks' good. 1355 to check if a tsh/msh file 'looks' good. 1351 1356 """ 1352 1357 from shallow_water import Domain 1353 1358 from pmesh2domain import pmesh_to_domain_instance 1354 import time, os 1355 from data_manager import get_dataobject 1356 from os import sep, path 1359 import time, os 1360 from data_manager import get_dataobject 1361 from os import sep, path 1357 1362 #from util import mean 1358 1363 1359 1364 if verbose == True:print 'Creating domain from', filename 1360 1365 domain = pmesh_to_domain_instance(filename, Domain) … … 1373 1378 if verbose == True: 1374 1379 print "Output written to " + domain.get_datadir() + sep + \ 1375 domain.filename + "." + domain.format 1380 domain.filename + "." + domain.format 1376 1381 sww = get_dataobject(domain) 1377 1382 sww.store_connectivity() … … 1385 1390 """ 1386 1391 """ 1387 1388 det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0) 1392 1393 det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0) 1389 1394 a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0) 1390 1395 a /= det 1391 1396 1392 1397 b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0) 1393 b /= det 1398 b /= det 1394 1399 1395 1400 return a, b … … 1411 1416 pass 1412 1417 1413 1414 1415 1416 1417 1418 #OBSOLETED STUFF 1418 1419 1420 1421 1422 1423 #OBSOLETED STUFF 1419 1424 def inside_polygon_old(point, polygon, closed = True, verbose = False): 1420 1425 #FIXME Obsoleted 1421 1426 """Determine whether points are inside or outside a polygon 1422 1427 1423 1428 Input: 1424 1429 point - Tuple of (x, y) coordinates, or list of tuples 1425 1430 polygon - list of vertices of polygon 1426 closed - (optional) determine whether points on boundary should be 1427 regarded as belonging to the polygon (closed = True) 1428 or not (closed = False) 1429 1431 closed - (optional) determine whether points on boundary should be 1432 regarded as belonging to the polygon (closed = True) 1433 or not (closed = False) 1434 1430 1435 Output: 1431 1436 If one point is considered, True or False is returned. 1432 If multiple points are passed in, the function returns indices 1433 of those points that fall inside the polygon 1434 1437 If multiple points are passed in, the function returns indices 1438 of those points that fall inside the polygon 1439 1435 1440 Examples: 1436 1441 U = [[0,0], [1,0], [1,1], [0,1]] #Unit square 1437 1442 inside_polygon( [0.5, 0.5], U) 1438 will evaluate to True as the point 0.5, 0.5 is inside the unit square 1439 1443 will evaluate to True as the point 0.5, 0.5 is inside the unit square 1444 1440 1445 inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U) 1441 will return the indices [0, 2] as only the first and the last point 1446 will return the indices [0, 2] as only the first and the last point 1442 1447 is inside the unit square 1443 1448 1444 1449 Remarks: 1445 1450 The vertices may be listed clockwise or counterclockwise and 1446 the first point may optionally be repeated. 1451 the first point may optionally be repeated. 1447 1452 Polygons do not need to be convex. 1448 Polygons can have holes in them and points inside a hole is 1449 regarded as being outside the polygon. 1450 1451 1452 Algorithm is based on work by Darel Finley, 1453 http://www.alienryderflex.com/polygon/ 1454 1455 """ 1453 Polygons can have holes in them and points inside a hole is 1454 regarded as being outside the polygon. 1455 1456 1457 Algorithm is based on work by Darel Finley, 1458 http://www.alienryderflex.com/polygon/ 1459 1460 """ 1456 1461 1457 1462 from Numeric import array, Float, reshape 1458 1459 1463 1464 1460 1465 #Input checks 1461 1466 try: 1462 point = array(point).astype(Float) 1467 point = array(point).astype(Float) 1463 1468 except: 1464 1469 msg = 'Point %s could not be converted to Numeric array' %point 1465 raise msg 1466 1470 raise msg 1471 1467 1472 try: 1468 polygon = array(polygon).astype(Float) 1473 polygon = array(polygon).astype(Float) 1469 1474 except: 1470 1475 msg = 'Polygon %s could not be converted to Numeric array' %polygon 1471 raise msg 1472 1476 raise msg 1477 1473 1478 assert len(polygon.shape) == 2,\ 1474 1479 'Polygon array must be a 2d array of vertices: %s' %polygon 1475 1480 1476 1481 1477 1482 assert polygon.shape[1] == 2,\ 1478 'Polygon array must have two columns: %s' %polygon 1483 'Polygon array must have two columns: %s' %polygon 1479 1484 1480 1485 … … 1482 1487 one_point = True 1483 1488 points = reshape(point, (1,2)) 1484 else: 1489 else: 1485 1490 one_point = False 1486 1491 points = point 1487 1492 1488 N = polygon.shape[0] #Number of vertices in polygon 1493 N = polygon.shape[0] #Number of vertices in polygon 1489 1494 M = points.shape[0] #Number of points 1490 1495 1491 1496 px = polygon[:,0] 1492 py = polygon[:,1] 1497 py = polygon[:,1] 1493 1498 1494 1499 #Used for an optimisation when points are far away from polygon 1495 1500 minpx = min(px); maxpx = max(px) 1496 minpy = min(py); maxpy = max(py) 1501 minpy = min(py); maxpy = max(py) 1497 1502 1498 1503 … … 1503 1508 if verbose: 1504 1509 if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M) 1505 1506 #for each point 1510 1511 #for each point 1507 1512 x = points[k, 0] 1508 y = points[k, 1] 1513 y = points[k, 1] 1509 1514 1510 1515 inside = False … … 1512 1517 #Optimisation 1513 1518 if x > maxpx or x < minpx: continue 1514 if y > maxpy or y < minpy: continue 1515 1516 #Check polygon 1519 if y > maxpy or y < minpy: continue 1520 1521 #Check polygon 1517 1522 for i in range(N): 1518 1523 j = (i+1)%N 1519 1520 #Check for case where point is contained in line segment 1524 1525 #Check for case where point is contained in line segment 1521 1526 if point_on_line(x, y, px[i], py[i], px[j], py[j]): 1522 1527 if closed: 1523 1528 inside = True 1524 else: 1529 else: 1525 1530 inside = False 1526 1531 break 1527 else: 1528 #Check if truly inside polygon 1532 else: 1533 #Check if truly inside polygon 1529 1534 if py[i] < y and py[j] >= y or\ 1530 1535 py[j] < y and py[i] >= y: 1531 1536 if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x: 1532 1537 inside = not inside 1533 1538 1534 1539 if inside: indices.append(k) 1535 1540 1536 if one_point: 1541 if one_point: 1537 1542 return inside 1538 1543 else: 1539 1544 return indices 1540 1545 1541 1546 1542 1547 #def remove_from(A, B): … … 1547 1552 # ind = search_sorted(A, B) 1548 1553 1549 1554 1550 1555 1551 1556 def outside_polygon_old(point, polygon, closed = True, verbose = False): 1552 1557 #OBSOLETED 1553 1558 """Determine whether points are outside a polygon 1554 1559 1555 1560 Input: 1556 1561 point - Tuple of (x, y) coordinates, or list of tuples 1557 1562 polygon - list of vertices of polygon 1558 closed - (optional) determine whether points on boundary should be 1559 regarded as belonging to the polygon (closed = True) 1560 or not (closed = False) 1561 1563 closed - (optional) determine whether points on boundary should be 1564 regarded as belonging to the polygon (closed = True) 1565 or not (closed = False) 1566 1562 1567 Output: 1563 1568 If one point is considered, True or False is returned. 1564 If multiple points are passed in, the function returns indices 1565 of those points that fall outside the polygon 1566 1569 If multiple points are passed in, the function returns indices 1570 of those points that fall outside the polygon 1571 1567 1572 Examples: 1568 U = [[0,0], [1,0], [1,1], [0,1]] #Unit square 1573 U = [[0,0], [1,0], [1,1], [0,1]] #Unit square 1569 1574 outside_polygon( [0.5, 0.5], U ) 1570 1575 will evaluate to False as the point 0.5, 0.5 is inside the unit square … … 1572 1577 ouside_polygon( [1.5, 0.5], U ) 1573 1578 will evaluate to True as the point 1.5, 0.5 is outside the unit square 1574 1579 1575 1580 outside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U ) 1576 will return the indices [1] as only the first and the last point 1581 will return the indices [1] as only the first and the last point 1577 1582 is inside the unit square 1578 1583 """ 1579 1584 1580 1585 #FIXME: This is too slow 1581 1586 1582 1587 res = inside_polygon(point, polygon, closed, verbose) 1583 1588 … … 1597 1602 1598 1603 C-wrapper 1599 """ 1604 """ 1600 1605 1601 1606 from Numeric import array, Float, reshape, zeros, Int 1602 1603 1604 if verbose: print 'Checking input to inside_polygon' 1607 1608 1609 if verbose: print 'Checking input to inside_polygon' 1605 1610 #Input checks 1606 1611 try: 1607 point = array(point).astype(Float) 1612 point = array(point).astype(Float) 1608 1613 except: 1609 1614 msg = 'Point %s could not be converted to Numeric array' %point 1610 raise msg 1611 1615 raise msg 1616 1612 1617 try: 1613 polygon = array(polygon).astype(Float) 1618 polygon = array(polygon).astype(Float) 1614 1619 except: 1615 1620 msg = 'Polygon %s could not be converted to Numeric array' %polygon 1616 raise msg 1617 1621 raise msg 1622 1618 1623 assert len(polygon.shape) == 2,\ 1619 1624 'Polygon array must be a 2d array of vertices: %s' %polygon 1620 1625 1621 1626 1622 1627 assert polygon.shape[1] == 2,\ 1623 'Polygon array must have two columns: %s' %polygon 1628 'Polygon array must have two columns: %s' %polygon 1624 1629 1625 1630 … … 1627 1632 one_point = True 1628 1633 points = reshape(point, (1,2)) 1629 else: 1634 else: 1630 1635 one_point = False 1631 1636 points = point … … 1637 1642 indices = zeros( points.shape[0], Int ) 1638 1643 1639 if verbose: print 'Calling C-version of inside poly' 1644 if verbose: print 'Calling C-version of inside poly' 1640 1645 count = inside_polygon(points, polygon, indices, 1641 1646 int(closed), int(verbose)) … … 1645 1650 else: 1646 1651 if verbose: print 'Got %d points' %count 1647 1652 1648 1653 return indices[:count] #Return those indices that were inside 1654
Note: See TracChangeset
for help on using the changeset viewer.
