Changeset 3411


Ignore:
Timestamp:
Jul 25, 2006, 12:14:42 PM (19 years ago)
Author:
duncan
Message:

adding sww to csv function

Location:
inundation/fit_interpolate
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • inundation/fit_interpolate/interpolate.py

    r3409 r3411  
    2121import os
    2222from warnings import warn
     23from math import sqrt
     24from csv import writer, DictWriter
    2325
    2426from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
     
    3335from utilities.numerical_tools import ensure_numeric, mean, INF
    3436from utilities.polygon import in_and_outside_polygon
    35 from geospatial_data.geospatial_data import Geospatial_data
     37from geospatial_data.geospatial_data import Geospatial_data, ensure_absolute
    3638from fit_interpolate.search_functions import search_tree_of_vertices
    3739from fit_interpolate.general_fit_interpolate import FitInterpolate
    38 
     40from pyvolution.util import file_function
    3941
    4042
     
    276278        return A
    277279
     280def 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
    278340class Interpolation_function:
    279341    """Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
     
    715777#-------------------------------------------------------------
    716778if __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  
    2323from data_manager import get_dataobject
    2424from geospatial_data.geospatial_data import Geospatial_data
     25from pmesh.mesh import Mesh
    2526
    2627def distance(x, y):
     
    13511352            self.failUnless( z[i,1] == answer[11,1], 'Fail!')
    13521353            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)
    13541412#-------------------------------------------------------------
    13551413if __name__ == "__main__":
    13561414
    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')
    13591417    runner = unittest.TextTestRunner(verbosity=1)
    13601418    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.