Changeset 1660

Timestamp:

Aug 1, 2005, 11:58:24 AM (19 years ago)

Author:

hitchman

Message:

introduced vector operations

Location:

inundation/ga/storm_surge/pyvolution

Files:

• data_manager.py (modified) (8 diffs)
• test_data_manager.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -1 @@
-"""datamanager.py - input output for AnuGA
-
-
-This module takes care of reading and writing datafiles such as topograhies,
-model output, etc
-
-
-Formats used within AnuGA:
+"""Functions to store and retrieve data for the Shallow Water Wave equation.
+There are two kinds of data
+
+  1: Constant data: Vertex coordinates and field values. Stored once
+  2: Variable data: Conserved quantities. Stored once per timestep.
+
+All data is assumed to reside at vertex locations.
+
+
+Formats used within ANUGA:
 
 .sww: Netcdf format for storing model output
@@ -207 +209 @@
 #Class for storing output to e.g. visualisation
 class Data_format_sww(Data_format):
-    """Interface to native NetCDF format (.sww) for storing model output
-
-    There are two kinds of data
-
-    1: Constant data: Vertex coordinates and field values. Stored once
-    2: Variable data: Conserved quantities. Stored once per timestep.
-
-    All data is assumed to reside at vertex locations.
+    """Interface to native NetCDF format (.sww)
     """
 
@@ -253 +248 @@
             if hasattr(domain, 'texture'):                  
                 fid.texture = domain.texture
-            else:
-                pass
-                #fid.texture = 'None'       
+            else:                           
+                fid.texture = 'None'        
 
             #Reference point
@@ -927 +921 @@
     return cls(domain, mode)
 
-
-def xya2pts(basename_in, basename_out=None, verbose=False,
-            #easting_min=None, easting_max=None,
-            #northing_min=None, northing_max=None,
-            attribute_name = 'elevation'):
-    """Read points data from ascii (.xya)
-
-    Example:
-
-              x(m)        y(m)        z(m)
-         0.00000e+00  0.00000e+00  1.3535000e-01
-         0.00000e+00  1.40000e-02  1.3535000e-01
-
-    
-
-    Convert to NetCDF pts format which is
-
-    points:  (Nx2) Float array
-    elevation: N Float array
-
-    Only lines that contain three numeric values are processed
-    """
-
-    import os 
-    from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import Float, arrayrange, concatenate
-
-    root, ext = os.path.splitext(basename_in)
-
-    if ext == '': ext = '.xya'
-
-    #Get NetCDF
-    infile = open(root + ext, 'r')  #Open existing xya file for read
-
-    if verbose: print 'Reading xya points from %s' %(root + ext)
-
-    points = []
-    attribute = []
-    for i, line in enumerate(infile.readlines()):
-        fields = line.split()
-
-        try:
-            assert len(fields) == 3
-        except:    
-            print 'WARNING: Line %d doesn\'t have 3 elements: %s' %(i, line) 
-
-        try:
-            x = float( fields[0] )
-            y = float( fields[1] )
-            z = float( fields[2] )            
-        except:
-            continue
-
-        points.append( [x, y] )
-        attribute.append(z)
-
-
-    #Get output file
-    if basename_out == None:
-        ptsname = root + '.pts'
-    else:
-        ptsname = basename_out + '.pts'
-
-    if verbose: print 'Store to NetCDF file %s' %ptsname
-    # NetCDF file definition
-    outfile = NetCDFFile(ptsname, 'w')
-
-    #Create new file
-    outfile.institution = 'Geoscience Australia'
-    outfile.description = 'NetCDF pts format for compact and '\
-                          'portable storage of spatial point data'
-    
-
-    #Georeferencing
-    from coordinate_transforms.geo_reference import Geo_reference
-    Geo_reference().write_NetCDF(outfile)
-    
-
-    outfile.createDimension('number_of_points', len(points))
-    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
-
-    # variable definitions
-    outfile.createVariable('points', Float, ('number_of_points',
-                                             'number_of_dimensions'))
-    outfile.createVariable(attribute_name, Float, ('number_of_points',))
-
-    # Get handles to the variables
-    nc_points = outfile.variables['points']
-    nc_attribute = outfile.variables[attribute_name]
-
-    #Store data
-    nc_points[:, :] = points
-    nc_attribute[:] = attribute
-
-    infile.close()
-    outfile.close()
-
-
 def dem2pts(basename_in, basename_out=None, verbose=False,
             easting_min=None, easting_max=None,
@@ -1048 +944 @@
     import os
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import Float, arrayrange, concatenate
+    from Numeric import Float, zeros
 
     root = basename_in
@@ -1117 +1013 @@
 
     # dimension definitions
-    nrows_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
-    ncols_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
+    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
+    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
     outfile.createDimension('number_of_points', nrows_in_bounding_box*ncols_in_bounding_box)
     outfile.createDimension('number_of_dimensions', 2) #This is 2d data
@@ -1132 +1028 @@
 
     #Store data
-    #FIXME: Could perhaps be faster using array operations
-    index = 0
+    #FIXME: Could perhaps be faster using array operations (Fixed 27/7/05)
+    global_index = 0
     for i in range(nrows):
         if verbose: print 'Processing row %d of %d' %(i, nrows)
+        lower_index = global_index
+        tpoints = zeros((ncols_in_bounding_box, 2), Float)
+        telev =  zeros(ncols_in_bounding_box, Float)
+        local_index = 0
 
         y = (nrows-i)*cellsize + yllcorner
         for j in range(ncols):
-            #index = i*ncols + j
 
             x = j*cellsize + xllcorner
             if easting_min <= x <= easting_max and \
                northing_min <= y <= northing_max:
-                points[index, :] = [x-easting_min,y-northing_min]
-                elevation[index] = dem_elevation[i, j]
-                index += 1
-
-    #print easting_min, easting_max
-    #print index, nrows_in_bounding_box*ncols_in_bounding_box
-    #print nrows, ncols
-    assert index == nrows_in_bounding_box*ncols_in_bounding_box, 'index not equal to number of points'
+                tpoints[local_index, :] = [x-easting_min,y-northing_min]
+                telev[local_index] = dem_elevation[i, j]
+                global_index += 1
+                local_index += 1
+
+        upper_index = global_index
+
+        if upper_index == lower_index + ncols_in_bounding_box:
+            points[lower_index:upper_index, :] = tpoints
+            elevation[lower_index:upper_index] = telev
+
+    assert global_index == nrows_in_bounding_box*ncols_in_bounding_box, 'index not equal to number of points'
 
     infile.close()
@@ -2379 +2282 @@
 
     volumes = fid.variables['volumes'] #Connectivity
+
+

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -1819 +1819 @@
 
 
+
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Data_Manager,'test')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2domain')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
     runner = unittest.TextTestRunner()
     runner.run(suite)