Changeset 1668
 Timestamp:
 Aug 1, 2005, 5:30:53 PM (19 years ago)
 Location:
 inundation/ga/storm_surge/pyvolution
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

inundation/ga/storm_surge/pyvolution/test_util.py
r1664 r1668 559 559 interpolation_points = interpolation_points) 560 560 assert allclose(domain.starttime, start+1) 561 562 561 domain.starttime = start 563 562 564 563 565 566 564 #Check linear interpolation in time 565 F = file_function(filename + '.sww', domain, 566 quantities = domain.conserved_quantities, 567 interpolation_points = interpolation_points) 567 568 #for id in range(len(interpolation_points)): 568 569 for id in [1]: … … 580 581 581 582 q = F(t, point_id=id) 583 #print i, k, t, q 584 #print ' ', q0 585 #print ' ', q1 586 #print 's', (k*q1 + (6k)*q0)/6 587 #print 582 588 assert allclose(q, (k*q1 + (6k)*q0)/6) 583 589 … … 589 595 590 596 t = 90 #Halfway between 60 and 120 591 q = F(t, point_id=id)597 q = F(t, point_id=id) 592 598 assert allclose( (q120+q60)/2, q ) 593 599 … … 1083 1089 # 1084 1090 if __name__ == "__main__": 1085 #suite = unittest.makeSuite(Test Case,'test_inside_polygon_main')1091 #suite = unittest.makeSuite(Test_Util,'test_spatio_temporal_file_function_time') 1086 1092 suite = unittest.makeSuite(Test_Util,'test') 1087 1093 runner = unittest.TextTestRunner() 
inundation/ga/storm_surge/pyvolution/util.py
r1664 r1668 5 5 """ 6 6 7 #FIXME: Move polygon stuff out 8 #FIXME: Finish Interpolation function and get rid of File_function_ASCII 9 10 7 11 def angle(v): 8 12 """Compute angle between e1 (the unit vector in the xdirection) … … 17 21 v2 = v[1]/l 18 22 23 19 24 theta = acos(v1) 20 try: 21 theta = acos(v1) 22 except ValueError, e: 23 print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e) 24 25 #FIXME, hack to avoid acos(1.0) Value error 26 # why is it happening? 27 # is it catching something we should avoid? 28 #FIXME (Ole): When did this happen? We need a unit test to 29 #reveal this condition or otherwise remove the hack 30 31 s = 1e6 32 if (v1+s > 1.0) and (v1s < 1.0) : 33 theta = 0.0 34 elif (v1+s > 1.0) and (v1s < 1.0): 35 theta = 3.1415926535897931 36 print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\ 37 %(v1, theta) 25 26 #try: 27 # theta = acos(v1) 28 #except ValueError, e: 29 # print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e) 30 # 31 # #FIXME, hack to avoid acos(1.0) Value error 32 # # why is it happening? 33 # # is it catching something we should avoid? 34 # #FIXME (Ole): When did this happen? We need a unit test to 35 # #reveal this condition or otherwise remove the hack 36 # 37 # s = 1e6 38 # if (v1+s > 1.0) and (v1s < 1.0) : 39 # theta = 0.0 40 # elif (v1+s > 1.0) and (v1s < 1.0): 41 # theta = 3.1415926535897931 42 # print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\ 43 # %(v1, theta) 38 44 39 45 if v2 < 0: … … 156 162 quantities = None 157 163 158 #Choose format 159 #FIXME: Maybe these can be merged later on 160 #FIXME: We shoulp probably not waste energy on ASCII files  161 #the different formats are a pain. E.g. the time field. 162 # 164 163 165 if line[:3] == 'CDF': 164 return File_function_NetCDF(filename, domain, quantities,165 interpolation_points, verbose = verbose)166 return get_netcdf_file_function(filename, domain, quantities, 167 interpolation_points, verbose = verbose) 166 168 else: 169 #FIXME Should be phased out 167 170 return File_function_ASCII(filename, domain, quantities, 168 171 interpolation_points) 169 170 172 173 174 175 176 177 178 179 180 181 182 def get_netcdf_file_function(filename, domain=None, quantity_names=None, 183 interpolation_points=None, verbose = False): 184 """Read time history of spatial data from NetCDF sww file and 185 return a callable object f(t,x,y) 186 which will return interpolated values based on the input file. 187 188 If domain is specified, model time (domain.starttime) 189 will be checked and possibly modified 190 191 All times are assumed to be in UTC 192 193 See Interpolation function for further documetation 194 195 196 """ 197 198 #verbose = True 199 200 #FIXME: Check that model origin is the same as file's origin 201 #(both in UTM coordinates) 202 #If not  modify those from file to match domain 203 #Take this code from e.g. dem2pts in data_manager.py 204 #FIXME: Use geo_reference to read and write xllcorner... 205 206 207 import time, calendar 208 from config import time_format 209 from Scientific.IO.NetCDF import NetCDFFile 210 from Numeric import array, zeros, Float, alltrue, concatenate, reshape 211 212 #Open NetCDF file 213 if verbose: print 'Reading', filename 214 fid = NetCDFFile(filename, 'r') 215 216 if quantity_names is None or len(quantity_names) < 1: 217 msg = 'ERROR: At least one independent value must be specified' 218 raise msg 219 220 missing = [] 221 for quantity in ['x', 'y', 'time'] + quantity_names: 222 if not fid.variables.has_key(quantity): 223 missing.append(quantity) 224 225 if len(missing) > 0: 226 msg = 'Quantities %s could not be found in file %s'\ 227 %(str(missing), filename) 228 raise msg 229 230 231 #Get first timestep 232 try: 233 starttime = fid.starttime[0] 234 except ValueError: 235 msg = 'Could not read starttime from file %s' %filename 236 raise msg 237 238 #Get origin 239 xllcorner = fid.xllcorner[0] 240 yllcorner = fid.yllcorner[0] 241 242 243 #Get variables 244 if verbose: print 'Get variables' 245 x = fid.variables['x'][:] 246 y = fid.variables['y'][:] 247 triangles = fid.variables['volumes'][:] 248 time = fid.variables['time'][:] 249 250 x = reshape(x, (len(x),1)) 251 y = reshape(y, (len(y),1)) 252 vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array 253 254 if domain is not None: 255 #Update domain.startime if it is *earlier* than starttime 256 if starttime > domain.starttime: 257 msg = 'WARNING: Start time as specified in domain (%f)'\ 258 %domain.starttime 259 msg += ' is earlier than the starttime of file %s (%f).'\ 260 %(filename, starttime) 261 msg += ' Modifying domain starttime accordingly.' 262 263 if verbose: print msg 264 domain.starttime = starttime #Modifying model time 265 if verbose: print 'Domain starttime is now set to %f'\ 266 %domain.starttime 267 268 269 #If domain.startime is *later* than starttime, 270 #move time back  relative to domain's time 271 if domain.starttime > starttime: 272 time = time  domain.starttime + starttime 273 274 275 if verbose: 276 print 'File_function data obtained from: %s' %filename 277 print ' References:' 278 #print ' Datum ....' #FIXME 279 print ' Lower left corner: [%f, %f]'\ 280 %(xllcorner, yllcorner) 281 print ' Start time: %f' %starttime 282 283 284 #Produce values for desired data points at 285 #each timestep 286 287 quantities = {} 288 for i, name in enumerate(quantity_names): 289 quantities[name] = fid.variables[name][:] 290 fid.close() 291 292 293 return Interpolation_function(vertex_coordinates, 294 triangles, 295 time, 296 quantities, 297 quantity_names, 298 interpolation_points, 299 alpha = 0, 300 verbose = verbose) 301 302 303 304 305 306 class Interpolation_function: 307 """Interpolation_function  creates callable object f(t, id) or f(t,x,y) 308 which is interpolated from time series defined at vertices of 309 triangular mesh 310 311 Input 312 vertex_coordinates: mx2 array of coordinates (Float) 313 triangles: nx3 array of indices into vertex_coordinates (Int) 314 time: px1 array of times (Float) 315 quantities: Dictionary of pxm values 316 quantity_names: 317 interpolation_points: 318 alpha: 319 verbose: 320 321 322 323 The quantities returned by the callable object are specified by 324 the list quantities which must contain the names of the 325 quantities to be returned and also reflect the order, e.g. for 326 the shallow water wave equation, on would have 327 quantities = ['stage', 'xmomentum', 'ymomentum'] 328 329 The parameter interpolation_points decides at which points interpolated 330 quantities are to be computed whenever object is called. 331 If None, return average value 332 """ 333 334 335 336 337 def __init__(self, 338 vertex_coordinates, 339 triangles, 340 time, 341 quantities, 342 quantity_names = None, 343 interpolation_points = None, 344 alpha = 0, 345 verbose = False): 346 """ 347 """ 348 349 from Numeric import array, zeros, Float, alltrue, concatenate,\ 350 reshape, ArrayType 351 352 from config import time_format 353 from least_squares import Interpolation 354 355 #Check 356 msg = 'Time must be a monotonuosly ' 357 msg += 'increasing sequence %s' %time 358 assert alltrue(time[1:]  time[:1] > 0 ), msg 359 360 361 #Check if quantities is a single array only 362 if type(quantities) == ArrayType: 363 quantity_names = 'Attribute' 364 quantities = {quantity_names: quantities} 365 366 #Use keys if no names specified 367 if quantity_names is not None: 368 self.quantity_names = quantity_names 369 else: 370 self.quantity_names = quantities.keys() 371 372 373 self.interpolation_points = interpolation_points 374 self.T = time[:] #Time assumed to be relative to starttime 375 self.index = 0 #Initial time index 376 self.values = {} 377 378 379 if interpolation_points is not None: 380 #Spatial interpolation 381 382 try: 383 interpolation_points = ensure_numeric(interpolation_points) 384 except: 385 msg = 'Interpolation points must be an N x 2 Numeric array or a list of points\n' 386 msg += 'I got: %s.' %( str(interpolation_points)[:60] + '...') 387 raise msg 388 389 390 for name in quantity_names: 391 self.values[name] = zeros( (len(self.T), 392 len(interpolation_points)), 393 Float) 394 395 #Build interpolator 396 interpol = Interpolation(vertex_coordinates, 397 triangles, 398 point_coordinates = interpolation_points, 399 alpha = 0, 400 verbose = verbose) 401 402 403 if verbose: print 'Interpolate' 404 for i, t in enumerate(self.T): 405 #Interpolate quantities at this timestep 406 if verbose: print ' time step %d of %d' %(i, len(self.T)) 407 for name in quantity_names: 408 self.values[name][i, :] =\ 409 interpol.interpolate(quantities[name][i,:]) 410 411 #Report 412 if verbose: 413 x = vertex_coordinates[:,0] 414 y = vertex_coordinates[:,1] 415 416 print '' 417 print 'Interpolation_function statistics:' 418 print ' Extent:' 419 print ' x in [%f, %f], len(x) == %d'\ 420 %(min(x.flat), max(x.flat), len(x.flat)) 421 print ' y in [%f, %f], len(y) == %d'\ 422 %(min(y.flat), max(y.flat), len(y.flat)) 423 print ' t in [%f, %f], len(t) == %d'\ 424 %(min(self.T), max(self.T), len(self.T)) 425 print ' Quantities:' 426 for name in quantity_names: 427 q = quantities[name][:].flat 428 print ' %s in [%f, %f]' %(name, min(q), max(q)) 429 430 431 print ' Interpolation points (xi, eta):'\ 432 'number of points == %d ' %interpolation_points.shape[0] 433 print ' xi in [%f, %f]' %(min(interpolation_points[:,0]), 434 max(interpolation_points[:,0])) 435 print ' eta in [%f, %f]' %(min(interpolation_points[:,1]), 436 max(interpolation_points[:,1])) 437 print ' Interpolated quantities (over all timesteps):' 438 for name in quantity_names: 439 q = self.values[name][:].flat 440 print ' %s at interpolation points in [%f, %f]'\ 441 %(name, min(q), max(q)) 442 print '' 443 else: 444 #Return an average, making this a time series 445 #FIXME: Needs a test 446 for name in quantity_names: 447 self.values[name] = avg(quantities[name][i,:]) 448 449 450 def __repr__(self): 451 return 'Interpolation function' 452 453 def __call__(self, t, point_id = None, x = None, y = None): 454 """Evaluate f(t), f(t, point_id) or f(t, x, y) 455 456 Inputs: 457 t: time  Model time. Must lie within existing timesteps 458 point_id: index of one of the preprocessed points. 459 If point_id is None all preprocessed points are computed 460 461 If no x,y info is present, point_id and x,y arguments are ignored 462 making f a function of time only. 463 464 FIXME: point_id could also be a slice 465 FIXME: One could allow arbitrary x, y coordinates 466 (requires computation of a new interpolator) 467 FIXME: What if x and y are vectors? 468 FIXME: Does point_id or x,y take precedence? 469 FIXME: What about f(x,y) without t 470 471 """ 472 473 from math import pi, cos, sin, sqrt 474 from Numeric import zeros, Float 475 476 477 msg = 'Time interval [%s:%s]' %(self.T[0], self.T[1]) 478 msg += ' does not match model time: %s\n' %t 479 if t < self.T[0]: raise msg 480 if t > self.T[1]: raise msg 481 482 oldindex = self.index #Time index 483 484 #Find current time slot 485 while t > self.T[self.index]: self.index += 1 486 while t < self.T[self.index]: self.index = 1 487 488 if t == self.T[self.index]: 489 #Protect against case where t == T[1] (last time) 490 #  also works in general when t == T[i] 491 ratio = 0 492 else: 493 #t is now between index and index+1 494 ratio = (t  self.T[self.index])/\ 495 (self.T[self.index+1]  self.T[self.index]) 496 497 #Compute interpolated values 498 q = zeros(len(self.quantity_names), Float) 499 500 for i, name in enumerate(self.quantity_names): 501 Q = self.values[name] 502 503 if point_id is None: 504 Q0 = Q[self.index] 505 506 try: 507 Q1 = Q[self.index+1] 508 except: 509 pass 510 else: 511 Q0 = Q[self.index, point_id] 512 513 try: 514 Q1 = Q[self.index+1, point_id] 515 except: 516 pass 517 518 519 #Linear temporal interpolation 520 if ratio > 0: 521 q[i] = Q0 + ratio*(Q1  Q0) 522 else: 523 q[i] = Q0 524 525 526 #if ratio > 0: 527 # q[i] = Q[self.index, point_id] +\ 528 # ratio*(Q[self.index+1, point_id] \ 529 # Q[self.index, point_id]) 530 #else: 531 # q[i] = Q[self.index, point_id] 532 533 #Return vector of interpolated values 534 return q 535 536 537 538 539 540 541 542 543 544 545 #FIXME Obsolete 171 546 class File_function_NetCDF: 172 547 """Read time history of spatial data from NetCDF sww file and … … 203 578 #If not  modify those from file to match domain 204 579 #Take this code from e.g. dem2pts in data_manager.py 580 #FIXME: Use geo_reference to read and write xllcorner... 205 581 206 582 … … 281 657 x = fid.variables['x'][:] 282 658 y = fid.variables['y'][:] 283 z = fid.variables['elevation'][:]659 #z = fid.variables['elevation'][:] 284 660 triangles = fid.variables['volumes'][:] 285 661 time = fid.variables['time'][:] … … 410 786 tau = self.domain.starttimeself.starttime+t 411 787 else: 412 #print 'DOMAIN IS NONE!!!!!!!!!'413 788 tau = t 414 789 415 #print 'D start', self.domain.starttime416 #print 'S start', self.starttime417 #print 't', t418 #print 'tau', tau419 790 420 791 msg = 'Time interval derived from file %s [%s:%s]'\ … … 464 835 return q 465 836 466 837 #FIXME Obsolete 467 838 class File_function_ASCII: 468 839 """Read time series from file and return a callable object:
Note: See TracChangeset
for help on using the changeset viewer.