Changeset 1919

Timestamp:

Oct 14, 2005, 11:26:21 AM (19 years ago)

Author:

ole

Message:

Implemented arbitrary expressions for sww2dem using code from changeset:1916

Location:

inundation/pyvolution

Files:

• data_manager.py (modified) (3 diffs)
• domain.py (modified) (1 diff)
• test_data_manager.py (modified) (1 diff)
• test_util.py (modified) (1 diff)
• util.py (modified) (1 diff)

inundation/pyvolution/data_manager.py

@@ -1248 +1248 @@
 
 
-    if quantity is given, out values from quantity otherwise default to
-    elevation
+    The parameter quantity must be the name of an existing quantity or
+    an expression involving existing quantities. The default is
+    'elevation'. 
 
     if timestep (an index) is given, output quantity at that timestep
@@ -1262 +1263 @@
     from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue 
     from utilities.polygon import inside_polygon
+    from util import apply_expression_to_dictionary
     
     msg = 'Format must be either asc or ers'
@@ -1315 +1317 @@
 
     #Get quantity and reduce if applicable
-    if verbose: print 'Reading quantity %s' %quantity
-
-    if quantity.lower() == 'depth':
-        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
-    else:
-        q = fid.variables[quantity][:]
+    if verbose: print 'Processing quantity %s' %quantity
+
+    #Turn NetCDF objects into Numeric arrays
+    quantity_dict = {}
+    for name in fid.variables.keys():
+
+        if name not in ['x', 'y', 'volumes', 'time']:
+            quantity_dict[name] = fid.variables[name][:]
+
+    #print quantity_dict     
+   
+
+    q = apply_expression_to_dictionary(quantity, quantity_dict)
 
 

inundation/pyvolution/domain.py

@@ -745 +745 @@
     """
 
-    #Replace quantity names with actual quantities
-    quantities = domain.quantities
-    for quantity_name in quantities.keys():
-        expression = expression.replace(quantity_name,
-                                        'quantities["%s"]' %quantity_name)
-
-    #Evaluate and return
-    return eval(expression)   
+    from util import apply_expression_to_dictionary
+    return apply_expression_to_dictionary(expression, domain.quantities)
     
 

inundation/pyvolution/test_data_manager.py

@@ -1351 +1351 @@
 
 
+    def test_sww2dem_asc_derived_quantity(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+
+        This tests the use of derived quantities
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+        self.domain.filename = 'datatest'
+
+        prjfile = self.domain.filename + '_depth.prj'
+        ascfile = self.domain.filename + '_depth.asc'
+        swwfile = self.domain.filename + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+        self.domain.set_quantity('stage', 0.0)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+
+        sww = get_dataobject(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep('stage')
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, 'r')
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+
+        #Export to ascii/prj files
+        sww2dem(self.domain.filename,
+                basename_out = 'datatest_depth',
+                quantity = 'stage - elevation',
+                cellsize = cellsize,
+                reduction = min,
+                format = 'asc',
+                verbose = False)
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'
+
+
+        #Check grid values (where applicable)
+        for j in range(5):
+            if j%2 == 0:
+                L = lines[6+j].strip().split()
+                jj = 4-j
+                for i in range(5):
+                    if i%2 == 0:
+                        index = jj/2 + i/2*3
+                        val0 = stage[0,index] - z[index]
+                        val1 = stage[1,index] - z[index]
+
+                        #print i, j, index, ':', L[i], val0, val1
+                        assert allclose(float(L[i]), min(val0, val1))
+
+
+        fid.close()
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        #os.remove(swwfile)
+
+
+
+
 
     def test_sww2dem_asc_missing_points(self):

inundation/pyvolution/test_util.py

@@ -972 +972 @@
 
 
-
+    def test_apply_expression_to_dictionary(self):
+
+        foo = array([[1,2,3],
+                     [4,5,6]])
+
+        bar = array([[-1,0,5],
+                     [6,1,1]])                  
+
+        D = {'X': foo, 'Y': bar}
+
+        Z = apply_expression_to_dictionary('X+Y', D)        
+        assert allclose(Z, foo+bar)
+
+        Z = apply_expression_to_dictionary('X*Y', D)        
+        assert allclose(Z, foo*bar)        
+
+        Z = apply_expression_to_dictionary('4*X+Y', D)        
+        assert allclose(Z, 4*foo+bar)        
+
+        #Check exceptions
+        try:
+            #Wrong name
+            Z = apply_expression_to_dictionary('4*X+A', D)        
+        except NameError:
+            pass
+        else:
+            msg = 'Should have raised a NameError Exception'
+            raise msg
+
+
+        try:
+            #Wrong order
+            Z = apply_expression_to_dictionary(D, '4*X+A')        
+        except AssertionError:
+            pass
+        else:
+            msg = 'Should have raised a AssertionError Exception'
+            raise msg        
+        
 
 

inundation/pyvolution/util.py

@@ -457 +457 @@
     sww.store_timestep('stage')
 
+
+
+
+
+
+
+def apply_expression_to_dictionary(expression, _):
+    """Apply arbitrary expression to values of dictionary _
+
+    expression: Arbitrary, e.g. arithmetric, expression relating keys
+                from dictionary _. A key must not be the single character '_'
+
+    _: Dictionary of variable names used in expression and their
+       corresponding values. The dictinary name is chosen as the
+       inconspicuous symbol _ in order to avoid single letter variable
+       to inabvertently substitute parts of the dictionary name
+    
+
+    Given an expression in terms of the keys, replace key by the
+    corresponding values and evaluate.
+
+    Values in dictionary must support operators given in expression
+    e.g. by overloading
+    """
+
+    #FIXME: Varibles may substitute into other variables e.g. x into xmomentum.
+    
+    import types
+
+    assert isinstance(expression, basestring)
+    assert type(_) == types.DictType
+
+    #Replace key names with values
+    for key in _.keys():
+        msg = 'Key must not be the single character "_"'
+        assert key != '_', msg
+        expression = expression.replace(key,
+                                        '_["%s"]' %key)
+
+
+    #Evaluate and return
+    try:
+        return eval(expression)
+    except NameError, e:
+        msg = 'Expression "%s" could not be evaluated: %s' %(expression, e)
+        raise NameError, msg
+    
+
+
+
+
+
+
+
 ####################################################################
 #Python versions of function that are also implemented in util_gateway.c