Changeset 656
- Timestamp:
- Dec 2, 2004, 12:47:16 PM (20 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/test_util.py
r655 r656 332 332 assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] ) 333 333 334 335 336 def test_polygon_function(self): 337 p1 = [[0,0], [10,0], [10,10], [0,10]] 338 p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]] 339 340 f = Polygon_function( [(p1, 1.0)] ) 341 z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1 342 assert allclose(z, [1,1,0,0]) 343 344 345 f = Polygon_function( [(p2, 2.0)] ) 346 z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2 347 assert allclose(z, [2,0,0,2]) 348 349 350 #Combined 351 f = Polygon_function( [(p1, 1.0), (p2, 2.0)] ) 352 z = f([5, 5, 27, 35], [5, 9, 8, -5]) 353 assert allclose(z, [2,1,0,2]) 354 355 356 334 357 def test_point_on_line(self): 335 358 … … 355 378 polygon = [[0,0], [1,0], [1,1], [0,1]] 356 379 357 x = y = 0.5 358 assert inside_polygon( (x, y), polygon ) 359 360 y = 1.5 361 assert not inside_polygon( (x, y), polygon ) 380 assert inside_polygon( (0.5, 0.5), polygon ) 381 assert not inside_polygon( (0.5, 1.5), polygon ) 382 assert not inside_polygon( (0.5, -0.5), polygon ) 383 assert not inside_polygon( (-0.5, 0.5), polygon ) 384 assert not inside_polygon( (1.5, 0.5), polygon ) 385 362 386 #Try point on borders 363 387 assert inside_polygon( (1., 0.5), polygon, closed=True) … … 381 405 assert not inside_polygon( (0.5, -0.5), polygon ) 382 406 383 384 385 407 #Now try the vector formulation returning indices 386 408 points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] … … 389 411 assert allclose( res, [0,1,2] ) 390 412 413 #Very convoluted polygon 414 polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]] 415 assert inside_polygon( (5, 5), polygon ) 416 assert inside_polygon( (17, 7), polygon ) 417 assert inside_polygon( (27, 2), polygon ) 418 assert inside_polygon( (35, -5), polygon ) 419 assert not inside_polygon( (15, 7), polygon ) 420 assert not inside_polygon( (24, 3), polygon ) 421 assert not inside_polygon( (25, -10), polygon ) 422 423 424 425 #Another combination (that failed) 426 polygon = [[0,0], [10,0], [10,10], [0,10]] 427 assert inside_polygon( (5, 5), polygon ) 428 assert inside_polygon( (7, 7), polygon ) 429 assert not inside_polygon( (-17, 7), polygon ) 430 assert not inside_polygon( (7, 17), polygon ) 431 assert not inside_polygon( (17, 7), polygon ) 432 assert not inside_polygon( (27, 8), polygon ) 433 assert not inside_polygon( (35, -5), polygon ) 434 435 436 437 438 #Multiple polygons 439 440 polygon = [[0,0], [1,0], [1,1], [0,1], [0,0], 441 [10,10], [11,10], [11,11], [10,11], [10,10]] 442 assert inside_polygon( (0.5, 0.5), polygon ) 443 assert inside_polygon( (10.5, 10.5), polygon ) 444 445 #FIXME: Fails 446 #assert not inside_polygon( (5.5, 5.5), polygon ) 447 448 #Polygon with a hole 449 polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1], 450 [0,0], [1,0], [1,1], [0,1], [0,0]] 451 452 assert inside_polygon( (0, -0.5), polygon ) 453 454 #FIXME: Fails 455 #assert not inside_polygon( (0.5, 0.5), polygon ) 456 457 points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] 391 458 392 459 -
inundation/ga/storm_surge/pyvolution/util.py
r655 r656 154 154 px = polygon[:,0] 155 155 py = polygon[:,1] 156 156 157 157 158 #Begin algorithm 158 159 indices = [] … … 161 162 x = points[k, 0] 162 163 y = points[k, 1] 163 164 164 165 165 inside = False 166 166 for i in range(N): … … 174 174 inside = False 175 175 break 176 177 #Checktruly inside polygon178 if px[i] < y and py[j] >= y or\179 py[j] < y and py[i] >= y:180 if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:181 inside = not inside176 else: 177 #Check if truly inside polygon 178 if py[i] < y and py[j] >= y or\ 179 py[j] < y and py[i] >= y: 180 if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x: 181 inside = not inside 182 182 183 183 if inside: indices.append(k) 184 185 #print x, y, polygon, inside 186 184 187 185 188 if one_point: … … 189 192 190 193 194 class Polygon_function: 195 """Create callable object f: x,y -> z, where a,y,z are vectors and 196 where f will return different values depending on whether x,y belongs 197 to specified polygons. 198 199 To instantiate: 200 201 Polygon_function(polygons) 202 203 where polygons is a dictionary of the form 204 205 {P0: v0, P1: v1, ...} 206 207 with Pi being lists of vertices defining polygons and vi either 208 constants or functions of x,y to be applied to points with the polygon. 209 210 The function takes an optional argument, default which is the value 211 (or function) to used for points not belonging to any polygon. 212 For example: 213 214 Polygon_function(polygons, default = 0.03) 215 216 If omitted the default value will be 0.0 217 218 Note: If two polygons overlap, the one last in the list takes precedence 219 220 Note: default can only be a constant. 221 """ 222 223 def __init__(self, regions, default = 0.0): 224 225 try: 226 len(regions) 227 except: 228 msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons 229 raise msg 230 231 232 T = regions[0] 233 try: 234 a = len(T) 235 except: 236 msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons 237 raise msg 238 239 assert a == 2, 'Must have two component each: %s' %T 240 241 self.regions = regions 242 self.default = default 243 244 245 def __call__(self, x, y): 246 from util import inside_polygon 247 from Numeric import ones, Float, concatenate, array, reshape, choose 248 249 x = array(x).astype(Float) 250 y = array(y).astype(Float) 251 x = reshape(x, (len(x), 1)) 252 y = reshape(y, (len(y), 1)) 253 254 points = concatenate( (x,y), axis=1 ) 255 256 z = ones(x.shape[0], Float) * self.default 257 for polygon, value in self.regions: 258 indices = inside_polygon(points, polygon) 259 260 for i in indices: 261 z[i] = value 262 263 return z 264 191 265 class File_function: 192 266 """Read time series from file and return a callable object:
Note: See TracChangeset
for help on using the changeset viewer.