Changeset 4503


Ignore:
Timestamp:
May 28, 2007, 4:34:31 PM (18 years ago)
Author:
duncan
Message:

inundation_damage now inputs a base sww name, then iterates over all sww files

Location:
anuga_core/source/anuga/damage_modelling
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/damage_modelling/inundation_damage.py

    r4326 r4503  
    44   Geoscience Australia, 2006
    55"""
    6 
     6import os
    77from math import sqrt
    88from Scientific.Functions.Interpolation import InterpolatingFunction
    … …  
    4141CONT_VALUE_LABEL = 'CONT_VALUE'
    4242
    43 def inundation_damage(sww_file, exposure_file_in,
     43def inundation_damage(sww_base_name, exposure_file_in,
    4444                      exposure_file_out=None,
    4545                      ground_floor_height=0.3,
    … …  
    5656    Note, structures outside of the sww file get the minimum inundation
    5757    (-ground_floor_height).
     58   
     59    These calculations are done over all the sww files with the sww_base_name
     60    in the specified directory.
    5861    """
    5962
    … …  
    6366    geospatial = csv.get_location()
    6467    geospatial = ensure_absolute(geospatial)
    65     max_depths, max_momentums = calc_max_depth_and_momentum(sww_file,
     68    max_depths, max_momentums = calc_max_depth_and_momentum(sww_base_name,
    6669                                    geospatial,
    6770                                    ground_floor_height=ground_floor_height,
    … …  
    8386    csv.save(exposure_file_out)
    8487   
    85 def add_depth_and_momentum2csv(sww_file, exposure_file_in,
     88def add_depth_and_momentum2csv(sww_base_name, exposure_file_in,
    8689                      exposure_file_out=None,
    8790                      overwrite=False, verbose=True,
    … …  
    9093    Calculate the maximum depth and momemtum in an sww file, for locations
    9194    specified in an csv exposure file.
     95   
     96    These calculations are done over all the sww files with the sww_base_name
     97    in the specified directory.
    9298    """
    9399
    94100    csv = Exposure_csv(exposure_file_in)
    95101    geospatial = csv.get_location()
    96     max_depths, max_momentums = calc_max_depth_and_momentum(sww_file,
     102    max_depths, max_momentums = calc_max_depth_and_momentum(sww_base_name,
    97103                                                          geospatial,
    98104                                                          verbose=verbose,
    … …  
    102108    csv.save(exposure_file_out)
    103109   
    104 def calc_max_depth_and_momentum(sww_file, points,
     110def calc_max_depth_and_momentum(sww_base_name, points,
    105111                                ground_floor_height=0.0,
    106112                                verbose=True,
    … …  
    112118    The inundation value is in the range -ground_floor_height to
    113119    overflow errors.
     120
     121    These calculations are done over all the sww files with the sww_base_name
     122    in the specified directory.
    114123    """
    115124    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
    … …  
    117126    points = ensure_absolute(points)
    118127    point_count = len(points)
    119     callable_sww = file_function(sww_file,
    120                                  quantities=quantities,
    121                                  interpolation_points=points,
    122                                  verbose=verbose,
    123                                  use_cache=use_cache)
     128
    124129    # initialise the max lists
    125130    max_depths = [-ground_floor_height]*point_count
    126131    max_momentums = [-ground_floor_height]*point_count
    127     for point_i, point in enumerate(points):
    128         for time in callable_sww.get_time():
    129             quantities = callable_sww(time,point_i)
    130             #print "quantities", quantities
    131            
    132             w = quantities[0]
    133             z = quantities[1]
    134             uh = quantities[2]
    135             vh = quantities[3]
    136 
    137             #  -ground_floor_height is the minimum value.
    138             depth = w - z - ground_floor_height
     132   
     133    # How many sww files are there?
     134    dir, base = os.path.split(sww_base_name)
     135    #print "basename_in",basename_in
     136    if base[-4:] == '.sww':
     137        base = base[:-4]
     138    #print "base",base
     139    if dir == "": dir = "." # Unix compatibility
     140    dir_ls = os.listdir(dir)
     141    interate_over = [x for x in dir_ls if base in x and x[-4:] == '.sww']
     142    #print "interate_over", interate_over
     143
     144    for this_sww_file in interate_over:
     145        callable_sww = file_function(this_sww_file,
     146                                     quantities=quantities,
     147                                     interpolation_points=points,
     148                                     verbose=verbose,
     149                                     use_cache=use_cache)
     150   
     151        for point_i, point in enumerate(points):
     152            for time in callable_sww.get_time():
     153                quantity_values = callable_sww(time,point_i)
     154           
     155                w = quantity_values[0]
     156                z = quantity_values[1]
     157                uh = quantity_values[2]
     158                vh = quantity_values[3]
     159               
     160                #  -ground_floor_height is the minimum value.
     161                depth = w - z - ground_floor_height
    139162             
    140             if depth > max_depths[point_i]:
    141                 max_depths[point_i] = depth
    142             if w == NAN or z == NAN or uh == NAN or vh == NAN:
    143                 continue
    144             momentum = sqrt(uh*uh + vh*vh)
    145             if momentum > max_momentums[point_i]:
    146                 max_momentums[point_i] = momentum
     163                if depth > max_depths[point_i]:
     164                    max_depths[point_i] = depth
     165                if w == NAN or z == NAN or uh == NAN or vh == NAN:
     166                    continue
     167                momentum = sqrt(uh*uh + vh*vh)
     168                if momentum > max_momentums[point_i]:
     169                    max_momentums[point_i] = momentum
    147170    return max_depths, max_momentums
    148171

  • anuga_core/source/anuga/damage_modelling/test_inundation_damage.py

    r4324 r4503  
    44import unittest
    55import tempfile
    6 import os
     6import os, sys
    77import time
    88import csv
    … …  
    2424
    2525class Test_inundation_damage(unittest.TestCase):
     26    # Class variable
     27    verbose = False
     28
     29    def set_verbose(self):
     30        Test_Data_Manager.verbose = True
     31       
    2632    def setUp(self):
    2733        #print "****set up****"
    … …  
    4450
    4551        #have  a big area covered.
    46 
    4752        mesh_file = tempfile.mktemp(".tsh")
    48        
    4953        points_lat_long = [[-33,152],[-35,152],[-35,150],[-33,150]]
    50        
    5154        spat = Geospatial_data(data_points=points_lat_long,
    5255                               points_are_lats_longs=True)
    5356        points_ab = spat.get_data_points( absolute = True)
    54 
    5557       
    5658        geo =  Geo_reference(56,400000,6000000)
    5759        spat.set_geo_reference(geo)
    58 
    5960        m = Mesh()
    6061        m.add_vertices(spat)
    … …  
    6263        m.generate_mesh(verbose=False)
    6364        m.export_mesh_file(mesh_file)
    64 
    6565       
    6666        #Create shallow water domain
    … …  
    7171        domain.default_order=2
    7272        domain.beta_h = 0
    73 
    74 
    7573        #Set some field values
    7674        #domain.set_quantity('stage', 1.0)
    … …  
    7876        domain.set_quantity('friction', 0.03)
    7977
    80 
    8178        ######################
    8279        # Boundary conditions
    … …  
    8481        domain.set_boundary( {'exterior': B})
    8582
    86 
    8783        ######################
    8884        #Initial condition - with jumps
    89 
    9085        bed = domain.quantities['elevation'].vertex_values
    9186        stage = zeros(bed.shape, Float)
    … …  
    10499        domain.distribute_to_vertices_and_edges()
    105100
    106 
    107101        self.domain = domain
    108102
    … …  
    112106
    113107        self.F = bed
    114 
    115        
     108     
    116109        #sww_file = tempfile.mktemp("")
    117         self.domain.set_name('datatest' + str(time.time()))
     110        self.domain.set_name('tid_P0')
    118111        self.domain.format = 'sww'
    119112        self.domain.smooth = True
    … …  
    127120        self.sww = sww # so it can be deleted
    128121       
     122        #Create another sww file
     123        mesh_file = tempfile.mktemp(".tsh")
     124        points_lat_long = [[-35,152],[-36,152],[-36,150],[-35,150]]
     125        spat = Geospatial_data(data_points=points_lat_long,
     126                               points_are_lats_longs=True)
     127        points_ab = spat.get_data_points( absolute = True)
     128       
     129        geo =  Geo_reference(56,400000,6000000)
     130        spat.set_geo_reference(geo)
     131        m = Mesh()
     132        m.add_vertices(spat)
     133        m.auto_segment()
     134        m.generate_mesh(verbose=False)
     135        m.export_mesh_file(mesh_file)
     136       
     137        #Create shallow water domain
     138        domain = Domain(mesh_file)
     139
     140        os.remove(mesh_file)
     141       
     142        domain.default_order=2
     143        domain.beta_h = 0
     144        #Set some field values
     145        #domain.set_quantity('stage', 1.0)
     146        domain.set_quantity('elevation', -40)
     147        domain.set_quantity('friction', 0.03)
     148
     149        ######################
     150        # Boundary conditions
     151        B = Transmissive_boundary(domain)
     152        domain.set_boundary( {'exterior': B})
     153
     154        ######################
     155        #Initial condition - with jumps
     156        bed = domain.quantities['elevation'].vertex_values
     157        stage = zeros(bed.shape, Float)
     158
     159        h = 30.
     160        for i in range(stage.shape[0]):
     161            if i % 2 == 0:
     162                stage[i,:] = bed[i,:] + h
     163            else:
     164                stage[i,:] = bed[i,:]
     165
     166        domain.set_quantity('stage', stage)
     167        domain.set_quantity('xmomentum', stage*22.0)
     168        domain.set_quantity('ymomentum', stage*55.0)
     169
     170        domain.distribute_to_vertices_and_edges()
     171
     172        self.domain2 = domain
     173
     174        C = domain.get_vertex_coordinates()
     175        self.X2 = C[:,0:6:2].copy()
     176        self.Y2 = C[:,1:6:2].copy()
     177
     178        self.F2 = bed
     179     
     180        #sww_file = tempfile.mktemp("")
     181        domain.set_name('tid_P1')
     182        domain.format = 'sww'
     183        domain.smooth = True
     184        domain.reduction = mean
     185
     186        sww = get_dataobject(domain)
     187        sww.store_connectivity()
     188        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     189        domain.time = 2.
     190        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
     191        self.swwII = sww # so it can be deleted
     192
     193        # print "sww.filename", sww.filename
    129194        #Create a csv file
    130195        self.csv_file = tempfile.mktemp(".csv")
    … …  
    137202        fd.close()
    138203
     204        #Create a csv file
     205        self.csv_fileII = tempfile.mktemp(".csv")
     206        fd = open(self.csv_file,'wb')
     207        writer = csv.writer(fd)
     208        writer.writerow(['LONGITUDE','LATITUDE',STR_VALUE_LABEL,CONT_VALUE_LABEL,'ROOF_TYPE',WALL_TYPE_LABEL, SHORE_DIST_LABEL])
     209        writer.writerow(['151.5','-34','199770','130000','Metal','Timber',20.])
     210        writer.writerow(['151','-34.5','150000','76000','Metal','Double Brick',200.])
     211        writer.writerow(['151','-34.25','150000','76000','Metal','Brick Veneer',200.])
     212        writer.writerow(['151.5','-35.5','199770','130000','Metal','Timber',20.])
     213        fd.close()
    139214       
    140215    def tearDown(self):
    … …  
    147222            # fails there are no sww files in the anuga directory
    148223            os.remove(self.sww.filename)
     224            os.remove(self.swwII.filename)
    149225        except OSError:
    150226            pass
    … …  
    387463        results_dic = edm.calc_damage_and_costs(verbose_csv=True)
    388464        #print "results_dic",results_dic
     465       
     466    def test_calc_max_depth_and_momentum(self):
     467        sww_file = "tid" # self.domain.get_name() + "." + self.domain.format
     468        points_lat_long = [[-34, 151.5],[-35.5, 151.5],[-50, 151]]
     469        spat = Geospatial_data(data_points=points_lat_long,
     470                               points_are_lats_longs=True)
     471        points_ab = spat.get_data_points( absolute = True)
     472        deps, _ = calc_max_depth_and_momentum(sww_file,
     473                                              points_ab,
     474                                              verbose=self.verbose,
     475                                              use_cache = False)
     476
     477        # Test values based on returned results, so not an excellent test
     478       
     479        assert allclose(deps[0],0.113204555211)
     480        assert allclose(deps[1],11.3215)
     481        assert allclose(deps[2],0.0) # this value is outside both sww files
     482       
    389483#-------------------------------------------------------------
    390484if __name__ == "__main__":
    391     #suite = unittest.makeSuite(Test_inundation_damage,'test_in_damage2')
     485    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
     486        Test_inundation_damage.verbose=True
     487        saveout = sys.stdout   
     488        filename = ".temp_verbose"
     489        fid = open(filename, 'w')
     490        sys.stdout = fid
     491    else:
     492        pass
     493    #suite = unittest.makeSuite(Test_inundation_damage,'test_calc_max_depth_and_momentum')
    392494    suite = unittest.makeSuite(Test_inundation_damage,'test')
    393495    runner = unittest.TextTestRunner()
    394496    runner.run(suite)
    395497
     498    # Cleaning up
     499    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
     500        sys.stdout = saveout
     501        fid.close()
     502        os.remove(filename)
Note: See TracChangeset for help on using the changeset viewer.