Changeset 4554

Timestamp:

Jun 21, 2007, 3:04:42 PM (17 years ago)

Author:

ole

Message:

Implemented get_maximum_inundation_elevation and
get_maximum_inundation_location for use with sww files, added tests and
updated user manual.

Location:

anuga_core

Files:

• documentation/user_manual/anuga_user_manual.tex (modified) (3 diffs)
• source/anuga/shallow_water/data_manager.py (modified) (5 diffs)
• source/anuga/shallow_water/shallow_water_domain.py (modified) (2 diffs)
• source/anuga/shallow_water/test_shallow_water_domain.py (modified) (5 diffs)

anuga_core/documentation/user_manual/anuga_user_manual.tex

@@ -2238 +2238 @@
 
 \subsection{Diagnostics}
+\label{sec:diagnostics}
 
 
@@ -2369 +2370 @@
 
 
+\section{Queries of SWW model output files}  
+After a model has been run, it is often useful to extract various information from the sww 
+output file (see Section \ref{sec:sww format}). This is typically more convenient than using the
+diagnostics described in Section \ref{sec:diagnostics} which rely on the model running - something 
+that can be very time consuming. The sww files are easy and quick to read and offer much information 
+about the model results such as runup heights, time histories of selected quantities, 
+flow through cross sections and much more.
+
+\begin{funcdesc}{get\_maximum\_inundation\_elevation}{filename, polygon=None, 
+    time_interval=None, verbose=False}
+  Module: \module{shallow\_water.data\_manager}
+
+  Return highest elevation where depth ($h$) > 0
+
+  Example to find maximum runup elevation:\\    
+  max_runup = get_maximum_inundation_elevation(filename,
+  polygon=None,
+  time_interval=None,
+  verbose=False)
+
+   
+  filename is a NetCDF sww file containing ANUGA model output.                                                       
+  Optional arguments polygon and time_interval restricts the maximum runup calculation
+  to a points that lie within the specified polygon and time interval.
+
+  If no inundation is found within polygon and time_interval the return value
+  is None signifying "No Runup" or "Everything is dry".
+
+  See doc string for general function get_maximum_inundation_data for details.
+\end{funcdesc}
+
+
+\begin{funcdesc}{get\_maximum\_inundation\_location}{filename, polygon=None, 
+    time_interval=None, verbose=False}
+  Module: \module{shallow\_water.data\_manager}
+
+  Return location (x,y) of highest elevation where depth ($h$) > 0
+
+  Example to find maximum runup location:\\    
+  max_runup_location = get_maximum_inundation_location(filename,
+  polygon=None,
+  time_interval=None,
+  verbose=False)
+
+   
+  filename is a NetCDF sww file containing ANUGA model output.                                                       
+  Optional arguments polygon and time_interval restricts the maximum runup calculation
+  to a points that lie within the specified polygon and time interval.
+
+  If no inundation is found within polygon and time_interval the return value
+  is None signifying "No Runup" or "Everything is dry".
+
+  See doc string for general function get_maximum_inundation_data for details.
+\end{funcdesc}
+
+
+\begin{funcdesc}{sww2time\_series}{}
+To appear
+\end{funcdesc}
+ 
+  
 
 \section{Other}
@@ -2491 +2553 @@
 
 \subsection{SWW and TMS Formats}
+\label{sec:sww format}
 
 The SWW and TMS formats are both NetCDF formats, and are of key

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -69 +69 @@
 from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \
      searchsorted, zeros, allclose, around, reshape, transpose, sort, \
-     NewAxis, ArrayType, compress, take, arange, argmax
+     NewAxis, ArrayType, compress, take, arange, argmax, alltrue
 from Scientific.IO.NetCDF import NetCDFFile
 #from shutil import copy
@@ -5382 +5382 @@
 
 
-def get_maximum_inundation_elevation(sww_filename, region=None, timesteps=None, verbose=False):
+def get_maximum_inundation_elevation(filename,
+                                     polygon=None,
+                                     time_interval=None,
+                                     verbose=False):
+    
+    """Return highest elevation where depth > 0
+    
+    Usage:
+    max_runup = get_maximum_inundation_elevation(filename,
+                                                 polygon=None,
+                                                 time_interval=None,
+                                                 verbose=False)
+
+    filename is a NetCDF sww file containing ANUGA model output.                                                       
+    Optional arguments polygon and time_interval restricts the maximum runup calculation
+    to a points that lie within the specified polygon and time interval.
+
+    If no inundation is found within polygon and time_interval the return value
+    is None signifying "No Runup" or "Everything is dry".
+
+    See general function get_maximum_inundation_data for details.
+    
+    """
+    
+    runup, _ = get_maximum_inundation_data(filename,
+                                           polygon=polygon,
+                                           time_interval=time_interval,
+                                           verbose=verbose)
+    return runup
+
+
+
+
+def get_maximum_inundation_location(filename,
+                                    polygon=None,
+                                    time_interval=None,
+                                    verbose=False):
+    """Return location of highest elevation where h > 0
+    
+    
+    Usage:
+    max_runup_location = get_maximum_inundation_location(filename,
+                                                         polygon=None,
+                                                         time_interval=None,
+                                                         verbose=False)
+
+    filename is a NetCDF sww file containing ANUGA model output.
+    Optional arguments polygon and time_interval restricts the maximum runup calculation
+    to a points that lie within the specified polygon and time interval.
+
+    If no inundation is found within polygon and time_interval the return value
+    is None signifying "No Runup" or "Everything is dry".
+
+    See general function get_maximum_inundation_data for details.
+    """
+    
+    _, max_loc = get_maximum_inundation_data(filename,
+                                             polygon=polygon,
+                                             time_interval=time_interval,
+                                             verbose=verbose)
+    return max_loc
+    
+
+
+def get_maximum_inundation_data(filename, polygon=None, time_interval=None,
+                                verbose=False):
     """Compute maximum run up height from sww file.
+
+
+    Usage:
+    runup, location = get_maximum_inundation_data(filename,
+                                                  polygon=None,
+                                                  time_interval=None,
+                                                  verbose=False)
+    
 
     Algorithm is as in get_maximum_inundation_elevation from shallow_water_domain
@@ -5389 +5462 @@
     runup height over multiple timesteps. 
     
-    Optional argument region restricts this to specified polygon.
-    Optional argument timesteps restricts this to a specified time step (an index)
-    or time interval (a list of indices).
+    Optional arguments polygon and time_interval restricts the maximum runup calculation
+    to a points that lie within the specified polygon and time interval.
+
+    If no inundation is found within polygon and time_interval the return value
+    is None signifying "No Runup" or "Everything is dry".
     """
 
@@ -5398 +5473 @@
     # Water depth below which it is considered to be 0 in the model
     # FIXME (Ole): Allow this to be specified as a keyword argument as well
-    
+
+    from anuga.utilities.polygon import inside_polygon    
     from anuga.config import minimum_allowed_height
 
-    if region is not None:
-        msg = 'region not yet implemented in get_maximum_inundation_elevation'
-        raise Exception, msg
-
-
-    root, extension = os.path.splitext(sww_filename)
+
+    root, extension = os.path.splitext(filename)
     if extension == '':
-        sww_filename += '.sww'
+        filename += '.sww'
     
     # Read sww file
     if verbose: 
-        print 'Reading from %s' %sww_filename
+        print 'Reading from %s' %filename
     
     from Scientific.IO.NetCDF import NetCDFFile
-    fid = NetCDFFile(sww_filename)
+    fid = NetCDFFile(filename)
 
     # Get extent and reference
+    volumes = fid.variables['volumes'][:]    
     x = fid.variables['x'][:]
     y = fid.variables['y'][:]
-    volumes = fid.variables['volumes'][:]
-
-
-    time = fid.variables['time'][:]
-    if timesteps is not None:
-        if type(timesteps) is list or type(timesteps) is ArrayType:  
-            timesteps = ensure_numeric(timesteps, Int)
-        elif type(timesteps) is type(0):
-            timesteps = [timesteps]
-        else:
-            msg = 'timesteps must be either an integer or a sequence of integers. '
-            msg += 'I got timesteps==%s, %s' %(timesteps, type(timesteps))
-            raise Exception, msg
-    else:
-        # Take them all
-        timesteps = arange(len(time))    
-    
 
     # Get the relevant quantities
@@ -5443 +5499 @@
 
 
+    # Spatial restriction
+    if polygon is not None:
+        msg = 'polygon must be a sequence of points.'
+        assert len(polygon[0]) == 2, msg
+        # FIXME (Ole): Make a generic polygon input check in polygon.py and call it here
+        
+        points = concatenate((x[:,NewAxis], y[:,NewAxis]), axis=1)
+
+        point_indices = inside_polygon(points, polygon)
+
+        # Restrict quantities to polygon
+        elevation = take(elevation, point_indices)
+        stage = take(stage, point_indices, axis=1)
+
+        # Get info for location of maximal runup
+        points_in_polygon = take(points, point_indices)
+        x = points_in_polygon[:,0]
+        y = points_in_polygon[:,1]        
+    else:
+        # Take all points
+        point_indices = arange(len(x))
+        
+
+    # Temporal restriction
+    time = fid.variables['time'][:]
+    all_timeindices = arange(len(time))        
+    if time_interval is not None:
+        
+        msg = 'time_interval must be a sequence of length 2.'
+        assert len(time_interval) == 2, msg
+        msg = 'time_interval %s must not be decreasing.' %(time_interval)        
+        assert time_interval[1] >= time_interval[0], msg
+        
+        msg = 'Specified time interval [%.8f:%.8f]' %tuple(time_interval)
+        msg += ' must does not match model time interval: [%.8f, %.8f]\n'\
+               %(time[0], time[-1])
+        if time_interval[1] < time[0]: raise ValueError(msg)
+        if time_interval[0] > time[-1]: raise ValueError(msg)
+
+        # Take time indices corresponding to interval (& is bitwise AND)
+        timesteps = compress((time_interval[0] <= time) & (time <= time_interval[1]),
+                             all_timeindices)
+
+
+        msg = 'time_interval %s did not include any model timesteps.' %(time_interval)        
+        assert not alltrue(timesteps == 0), msg
+
+
+    else:
+        # Take them all
+        timesteps = all_timeindices
+    
+
+    fid.close()
+
     # Compute maximal runup for each timestep
     maximal_runup = None
+    maximal_runup_location = None
     for i in timesteps:
-        #print 'time', time[i]
-        depth = stage[i,:] - elevation[:]
+        depth = stage[i,:] - elevation 
     
         # Get wet nodes i.e. nodes with depth>0 within given region and timesteps
         wet_nodes = compress(depth > minimum_allowed_height, arange(len(depth)))
 
-        # Find maximum elevation among wet nodes and return
-        #runup_index = argmax(take(elevation, wet_nodes)) #FIXME Maybe get loc as well
-        runup = max(take(elevation, wet_nodes))
+        if alltrue(wet_nodes == 0):
+            runup = None
+        else:    
+            # Find maximum elevation among wet nodes
+            wet_elevation = take(elevation, wet_nodes)
+
+            runup_index = argmax(wet_elevation)
+            runup = max(wet_elevation)
+            assert wet_elevation[runup_index] == runup # Must always be True
 
         if runup > maximal_runup:
-            maximal_runup = runup  # This works even if maximal_runup is None
-                
-
-    return maximal_runup
+            maximal_runup = runup      # This works even if maximal_runup is None
+
+            # Record location
+            wet_x = take(x, wet_nodes)
+            wet_y = take(y, wet_nodes)            
+            maximal_runup_location = [wet_x[runup_index], wet_y[runup_index]]                
+
+    return maximal_runup, maximal_runup_location
 
 

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -342 +342 @@
 
     def get_maximum_inundation_location(self, indices=None):
-        """Return highest elevation where h > 0
+        """Return location of highest elevation where h > 0
 
         Optional argument:
@@ -348 +348 @@
 
         Usage:
-            q = get_maximum_inundation_elevation()
+            q = get_maximum_inundation_location()
 
         Note, centroid values are used for this operation            

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -7 +7 @@
 from Numeric import allclose, alltrue, array, zeros, ones, Float, take
 from anuga.utilities.numerical_tools import mean
+from anuga.utilities.polygon import is_inside_polygon
 
 from shallow_water_domain import *
@@ -951 +952 @@
         """test_get_maximum_inundation_from_sww(self)
 
-        Test of get_maximum_inundation_elevation(sww_filename) from data_manager.py
+        Test of get_maximum_inundation_elevation()
+        and get_maximum_inundation_location() from data_manager.py
         
         This is based on test_get_maximum_inundation_3(self) but works with the
         stored results instead of with the internal data structure.
-        
+
+        This test uses the underlying get_maximum_inundation_data for tests
         """
 
         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
-        from data_manager import get_maximum_inundation_elevation
+        from data_manager import get_maximum_inundation_elevation, get_maximum_inundation_location, get_maximum_inundation_data
         
 
@@ -1033 +1036 @@
         assert allclose(q, initial_runup_height, rtol = 1.0/N), msg 
 
-        
-
+
+        # Test error condition if time interval is out
+        try:
+            q = get_maximum_inundation_elevation('runup_test.sww',
+                                                 time_interval=[2.0, 3.0])
+        except ValueError:
+            pass
+        else:
+            msg = 'should have caught wrong time interval'
+            raise Exception, msg
+
+        # Check correct time interval
+        q, loc = get_maximum_inundation_data('runup_test.sww',
+                                             time_interval=[0.0, 3.0])        
+        msg = 'We got %f, should have been %f' %(q, initial_runup_height)
+        assert allclose(q, initial_runup_height, rtol = 1.0/N), msg
+        assert allclose(-loc[0]/2, q) # From topography formula         
+        
 
         #--------------------------------------------------------------
@@ -1046 +1065 @@
         #--------------------------------------------------------------
         q_max = None
-        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
+        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
+                               skip_initial_step = True): 
             q = domain.get_maximum_inundation_elevation()
             if q > q_max: q_max = q
@@ -1064 +1084 @@
         assert allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy
 
-        q = get_maximum_inundation_elevation('runup_test.sww', timesteps=31)
+        q, loc = get_maximum_inundation_data('runup_test.sww', time_interval=[3.0, 3.0])
         msg = 'We got %f, should have been %f' %(q, final_runup_height)
-        assert allclose(q, final_runup_height, rtol=1.0/N), msg        
-
+        assert allclose(q, final_runup_height, rtol=1.0/N), msg
+        #print 'loc', loc, q        
+        assert allclose(-loc[0]/2, q) # From topography formula         
 
         q = get_maximum_inundation_elevation('runup_test.sww')
+        loc = get_maximum_inundation_location('runup_test.sww')        
         msg = 'We got %f, should have been %f' %(q, q_max)
         assert allclose(q, q_max, rtol=1.0/N), msg
-
-        q = get_maximum_inundation_elevation('runup_test.sww', timesteps=range(32))
+        #print 'loc', loc, q
+        assert allclose(-loc[0]/2, q) # From topography formula 
+
+        
+
+        q = get_maximum_inundation_elevation('runup_test.sww', time_interval=[0, 3])
         msg = 'We got %f, should have been %f' %(q, q_max)
-        assert allclose(q, q_max, rtol=1.0/N), msg                
-
-
-
+        assert allclose(q, q_max, rtol=1.0/N), msg
+
+
+        # Check polygon mode
+        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]] # Runup region
+        q = get_maximum_inundation_elevation('runup_test.sww',
+                                             polygon = polygon,
+                                             time_interval=[0, 3])
+        msg = 'We got %f, should have been %f' %(q, q_max)
+        assert allclose(q, q_max, rtol=1.0/N), msg
+
+        
+        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]] # Offshore region
+        q, loc = get_maximum_inundation_data('runup_test.sww',
+                                             polygon = polygon,
+                                             time_interval=[0, 3])
+        msg = 'We got %f, should have been %f' %(q, -0.475)
+        assert allclose(q, -0.475, rtol=1.0/N), msg
+        assert is_inside_polygon(loc, polygon)
+        assert allclose(-loc[0]/2, q) # From topography formula         
+
+
+        polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]] # Dry region
+        q, loc = get_maximum_inundation_data('runup_test.sww',
+                                             polygon = polygon,
+                                             time_interval=[0, 3])
+        msg = 'We got %s, should have been None' %(q)
+        assert q is None, msg
+        msg = 'We got %s, should have been None' %(loc)
+        assert loc is None, msg        
+
+        # Check what happens if no time point is within interval
+        try:
+            q = get_maximum_inundation_elevation('runup_test.sww', time_interval=[2.8, 2.8])
+        except AssertionError:
+            pass
+        else:
+            msg = 'Time interval should have raised an exception'
+            raise msg
+            
 
     def test_another_runup_example(self):