Changeset 3411
- Timestamp:
- Jul 25, 2006, 12:14:42 PM (19 years ago)
- Location:
- inundation/fit_interpolate
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/fit_interpolate/interpolate.py
r3409 r3411 21 21 import os 22 22 from warnings import warn 23 from math import sqrt 24 from csv import writer, DictWriter 23 25 24 26 from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \ … … 33 35 from utilities.numerical_tools import ensure_numeric, mean, INF 34 36 from utilities.polygon import in_and_outside_polygon 35 from geospatial_data.geospatial_data import Geospatial_data 37 from geospatial_data.geospatial_data import Geospatial_data, ensure_absolute 36 38 from fit_interpolate.search_functions import search_tree_of_vertices 37 39 from fit_interpolate.general_fit_interpolate import FitInterpolate 38 40 from pyvolution.util import file_function 39 41 40 42 … … 276 278 return A 277 279 280 def interpolate_sww2csv(sww_file, points, 281 depth_file, 282 velocity_file, 283 #quantities = ['depth', 'velocity'], 284 verbose=True, 285 use_cache = True): 286 """ 287 Interpolate the quantities at a given set of locations, given 288 an sww file. 289 The results are written to a csv file. 290 291 In the future let points be a csv or xya file. 292 And the user choose the quantities. 293 """ 294 quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] 295 #print "points",points 296 points = ensure_absolute(points) 297 point_count = len(points) 298 callable_sww = file_function(sww_file, 299 quantities=quantities, 300 interpolation_points=points, 301 verbose=verbose, 302 use_cache=use_cache) 303 304 depth_writer = writer(file(depth_file, "wb")) 305 velocity_writer = writer(file(velocity_file, "wb")) 306 # Write heading 307 depth_writer.writerow([str(x[0])+ ':' + str(x[1]) for x in points]) 308 velocity_writer.writerow([str(x[0])+ ':' + str(x[1]) for x in points]) 309 #for row in someiterable: 310 #csvwriter.writerow(row) 311 for time in callable_sww.get_time(): 312 depths = [] 313 velocitys = [] 314 for point_i, point in enumerate(points): 315 quantities = callable_sww(time,point_i) 316 #print "quantities", quantities 317 318 w = quantities[0] 319 z = quantities[1] 320 uh = quantities[2] 321 vh = quantities[3] 322 depth = w - z 323 324 if w == INF or z == INF or uh == INF or vh == INF: 325 velocity = INF 326 else: 327 momentum = sqrt(uh*uh + vh*vh) 328 if depth > 1.e-30: # use epsilon 329 velocity = momentum / depth #Absolute velocity 330 else: 331 velocity = 0 332 333 print "depth",depth 334 print "velocity",velocity 335 depths.append(depth) 336 velocitys.append(velocity) 337 depth_writer.writerow(depths) 338 velocity_writer.writerow(velocitys) 339 278 340 class Interpolation_function: 279 341 """Interpolation_interface - creates callable object f(t, id) or f(t,x,y) … … 715 777 #------------------------------------------------------------- 716 778 if __name__ == "__main__": 717 a = [0.0, 0.0] 718 b = [0.0, 2.0] 719 c = [2.0,0.0] 720 points = [a, b, c] 721 vertices = [ [1,0,2] ] #bac 722 723 data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point 724 725 interp = Interpolate(points, vertices) #, data) 726 A = interp._build_interpolation_matrix_A(data, verbose=True) 727 A = A.todense() 728 print "A",A 729 assert allclose(A, [[1./3, 1./3, 1./3]]) 730 print "finished" 779 names = ["x","y"] 780 someiterable = [[1,2],[3,4]] 781 csvwriter = writer(file("some.csv", "wb")) 782 csvwriter.writerow(names) 783 for row in someiterable: 784 csvwriter.writerow(row) -
inundation/fit_interpolate/test_interpolate.py
r3019 r3411 23 23 from data_manager import get_dataobject 24 24 from geospatial_data.geospatial_data import Geospatial_data 25 from pmesh.mesh import Mesh 25 26 26 27 def distance(x, y): … … 1351 1352 self.failUnless( z[i,1] == answer[11,1], 'Fail!') 1352 1353 self.failUnless( z[i,0] == answer[11,0], 'Fail!') 1353 1354 1355 def ztest_interpolate_sww2csv(self): 1356 1357 def elevation_function(x, y): 1358 return -x 1359 1360 # create mesh 1361 mesh_file = tempfile.mktemp(".tsh") 1362 points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]] 1363 m = Mesh() 1364 m.add_vertices(points) 1365 m.autoSegment() 1366 m.generate_mesh(verbose=False) 1367 m.export_mesh_file(mesh_file) 1368 1369 #Create shallow water domain 1370 domain = Domain(mesh_file) 1371 os.remove(mesh_file) 1372 1373 domain.default_order=2 1374 domain.beta_h = 0 1375 1376 #Set some field values 1377 domain.set_quantity('elevation', elevation_function) 1378 domain.set_quantity('friction', 0.03) 1379 domain.set_quantity('xmomentum', 10.0) 1380 domain.set_quantity('ymomentum', 100.0) 1381 1382 ###################### 1383 # Boundary conditions 1384 B = Transmissive_boundary(domain) 1385 domain.set_boundary( {'exterior': B}) 1386 1387 # This call mangles the stage values. 1388 domain.distribute_to_vertices_and_edges() 1389 domain.set_quantity('stage', 1.0) 1390 1391 #sww_file = tempfile.mktemp("") 1392 domain.filename = 'datatest' + str(time.time()) 1393 #domain.filename = sww_file 1394 #print "domain.filename",domain.filename 1395 domain.format = 'sww' 1396 domain.smooth = True 1397 domain.reduction = mean 1398 1399 sww = get_dataobject(domain) 1400 sww.store_connectivity() 1401 sww.store_timestep(['stage', 'xmomentum', 'ymomentum']) 1402 domain.set_quantity('stage', -1.0) 1403 domain.time = 2. 1404 sww.store_timestep(['stage', 'xmomentum', 'ymomentum']) 1405 1406 # test the function 1407 points = [[5.0,1.],[0.5,2.]] 1408 interpolate_sww2csv(sww.filename, points, "aa.csv", "bb.csv") 1409 1410 print "sww.filename",sww.filename 1411 #os.remove(sww.filename) 1354 1412 #------------------------------------------------------------- 1355 1413 if __name__ == "__main__": 1356 1414 1357 suite = unittest.makeSuite(Test_Interpolate,'test')1358 #suite = unittest.makeSuite(Test_Interpolate,'test_interpolation_precompute_points')1415 #suite = unittest.makeSuite(Test_Interpolate,'test') 1416 suite = unittest.makeSuite(Test_Interpolate,'test_interpolate_sww2csv') 1359 1417 runner = unittest.TextTestRunner(verbosity=1) 1360 1418 runner.run(suite)
Note: See TracChangeset
for help on using the changeset viewer.