Changeset 2891 for inundation/pyvolution

Timestamp:

May 17, 2006, 12:11:37 PM (18 years ago)

Author:

ole

Message:

Added sww2pts and a test. This can be used to export values at centroids.

Location:

inundation/pyvolution

Files:

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

inundation/pyvolution/data_manager.py

@@ -1665 +1665 @@
 
 
-    #Get quantity and reduce if applicable
+    # Get quantity and reduce if applicable
     if verbose: print 'Processing quantity %s' %quantity
 
-    #Turn NetCDF objects into Numeric arrays
+    # Turn NetCDF objects into Numeric arrays
     quantity_dict = {}
     for name in fid.variables.keys():
@@ -1674 +1674 @@
 
 
+
+    # Convert quantity expression to quantities found in sww file    
     q = apply_expression_to_dictionary(quantity, quantity_dict)
 
@@ -1909 +1911 @@
             verbose = verbose,
             origin = origin,
-        datum = datum,
-        format = 'ers')
+            datum = datum,
+            format = 'ers')
 ################################# END COMPATIBILITY ##############
+
+
+
+def sww2pts(basename_in, basename_out=None,
+            data_points=None,
+            quantity=None,
+            timestep=None,
+            reduction=None,
+            NODATA_value=-9999,
+            verbose=False,
+            origin=None):
+            #datum = 'WGS84')
+
+
+    """Read SWW file and convert to interpolated values at selected points
+
+    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
+
+    if reduction is given use that to reduce quantity over all timesteps.
+
+    data_points (Nx2 array) give locations of points where quantity is to be computed.
+    
+    """
+
+    import sys
+    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
+    from Numeric import array2string
+
+    from utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
+    from util import apply_expression_to_dictionary
+
+    from geospatial_data import Geospatial_data
+
+    if quantity is None:
+        quantity = 'elevation'
+
+    if reduction is None:
+        reduction = max
+
+    if basename_out is None:
+        basename_out = basename_in + '_%s' %quantity
+
+    swwfile = basename_in + '.sww'
+    ptsfile = basename_out + '.pts'
+
+    # Read sww file
+    if verbose: print 'Reading from %s' %swwfile
+    from Scientific.IO.NetCDF import NetCDFFile
+    fid = NetCDFFile(swwfile)
+
+    # Get extent and reference
+    x = fid.variables['x'][:]
+    y = fid.variables['y'][:]
+    volumes = fid.variables['volumes'][:]
+
+    number_of_timesteps = fid.dimensions['number_of_timesteps']
+    number_of_points = fid.dimensions['number_of_points']
+    if origin is None:
+
+        # Get geo_reference
+        # sww files don't have to have a geo_ref
+        try:
+            geo_reference = Geo_reference(NetCDFObject=fid)
+        except AttributeError, e:
+            geo_reference = Geo_reference() #Default georef object
+
+        xllcorner = geo_reference.get_xllcorner()
+        yllcorner = geo_reference.get_yllcorner()
+        zone = geo_reference.get_zone()
+    else:
+        zone = origin[0]
+        xllcorner = origin[1]
+        yllcorner = origin[2]
+
+
+
+    # FIXME: Refactor using code from file_function.statistics
+    # Something like print swwstats(swwname)
+    if verbose:
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        times = fid.variables['time'][:]
+        print '------------------------------------------------'
+        print 'Statistics of SWW file:'
+        print '  Name: %s' %swwfile
+        print '  Reference:'
+        print '    Lower left corner: [%f, %f]'\
+              %(xllcorner, yllcorner)
+        print '    Start time: %f' %fid.starttime[0]
+        print '  Extent:'
+        print '    x [m] in [%f, %f], len(x) == %d'\
+              %(min(x.flat), max(x.flat), len(x.flat))
+        print '    y [m] in [%f, %f], len(y) == %d'\
+              %(min(y.flat), max(y.flat), len(y.flat))
+        print '    t [s] in [%f, %f], len(t) == %d'\
+              %(min(times), max(times), len(times))
+        print '  Quantities [SI units]:'
+        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
+            q = fid.variables[name][:].flat
+            print '    %s in [%f, %f]' %(name, min(q), max(q))
+
+
+
+    # Get quantity and reduce if applicable
+    if verbose: print 'Processing quantity %s' %quantity
+
+    # Turn NetCDF objects into Numeric arrays
+    quantity_dict = {}
+    for name in fid.variables.keys():
+        quantity_dict[name] = fid.variables[name][:]
+
+
+
+    # Convert quantity expression to quantities found in sww file    
+    q = apply_expression_to_dictionary(quantity, quantity_dict)
+
+
+
+    if len(q.shape) == 2:
+        # q has a time component and needs to be reduced along
+        # the temporal dimension
+        if verbose: print 'Reducing quantity %s' %quantity
+        q_reduced = zeros( number_of_points, Float )
+
+        for k in range(number_of_points):
+            q_reduced[k] = reduction( q[:,k] )
+
+        q = q_reduced
+
+    # Post condition: Now q has dimension: number_of_points
+    assert len(q.shape) == 1
+    assert q.shape[0] == number_of_points
+
+
+    if verbose:
+        print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
+
+
+    # Create grid and update xll/yll corner and x,y
+    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
+    assert len(vertex_points.shape) == 2
+
+    # Interpolate
+    from fit_interpolate.interpolate import Interpolate
+    interp = Interpolate(vertex_points, volumes, verbose = verbose)
+
+    # Interpolate using quantity values
+    if verbose: print 'Interpolating'
+    interpolated_values = interp.interpolate(q, data_points).flat
+
+
+    if verbose:
+        print 'Interpolated values are in [%f, %f]' %(min(interpolated_values),
+                                                      max(interpolated_values))
+
+    # Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
+    P = interp.mesh.get_boundary_polygon()
+    outside_indices = outside_polygon(data_points, P, closed=True)
+
+    for i in outside_indices:
+        interpolated_values[i] = NODATA_value
+
+
+    # Store results    
+    G = Geospatial_data(data_points=data_points,
+                        attributes=interpolated_values)
+
+    G.export_points_file(ptsfile, absolute = True)
 
 

inundation/pyvolution/test_data_manager.py

@@ -1996 +1996 @@
 
         #Export to ers files
-        #sww2ers(self.domain.filename,
-        #        quantity = 'elevation',
-        #        cellsize = cellsize,
-        #        verbose = False)
-
         sww2dem(self.domain.filename,
                 quantity = 'elevation',
@@ -2049 +2044 @@
         os.remove(self.domain.filename + '_elevation')
         os.remove(self.domain.filename + '_elevation.ers')
+
+
+
+    def test_sww2pts_centroids(self):
+        """Test that sww information can be converted correctly to pts data at specified coordinates
+        - in this case, the centroids.
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate, NewAxis
+        from Scientific.IO.NetCDF import NetCDFFile
+        from geospatial_data import Geospatial_data
+
+        # Used for points that lie outside mesh
+        NODATA_value = 1758323
+
+        # Setup
+        self.domain.filename = 'datatest'
+
+        ptsfile = self.domain.filename + '_elevation.pts'
+        swwfile = self.domain.filename + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.smooth = True #self.set_store_vertices_uniquely(False)
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+
+        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')
+
+        # Check contents in NetCDF
+        fid = NetCDFFile(sww.filename, 'r')
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        elevation = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+        volumes = fid.variables['volumes'][:]
+
+
+        # Invoke interpolation for vertex points       
+        points = concatenate( (x[:,NewAxis],y[:,NewAxis]), axis=1 )
+        sww2pts(self.domain.filename,
+                quantity = 'elevation',
+                data_points = points,
+                NODATA_value = NODATA_value,
+                verbose = False)
+        ref_point_values = elevation
+        point_values = Geospatial_data(ptsfile).get_attributes()
+        #print 'P', point_values
+        #print 'Ref', ref_point_values        
+        assert allclose(point_values, ref_point_values)        
+
+
+
+        # Invoke interpolation for centroids
+        points = self.domain.get_centroid_coordinates()
+        #print points
+        sww2pts(self.domain.filename,
+                quantity = 'elevation',
+                data_points = points,
+                NODATA_value = NODATA_value,
+                verbose = False)
+        ref_point_values = [-0.5, -0.5, -1, -1, -1, -1, -1.5, -1.5]   #At centroids
+
+        
+        point_values = Geospatial_data(ptsfile).get_attributes()
+        #print 'P', point_values
+        #print 'Ref', ref_point_values        
+        assert allclose(point_values, ref_point_values)        
+
+
+
+        fid.close()
+
+        #Cleanup
+        os.remove(sww.filename)
+        os.remove(ptsfile)
+
 
 

inundation/pyvolution/util.py

@@ -570 +570 @@
       Name = gauges_timeseries followed by gauge name 
     - latex file will be generated in same directory as where script is
-      run (usually production scenario direcotry.
+      run (usually production scenario directory.
       Name = latexoutputlabel_id.tex