Changeset 8970

Timestamp:

Sep 11, 2013, 7:43:43 AM (11 years ago)

Author:

steve

Message:

Lots of changes, most important lowering the minimum allowed height

Location:

trunk/anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/quantity.py (modified) (4 diffs)
• abstract_2d_finite_volumes/test_quantity.py (modified) (3 diffs)
• anuga_exceptions.py (modified) (1 diff)
• config.py (modified) (1 diff)
• fit_interpolate/fit.py (modified) (1 diff)
• fit_interpolate/interpolate2d.py (modified) (1 diff)
• shallow_water/test_data_manager.py (modified) (3 diffs)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -16 +16 @@
 
 import types
+import os.path
 
 from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
@@ -349 +350 @@
         import numpy as np
 
-
+        plt.clf()
         plt.tripcolor(X, Y, V, A)
         plt.colorbar()
@@ -1042 +1043 @@
         msg = 'Filename must be a text string'
         assert isinstance(filename, basestring), msg
+        
+        msg = 'Extension should be .pts .dem, .csv, or txt'
+        assert os.path.splitext(filename)[1] in ['.pts', '.dem', '.csv', '.txt'], msg
+        
 
         if location != 'vertices':
@@ -1083 +1088 @@
                                    indices, use_cache=use_cache,
                                    verbose=verbose)
-
+        
+        
+        
+    def set_values_from_asc_file(self,
+                             filename,
+                             location='vertices',
+                             indices=None,
+                             verbose=False):
+        
+        """Read Digital Elevation model from the following ASCII format (.asc or .grd)
+    
+        Example:
+        ncols         3121
+        nrows         1800
+        xllcorner     722000
+        yllcorner     5893000
+        cellsize      25
+        NODATA_value  -9999
+        138.3698 137.4194 136.5062 135.5558 ..........
+    
+    
+        An accompanying file with same basename but extension .prj must exist
+        and is used to fix the UTM zone, datum, false northings and eastings.
+    
+        The prj format is assumed to be as
+    
+        Projection    UTM
+        Zone          56
+        Datum         WGS84
+        Zunits        NO
+        Units         METERS
+        Spheroid      WGS84
+        Xshift        0.0000000000
+        Yshift        10000000.0000000000
+        Parameters
+        """
+
+
+
+
+        msg = 'Filename must be a text string'
+        assert isinstance(filename, basestring), msg
+        
+        msg = 'Extension should be .grd or asc'
+        assert os.path.splitext(filename)[1] in ['.grd', '.asc'], msg
+        
+
+        msg = "set_values_from_grd_file is only defined for location='vertices' or 'centroids'"
+        assert location in ['vertices', 'centroids'], msg
+
+            
+        if location == 'centroids':
+            points = self.domain.centroid_coordinates
+        else:
+            points = self.domain.vertex_coordinates
+            
+    
+        root = filename[:-4]
+    
+#         # Read Meta data
+#         if verbose: log.critical('Reading METADATA from %s' % (root + '.prj'))
+#     
+#         metadatafile = open(root + '.prj')
+#         metalines = metadatafile.readlines()
+#         metadatafile.close()
+#     
+#         L = metalines[0].strip().split()
+#         assert L[0].strip().lower() == 'projection'
+#         projection = L[1].strip()                   #TEXT
+#     
+#         L = metalines[1].strip().split()
+#         assert L[0].strip().lower() == 'zone'
+#         zone = int(L[1].strip())
+#     
+#         L = metalines[2].strip().split()
+#         assert L[0].strip().lower() == 'datum'
+#         datum = L[1].strip()                        #TEXT
+#     
+#         L = metalines[3].strip().split()
+#         assert L[0].strip().lower() == 'zunits'     #IGNORE
+#         zunits = L[1].strip()                       #TEXT
+#     
+#         L = metalines[4].strip().split()
+#         assert L[0].strip().lower() == 'units'
+#         units = L[1].strip()                        #TEXT
+#     
+#         L = metalines[5].strip().split()
+#         assert L[0].strip().lower() == 'spheroid'   #IGNORE
+#         spheroid = L[1].strip()                     #TEXT
+#     
+#         L = metalines[6].strip().split()
+#         assert L[0].strip().lower() == 'xshift'
+#         false_easting = float(L[1].strip())
+#     
+#         L = metalines[7].strip().split()
+#         assert L[0].strip().lower() == 'yshift'
+#         false_northing = float(L[1].strip())
+    
+    
+        #Read DEM data
+        datafile = open(filename)
+    
+        if verbose: log.critical('Reading data from %s' % (filename))
+    
+        lines = datafile.readlines()
+        datafile.close()
+    
+        if verbose: log.critical('Got %d lines' % len(lines))
+    
+
+        ncols = int(lines[0].split()[1].strip())
+        nrows = int(lines[1].split()[1].strip())
+
+    
+        # Do cellsize (line 4) before line 2 and 3
+        cellsize = float(lines[4].split()[1].strip())
+    
+        # Checks suggested by Joaquim Luis
+        # Our internal representation of xllcorner
+        # and yllcorner is non-standard.
+        xref = lines[2].split()
+        if xref[0].strip() == 'xllcorner':
+            xllcorner = float(xref[1].strip()) # + 0.5*cellsize # Correct offset
+        elif xref[0].strip() == 'xllcenter':
+            xllcorner = float(xref[1].strip())
+        else:
+            msg = 'Unknown keyword: %s' % xref[0].strip()
+            raise Exception, msg
+    
+        yref = lines[3].split()
+        if yref[0].strip() == 'yllcorner':
+            yllcorner = float(yref[1].strip()) # + 0.5*cellsize # Correct offset
+        elif yref[0].strip() == 'yllcenter':
+            yllcorner = float(yref[1].strip())
+        else:
+            msg = 'Unknown keyword: %s' % yref[0].strip()
+            raise Exception, msg
+    
+        NODATA_value = int(float(lines[5].split()[1].strip()))
+    
+        assert len(lines) == nrows + 6
+    
+  
+        #Store data
+        import numpy
+    
+        datafile = open(filename)
+        Z = numpy.loadtxt(datafile, skiprows=6)
+        datafile.close()
+        
+        #print Z.shape, Z
+    
+        x = num.linspace(xllcorner, xllcorner+cellsize*(nrows-1), nrows)
+        y = num.linspace(yllcorner, yllcorner+cellsize*(ncols-1), ncols)
+        
+        #print x.shape, x
+        #print y.shape, y
+        
+        from  anuga.fit_interpolate.interpolate2d import interpolate2d 
+        
+        #print points 
+        values = interpolate2d(x, y, Z, points, mode='linear', bounds_error=False)
+        
+        #print values
+
+        # Call underlying method using array values
+        if verbose:
+            log.critical('Applying fitted data to quantity')
+            
+        
+        if location == 'centroids':
+            if indices is None:
+                msg = 'Number of values must match number of elements'
+                #assert values.shape[0] == N, msg
+
+                self.centroid_values[:] = values
+            else:
+                msg = 'Number of values must match number of indices'
+                assert values.shape[0] == indices.shape[0], msg
+
+                # Brute force
+                self.centroid_values[indices] = values
+        else:
+            if indices is None:
+                msg = 'Number of values must match number of elements'
+                #assert values.shape[0] == N, msg
+
+                self.vertex_values[:] = values
+            else:
+                msg = 'Number of values must match number of indices'
+                assert values.shape[0] == indices.shape[0], msg
+
+                # Brute force
+                self.vertex_values[indices] = values
+                
     def get_extremum_index(self, mode=None, indices=None):
         """Return index for maximum or minimum value of quantity (on centroids)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -22 +22 @@
 def linear_function(point):
     point = num.array(point)
-    return point[:,0]+point[:,1]
+    return point[:,0]+3*point[:,1]
+    #return point[:,1]
+
+
+def axes2points(x, y):
+    """Generate all combinations of grid point coordinates from x and y axes
+
+    Args:
+        * x: x coordinates (array)
+        * y: y coordinates (array)
+
+    Returns:
+        * P: Nx2 array consisting of coordinates for all
+             grid points defined by x and y axes. The x coordinate
+             will vary the fastest to match the way 2D numpy
+             arrays are laid out by default ('C' order). That way,
+             the x and y coordinates will match a corresponding
+             2D array A when flattened (A.flat[:] or A.reshape(-1))
+
+    Note:
+        Example
+
+        x = [1, 2, 3]
+        y = [10, 20]
+
+        P = [[1, 10],
+             [2, 10],
+             [3, 10],
+             [1, 20],
+             [2, 20],
+             [3, 20]]
+    """
+    import numpy 
+    
+    # Reverse y coordinates to have them start at bottom of array
+    y = numpy.flipud(y)
+
+    # Repeat x coordinates for each y (fastest varying)
+    X = numpy.kron(numpy.ones(len(y)), x)
+
+    # Repeat y coordinates for each x (slowest varying)
+    Y = numpy.kron(y, numpy.ones(len(x)))
+
+    # Check
+    N = len(X)
+    assert len(Y) == N
+
+    # Create Nx2 array of x and y coordinates
+    X = numpy.reshape(X, (N, 1))
+    Y = numpy.reshape(Y, (N, 1))
+    P = numpy.concatenate((X, Y), axis=1)
+
+    # Return
+    return P
+
 
 
@@ -1333 +1387 @@
 
 
+    def test_set_values_from_asc_vertices(self):
+
+        quantity= Quantity(self.mesh_onslow)
+
+
+        """ Format of asc file 
+        ncols         11
+        nrows         12
+        xllcorner     240000
+        yllcorner     7620000
+        cellsize      6000
+        NODATA_value  -9999
+        """
+        
+        # UTM round Onslow 
+        
+
+
+        ncols = 11  # Nx
+        nrows = 12  # Ny
+        xllcorner = 240000
+        yllcorner = 7620000
+        cellsize  = 6000
+        NODATA_value =  -9999
+        
+        #xllcorner = 0
+        #yllcorner = 100
+        #cellsize  = 10
+        #NODATA_value =  -9999
+        
+        #Create .asc file
+        #txt_file = tempfile.mktemp(".asc")
+        txt_file = 'test_asc.asc'
+        datafile = open(txt_file,"w")
+        datafile.write('ncols '+str(ncols)+"\n")
+        datafile.write('nrows '+str(nrows)+"\n")
+        datafile.write('xllcorner '+str(xllcorner)+"\n")
+        datafile.write('yllcorner '+str(yllcorner)+"\n")
+        datafile.write('cellsize '+str(cellsize)+"\n")
+        datafile.write('NODATA_value '+str(NODATA_value)+"\n")
+        
+        x = num.linspace(xllcorner, xllcorner+(ncols-1)*cellsize, ncols)
+        y = num.linspace(yllcorner, yllcorner+(nrows-1)*cellsize, nrows)
+        points = axes2points(x, y)
+        
+        #print points
+        #print x.shape, x
+        #print y.shape, y
+        
+        datavalues = linear_function(points)
+        #print datavalues 
+        
+        datavalues = datavalues.reshape(nrows,ncols)
+
+        #print datavalues
+        #print datavalues.shape
+        for row in datavalues:
+            #print row
+            datafile.write(" ".join(str(elem) for elem in row) + "\n")         
+        datafile.close()
+
+        #print quantity.vertex_values
+        #print quantity.centroid_values 
+
+        quantity.set_values_from_asc_file(txt_file,
+                             location='vertices',
+                             indices=None,
+                             verbose=False)
+
+        # check order of vertices
+        answer =  [ 23298000. , 23358000. , 23118000. ]
+        answer = [ 23298000. ,  23118000. , 23358000 ]
+        print quantity.vertex_values
+        assert num.allclose(quantity.vertex_values, answer)
+        
+        
+        #print quantity.vertex_values
+        #print quantity.centroid_values 
+        quantity.set_values(0.0)
+        
+        
+        #print quantity.vertex_values
+        #print quantity.centroid_values 
+        quantity.set_values_from_asc_file(txt_file,
+                             location='centroids',
+                             indices=None,
+                             verbose=False)
+        
+        #print quantity.vertex_values
+        #print quantity.centroid_values      
+        
+        answer =  [ 23258000. ]
+        assert num.allclose(quantity.centroid_values, answer)
+        #Cleanup
+        #import os
+        #os.remove(txt_file)        
+        
 
 
@@ -2571 +2722 @@
 
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Quantity, 'test')
+    suite = unittest.makeSuite(Test_Quantity, 'test') #_set_values_from_asc')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)

trunk/anuga_core/source/anuga/anuga_exceptions.py

@@ -18 +18 @@
 class ANUGAError(Exception):
     """ Generic ANUGA error. """
-    def __init__(self, args=None):
-        self.args = args
+    #def __init__(self, args=None):
+    #self.args = args
+    pass
 
 class DataMissingValuesError(exceptions.Exception):

trunk/anuga_core/source/anuga/config.py

@@ -178 +178 @@
 
 # Water depth below which it is considered to be 0 in the model
-minimum_allowed_height = 1.0e-03 
+minimum_allowed_height = 1.0e-05 
 
 # Water depth below which it is *stored* as 0
-minimum_storable_height = 1.0e-05
+minimum_storable_height = 1.0e-03
 
 # FIXME (Ole): Redefine this parameter to control maximal speeds in general

trunk/anuga_core/source/anuga/fit_interpolate/fit.py

@@ -67 +67 @@
         Padarn Note 05/12/12: This documentation should probably
         be updated to account for the fact that the fitting is now
-        done in C. I wasn't sure what details were neccesary though.
+        done in C. I wasn't sure what details were necessary though.
 
         Fit data at points to the vertices of a mesh.

trunk/anuga_core/source/anuga/fit_interpolate/interpolate2d.py

@@ -261 +261 @@
     Nx = len(x)
     Ny = len(y)
-    msg = ('Input array Z must have dimensions %i x %i corresponding to the '
-           'lengths of the input coordinates x and y. However, '
-           'Z has dimensions %i x %i.' % (Nx, Ny, m, n))
+
+    msg =  ('Input array Z must have dimensions %i x %i corresponding to the '
+            'lengths of the input coordinates x and y. However, '
+            'Z has dimensions %i x %i.' % (Nx, Ny, m, n))
     if not (Nx == m and Ny == n):
         raise ANUGAError(msg)

trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -475 +475 @@
 
         A = stage[0,:]
-        #print A[0], (Q2[0,0] + Q1[1,0])/2
-        assert num.allclose(A[0], (Q2[0] + Q1[1])/2)
-        assert num.allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
+        
+#         print stage[0,:]
+#         print A[0], (Q2[0] + Q1[1])/2
+#         print A[1], (Q0[1] + Q1[3] + Q2[2])/3
+#         print A[2], Q0[3]
+#         print A[3], (Q0[0] + Q1[5] + Q2[4])/3
+        assert num.allclose(A[0], (Q2[0] + Q1[1])/2, rtol=1.0e-5, atol=1.0e-5)
+        
+        assert num.allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3, rtol=1.0e-5, atol=1.0e-5)
         assert num.allclose(A[2], Q0[3])
-        assert num.allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
+        assert num.allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3, rtol=1.0e-5, atol=1.0e-5)
 
         #Center point
         assert num.allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
-                                   Q0[5] + Q2[6] + Q1[7])/6)
+                                   Q0[5] + Q2[6] + Q1[7])/6, rtol=1.0e-5, atol=1.0e-5)
         
         fid.close()
@@ -527 +533 @@
 
         A = stage[1,:]
-        assert num.allclose(A[0], (Q2[0] + Q1[1])/2)
-        assert num.allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
-        assert num.allclose(A[2], Q0[3])
-        assert num.allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
+        assert num.allclose(A[0], (Q2[0] + Q1[1])/2, rtol=1.0e-5, atol=1.0e-5)
+        assert num.allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3, rtol=1.0e-5, atol=1.0e-5)
+        assert num.allclose(A[2], Q0[3], rtol=1.0e-5, atol=1.0e-5)
+        assert num.allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3, rtol=1.0e-5, atol=1.0e-5)
 
         #Center point
         assert num.allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
-                                   Q0[5] + Q2[6] + Q1[7])/6)
+                                   Q0[5] + Q2[6] + Q1[7])/6, rtol=1.0e-5, atol=1.0e-5)
 
 
@@ -577 +583 @@
             
             if t == 0.0:
-                assert num.allclose(stage, self.initial_stage)
-                assert num.allclose(stage_file[:], stage.flatten())
+                assert num.allclose(stage, self.initial_stage, rtol=1.0e-5, atol=1.0e-5)
+                assert num.allclose(stage_file[:], stage.flatten(),rtol=1.0e-5, atol=1.0e-5)
             else:
-                assert not num.allclose(stage, self.initial_stage)
-                assert not num.allclose(stage_file[:], stage.flatten())
+                assert not num.allclose(stage, self.initial_stage, rtol=1.0e-5, atol=1.0e-5)
+                assert not num.allclose(stage_file[:], stage.flatten(), rtol=1.0e-5, atol=1.0e-5)
 
             fid.close()