Changeset 6157 for anuga_core/source/anuga/shallow_water/data_manager.py

Timestamp:

Jan 13, 2009, 5:42:30 PM (15 years ago)

Author:

rwilson

Message:

Change Numeric imports to general form - ready to change to NumPy?.

File:

• anuga_core/source/anuga/shallow_water/data_manager.py (modified) (119 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -61 +61 @@
 from os import sep, path, remove, mkdir, access, F_OK, W_OK, getcwd
 
-from Numeric import concatenate, array, Float, Int, Int32, resize, \
-     sometrue, searchsorted, zeros, allclose, around, reshape, \
-     transpose, sort, NewAxis, ArrayType, compress, take, arange, \
-     argmax, alltrue, shape, Float32, size
+import Numeric as num
 
 from Scientific.IO.NetCDF import NetCDFFile
@@ -345 +342 @@
     def __init__(self, domain, mode=netcdf_mode_w, max_size=2000000000, recursion=False):
         from Scientific.IO.NetCDF import NetCDFFile
-        from Numeric import Int, Float, Float32
-
-        self.precision = Float32 #Use single precision for quantities
+
+        self.precision = num.Float32 #Use single precision for quantities
         self.recursion = recursion
         self.mode = mode
@@ -425 +421 @@
 
         from Scientific.IO.NetCDF import NetCDFFile
-        from Numeric import concatenate, Int
 
         domain = self.domain
@@ -443 +438 @@
 
         # store the connectivity data
-        points = concatenate( (X[:,NewAxis],Y[:,NewAxis]), axis=1 )
+        points = num.concatenate( (X[:,num.NewAxis],Y[:,num.NewAxis]), axis=1 )
         self.writer.store_triangulation(fid,
                                         points,
 #                                        V.astype(volumes.typecode()),
-                                        V.astype(Float32),
+                                        V.astype(num.Float32),
                                         Z,
                                         points_georeference=\
@@ -466 +461 @@
         from time import sleep
         from os import stat
-        from Numeric import choose
 
         if names is None:
@@ -563 +557 @@
 
                 storable_indices = (A-z[:] >= self.minimum_storable_height)
-                stage = choose(storable_indices, (z[:], A))
+                stage = num.choose(storable_indices, (z[:], A))
 
                 # Define a zero vector of same size and type as A
                 # for use with momenta
-                null = zeros(size(A), A.typecode())
+                null = num.zeros(num.size(A), A.typecode())
 
                 # Get xmomentum where depth exceeds minimum_storable_height
@@ -573 +567 @@
                 xmom, _ = Q.get_vertex_values(xy=False,
                                               precision=self.precision)
-                xmomentum = choose(storable_indices, (null, xmom))
+                xmomentum = num.choose(storable_indices, (null, xmom))
 
 
@@ -580 +574 @@
                 ymom, _ = Q.get_vertex_values(xy=False,
                                               precision=self.precision)
-                ymomentum = choose(storable_indices, (null, ymom))
+                ymomentum = num.choose(storable_indices, (null, ymom))
 
                 # Write quantities to underlying data  file
@@ -627 +621 @@
     def __init__(self, domain, mode=netcdf_mode_w):
         from Scientific.IO.NetCDF import NetCDFFile
-        from Numeric import Int, Float, Float
-
-        self.precision = Float #Use full precision
+
+        self.precision = num.Float #Use full precision
 
         Data_format.__init__(self, domain, 'sww', mode)
@@ -657 +650 @@
 
 
-            fid.createVariable('volumes', Int, ('number_of_volumes',
-                                                'number_of_vertices'))
+            fid.createVariable('volumes', num.Int, ('number_of_volumes',
+                                                    'number_of_vertices'))
 
             fid.createVariable('time', self.precision, ('number_of_timesteps',))
@@ -680 +673 @@
 
         from Scientific.IO.NetCDF import NetCDFFile
-        from Numeric import concatenate
 
         domain = self.domain
@@ -1151 +1143 @@
     Y = {}
     for key in X.keys():
-        Y[key] = array([float(x) for x in X[key]])
+        Y[key] = num.array([float(x) for x in X[key]])
 
     return Y
@@ -1269 +1261 @@
 
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import Float, zeros
 
     # Get NetCDF
@@ -1284 +1275 @@
 
     M = size  #Number of lines
-    xx = zeros((M,3), Float)
-    yy = zeros((M,3), Float)
-    zz = zeros((M,3), Float)
+    xx = num.zeros((M,3), num.Float)
+    yy = num.zeros((M,3), num.Float)
+    zz = num.zeros((M,3), num.Float)
 
     for i in range(M):
@@ -1324 +1315 @@
     import glob, os
     from anuga.config import data_dir
-    from Numeric import zeros, Float
 
     # Get bathymetry and x,y's
@@ -1330 +1320 @@
 
     M = len(lines)  #Number of lines
-    x = zeros((M,3), Float)
-    y = zeros((M,3), Float)
-    z = zeros((M,3), Float)
+    x = num.zeros((M,3), num.Float)
+    y = num.zeros((M,3), num.Float)
+    z = num.zeros((M,3), num.Float)
 
     for i, line in enumerate(lines):
@@ -1524 +1514 @@
     import os
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import Float, zeros, reshape, sum
 
     root = basename_in
@@ -1591 +1580 @@
     outfile.nrows = nrows
 
-    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
+    dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols))
     totalnopoints = nrows*ncols
 
@@ -1640 +1629 @@
 
     # Variable definitions
-    outfile.createVariable('points', Float, ('number_of_points',
-                                             'number_of_dimensions'))
-    outfile.createVariable('elevation', Float, ('number_of_points',))
+    outfile.createVariable('points', num.Float, ('number_of_points',
+                                                 'number_of_dimensions'))
+    outfile.createVariable('elevation', num.Float, ('number_of_points',))
 
     # Get handles to the variables
@@ -1660 +1649 @@
 
         v = dem_elevation_r[i,index1:index2+1]
-        no_NODATA = sum(v == NODATA_value)
+        no_NODATA = num.sum(v == NODATA_value)
         if no_NODATA > 0:
             newcols = lenv - no_NODATA  # ncols_in_bounding_box - no_NODATA
@@ -1666 +1655 @@
             newcols = lenv              # ncols_in_bounding_box
 
-        telev = zeros(newcols, Float)
-        tpoints = zeros((newcols, 2), Float)
+        telev = num.zeros(newcols, num.Float)
+        tpoints = num.zeros((newcols, 2), num.Float)
 
         local_index = 0
@@ -1789 +1778 @@
     import os
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import Float, zeros, reshape
 
     root = basename_in
@@ -2133 +2121 @@
 
     import sys
-    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape
-    from Numeric import array2string, sometrue
 
     from anuga.utilities.polygon import inside_polygon, outside_polygon, \
@@ -2254 +2240 @@
         #q has a time component, must be reduced alongthe temporal dimension
         if verbose: print 'Reducing quantity %s' %quantity
-        q_reduced = zeros(number_of_points, Float)
+        q_reduced = num.zeros(number_of_points, num.Float)
 
         if timestep is not None:
@@ -2314 +2300 @@
     y = y + yllcorner - newyllcorner
 
-    vertex_points = concatenate ((x[:,NewAxis], y[:,NewAxis]), axis=1)
+    vertex_points = num.concatenate ((x[:,num.NewAxis], y[:,num.NewAxis]), axis=1)
     assert len(vertex_points.shape) == 2
 
-    grid_points = zeros ((ncols*nrows, 2), Float)
+    grid_points = num.zeros ((ncols*nrows, 2), num.Float)
 
     for i in xrange(nrows):
@@ -2359 +2345 @@
     if format.lower() == 'ers':
         # setup ERS header information
-        grid_values = reshape(grid_values, (nrows, ncols))
+        grid_values = num.reshape(grid_values, (nrows, ncols))
         header = {}
         header['datum'] = '"' + datum + '"'
@@ -2426 +2412 @@
             slice = grid_values[base_index:base_index+ncols]
             #s = array2string(slice, max_line_width=sys.maxint)
-            s = array2string(slice, max_line_width=sys.maxint,
-                             precision=number_of_decimal_places)
+            s = num.array2string(slice, max_line_width=sys.maxint,
+                                 precision=number_of_decimal_places)
             ascid.write(s[1:-1] + '\n')
 
@@ -2543 +2529 @@
 
     import sys
-    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape
-    from Numeric import array2string, sometrue
     from anuga.utilities.polygon import inside_polygon, outside_polygon, \
              separate_points_by_polygon
@@ -2631 +2615 @@
         if verbose: print 'Reducing quantity %s' % quantity
 
-        q_reduced = zeros(number_of_points, Float)
+        q_reduced = num.zeros(number_of_points, num.Float)
         for k in range(number_of_points):
             q_reduced[k] = reduction(q[:,k])
@@ -2645 +2629 @@
 
     # Create grid and update xll/yll corner and x,y
-    vertex_points = concatenate((x[:, NewAxis], y[:, NewAxis]), axis=1)
+    vertex_points = num.concatenate((x[:, num.NewAxis], y[:, num.NewAxis]), axis=1)
     assert len(vertex_points.shape) == 2
 
@@ -2748 +2732 @@
     import os
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import Float, array
 
     root = basename_in
@@ -2866 +2849 @@
 
     # variable definitions
-    fid.createVariable('elevation', Float, ('number_of_rows',
-                                            'number_of_columns'))
+    fid.createVariable('elevation', num.Float, ('number_of_rows',
+                                                'number_of_columns'))
 
     # Get handles to the variables
@@ -2878 +2861 @@
         if verbose and i % ((n+10)/10) == 0:
             print 'Processing row %d of %d' % (i, nrows)
-        elevation[i, :] = array([float(x) for x in fields])
+        elevation[i, :] = num.array([float(x) for x in fields])
 
     fid.close()
@@ -2946 +2929 @@
     import os
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import Float, Int, Int32, searchsorted, zeros, array
-    from Numeric import allclose, around
-
-    precision = Float
+
+    precision = num.Float
 
     msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
@@ -3022 +3003 @@
 
     # Precision used by most for lat/lon is 4 or 5 decimals
-    e_lat = around(file_e.variables[dim_e_latitude][:], 5)
-    e_lon = around(file_e.variables[dim_e_longitude][:], 5)
+    e_lat = num.around(file_e.variables[dim_e_latitude][:], 5)
+    e_lon = num.around(file_e.variables[dim_e_longitude][:], 5)
 
     # Check that files are compatible
-    assert allclose(latitudes, file_u.variables[dim_u_latitude])
-    assert allclose(latitudes, file_v.variables[dim_v_latitude])
-    assert allclose(latitudes, e_lat)
-
-    assert allclose(longitudes, file_u.variables[dim_u_longitude])
-    assert allclose(longitudes, file_v.variables[dim_v_longitude])
-    assert allclose(longitudes, e_lon)
+    assert num.allclose(latitudes, file_u.variables[dim_u_latitude])
+    assert num.allclose(latitudes, file_v.variables[dim_v_latitude])
+    assert num.allclose(latitudes, e_lat)
+
+    assert num.allclose(longitudes, file_u.variables[dim_u_longitude])
+    assert num.allclose(longitudes, file_v.variables[dim_v_longitude])
+    assert num.allclose(longitudes, e_lon)
 
     if mint is None:
@@ -3038 +3019 @@
         mint = times[0]
     else:
-        jmin = searchsorted(times, mint)
+        jmin = num.searchsorted(times, mint)
 
     if maxt is None:
@@ -3044 +3025 @@
         maxt = times[-1]
     else:
-        jmax = searchsorted(times, maxt)
+        jmax = num.searchsorted(times, maxt)
 
     kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],
@@ -3089 +3070 @@
 
     #Cleanup
-    from Numeric import sometrue
-
     missing = (amplitudes == nan_ha)
-    if sometrue (missing):
+    if num.sometrue (missing):
         if fail_on_NaN:
             msg = 'NetCDFFile %s contains missing values' \
@@ -3101 +3080 @@
 
     missing = (uspeed == nan_ua)
-    if sometrue (missing):
+    if num.sometrue (missing):
         if fail_on_NaN:
             msg = 'NetCDFFile %s contains missing values' \
@@ -3110 +3089 @@
 
     missing = (vspeed == nan_va)
-    if sometrue (missing):
+    if num.sometrue (missing):
         if fail_on_NaN:
             msg = 'NetCDFFile %s contains missing values' \
@@ -3119 +3098 @@
 
     missing = (elevations == nan_e)
-    if sometrue (missing):
+    if num.sometrue (missing):
         if fail_on_NaN:
             msg = 'NetCDFFile %s contains missing values' \
@@ -3187 +3166 @@
     sww.store_header(outfile, times, number_of_volumes,
                      number_of_points, description=description,
-                     verbose=verbose, sww_precision=Float)
+                     verbose=verbose, sww_precision=num.Float)
 
     # Store
     from anuga.coordinate_transforms.redfearn import redfearn
-    x = zeros(number_of_points, Float)  #Easting
-    y = zeros(number_of_points, Float)  #Northing
+    x = num.zeros(number_of_points, num.Float)  #Easting
+    y = num.zeros(number_of_points, num.Float)  #Northing
 
     if verbose: print 'Making triangular grid'
@@ -3226 +3205 @@
             volumes.append([v4,v3,v2]) #Lower element
 
-    volumes = array(volumes)
+    volumes = num.array(volumes)
 
     if origin is None:
@@ -3242 +3221 @@
 
     #FIXME use the Write_sww instance(sww) to write this info
-    from Numeric import resize
-    z = resize(z, outfile.variables['z'][:].shape)
+    z = num.resize(z, outfile.variables['z'][:].shape)
     outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
     outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
     outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat.
     outfile.variables['elevation'][:] = z
-    outfile.variables['volumes'][:] = volumes.astype(Int32) #For Opteron 64
+    outfile.variables['volumes'][:] = volumes.astype(num.Int32) #For Opteron 64
 
     #Time stepping
@@ -3331 +3309 @@
 
     import time, calendar
-    from Numeric import array
     from anuga.config import time_format
     from anuga.utilities.numerical_tools import ensure_numeric
@@ -3371 +3348 @@
 
     # Read times proper
-    from Numeric import zeros, Float, alltrue
     from anuga.config import time_format
     import time, calendar
@@ -3382 +3358 @@
     d = len(q)
 
-    T = zeros(N, Float)       # Time
-    Q = zeros((N, d), Float)  # Values
+    T = num.zeros(N, num.Float)       # Time
+    Q = num.zeros((N, d), num.Float)  # Values
 
     for i, line in enumerate(lines):
@@ -3398 +3374 @@
     msg = 'File %s must list time as a monotonuosly ' % filename
     msg += 'increasing sequence'
-    assert alltrue(T[1:] - T[:-1] > 0), msg
+    assert num.alltrue(T[1:] - T[:-1] > 0), msg
 
     #Create NetCDF file
@@ -3419 +3395 @@
     fid.createDimension('number_of_timesteps', len(T))
 
-    fid.createVariable('time', Float, ('number_of_timesteps',))
+    fid.createVariable('time', num.Float, ('number_of_timesteps',))
 
     fid.variables['time'][:] = T
@@ -3429 +3405 @@
             name = 'Attribute%d' % i
 
-        fid.createVariable(name, Float, ('number_of_timesteps',))
+        fid.createVariable(name, num.Float, ('number_of_timesteps',))
         fid.variables[name][:] = Q[:,i]
 
@@ -3488 +3464 @@
     from Scientific.IO.NetCDF import NetCDFFile
     from shallow_water import Domain
-    from Numeric import asarray, transpose, resize
 
     # initialise NaN.
@@ -3511 +3486 @@
     starttime = fid.starttime[0]
     volumes = fid.variables['volumes'][:]       # Connectivity
-    coordinates = transpose(asarray([x.tolist(), y.tolist()]))
+    coordinates = num.transpose(num.asarray([x.tolist(), y.tolist()]))
     # FIXME (Ole): Something like this might be better:
     #                 concatenate((x, y), axis=1)
-    # or              concatenate((x[:,NewAxis], x[:,NewAxis]), axis=1)
+    # or              concatenate((x[:,num.NewAxis], x[:,num.NewAxis]), axis=1)
 
     conserved_quantities = []
@@ -3608 +3583 @@
                 X = (X*data) + (data==0)*NaN_filler
         if unique:
-            X = resize(X, (len(X)/3, 3))
+            X = num.resize(X, (len(X)/3, 3))
         domain.set_quantity(quantity, X)
     #
@@ -3633 +3608 @@
                 X = (X*data) + (data==0)*NaN_filler
         if unique:
-            X = resize(X, (X.shape[0]/3, 3))
+            X = num.resize(X, (X.shape[0]/3, 3))
         domain.set_quantity(quantity, X)
 
@@ -3702 +3677 @@
 # @param boundary 
 def weed(coordinates, volumes, boundary=None):
-    if type(coordinates) == ArrayType:
+    if type(coordinates) == num.ArrayType:
         coordinates = coordinates.tolist()
-    if type(volumes) == ArrayType:
+    if type(volumes) == num.ArrayType:
         volumes = volumes.tolist()
 
@@ -3781 +3756 @@
     import os
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import Float, zeros, sum, reshape, equal
 
     root = basename_in
@@ -3851 +3825 @@
 
     # variable definition
-    outfile.createVariable('elevation', Float, ('number_of_points',))
+    outfile.createVariable('elevation', num.Float, ('number_of_points',))
 
     # Get handle to the variable
     elevation = outfile.variables['elevation']
 
-    dem_elevation_r = reshape(dem_elevation, (nrows, ncols))
+    dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols))
 
     #Store data
@@ -3864 +3838 @@
 
         lower_index = global_index
-        telev = zeros(ncols_new, Float)
+        telev = num.zeros(ncols_new, num.Float)
         local_index = 0
         trow = i * cellsize_ratio
@@ -3876 +3850 @@
             #decimated dem to NODATA_value, else compute decimated
             #value using stencil
-            if sum(sum(equal(tmp, NODATA_value))) > 0:
+            if num.sum(num.sum(num.equal(tmp, NODATA_value))) > 0:
                 telev[local_index] = NODATA_value
             else:
-                telev[local_index] = sum(sum(tmp * stencil))
+                telev[local_index] = num.sum(num.sum(tmp * stencil))
 
             global_index += 1
@@ -3992 +3966 @@
     from anuga.coordinate_transforms.redfearn import redfearn
 
-    precision = Float # So if we want to change the precision its done here
+    precision = num.Float # So if we want to change the precision its done here
 
     # go in to the bath dir and load the only file,
@@ -4093 +4067 @@
     #################################
 
-    outfile.createVariable('volumes', Int, ('number_of_volumes',
-                                            'number_of_vertices'))
+    outfile.createVariable('volumes', num.Int, ('number_of_volumes',
+                                                'number_of_vertices'))
 
     outfile.createVariable('time', precision, ('number_of_timesteps',))
@@ -4110 +4084 @@
     from anuga.coordinate_transforms.redfearn import redfearn
 
-    x = zeros(number_of_points, Float)  #Easting
-    y = zeros(number_of_points, Float)  #Northing
+    x = num.zeros(number_of_points, num.Float)  #Easting
+    y = num.zeros(number_of_points, num.Float)  #Northing
 
     if verbose: print 'Making triangular grid'
@@ -4147 +4121 @@
             volumes.append([v4,v2,v3]) #Lower element
 
-    volumes = array(volumes)
+    volumes = num.array(volumes)
 
     geo_ref = Geo_reference(refzone, min(x), min(y))
@@ -4165 +4139 @@
         print 'geo_ref: ', geo_ref
 
-    z = resize(bath_grid,outfile.variables['z'][:].shape)
+    z = num.resize(bath_grid,outfile.variables['z'][:].shape)
     outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
     outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
@@ -4172 +4146 @@
     outfile.variables['z'][:] = z
     outfile.variables['elevation'][:] = z
-    outfile.variables['volumes'][:] = volumes.astype(Int32) # On Opteron 64
+    outfile.variables['volumes'][:] = volumes.astype(num.Int32) # On Opteron 64
 
     stage = outfile.variables['stage']
@@ -4198 +4172 @@
         # handle missing values
         missing = (elevation_grid == elevation_meta['NODATA_value'])
-        if sometrue (missing):
+        if num.sometrue (missing):
             if fail_on_NaN:
                 msg = 'File %s contains missing values' \
@@ -4254 +4228 @@
     longitudes = ensure_numeric(longitudes)
 
-    assert allclose(sort(longitudes), longitudes)
+    assert num.allclose(num.sort(longitudes), longitudes)
 
     #print latitudes[0],longitudes[0]
@@ -4261 +4235 @@
 
     lat_ascending = True
-    if not allclose(sort(latitudes), latitudes):
+    if not num.allclose(num.sort(latitudes), latitudes):
         lat_ascending = False
         # reverse order of lat, so it's in ascending order
         latitudes = latitudes[::-1]
-        assert allclose(sort(latitudes), latitudes)
+        assert num.allclose(num.sort(latitudes), latitudes)
 
     largest_lat_index = len(latitudes)-1
@@ -4273 +4247 @@
         lat_min_index = 0
     else:
-        lat_min_index = searchsorted(latitudes, minlat)-1
+        lat_min_index = num.searchsorted(latitudes, minlat)-1
         if lat_min_index <0:
             lat_min_index = 0
@@ -4280 +4254 @@
         lat_max_index = largest_lat_index #len(latitudes)
     else:
-        lat_max_index = searchsorted(latitudes, maxlat)
+        lat_max_index = num.searchsorted(latitudes, maxlat)
         if lat_max_index > largest_lat_index:
             lat_max_index = largest_lat_index
@@ -4287 +4261 @@
         lon_min_index = 0
     else:
-        lon_min_index = searchsorted(longitudes, minlon)-1
+        lon_min_index = num.searchsorted(longitudes, minlon)-1
         if lon_min_index <0:
             lon_min_index = 0
@@ -4294 +4268 @@
         lon_max_index = len(longitudes)
     else:
-        lon_max_index = searchsorted(longitudes, maxlon)
+        lon_max_index = num.searchsorted(longitudes, maxlon)
 
     # Reversing the indexes, if the lat array is decending
@@ -4349 +4323 @@
         cells = line.split()
         assert len(cells) == ncols
-        grid.append(array([float(x) for x in cells]))
-    grid = array(grid)
+        grid.append(num.array([float(x) for x in cells]))
+    grid = num.array(grid)
 
     return {'xllcorner':xllcorner,
@@ -4364 +4338 @@
 lat_name = 'LAT'
 time_name = 'TIME'
-precision = Float # So if we want to change the precision its done here
+precision = num.Float # So if we want to change the precision its done here
 
 ##
@@ -4645 +4619 @@
     lonlatdep = p_array.array('f')
     lonlatdep.read(mux_file, columns * points_num)
-    lonlatdep = array(lonlatdep, typecode=Float)
-    lonlatdep = reshape(lonlatdep, (points_num, columns))
+    lonlatdep = num.array(lonlatdep, typecode=num.Float)
+    lonlatdep = num.reshape(lonlatdep, (points_num, columns))
 
     lon, lat, depth = lon_lat2grid(lonlatdep)
@@ -4674 +4648 @@
         hz_p_array = p_array.array('f')
         hz_p_array.read(mux_file, points_num)
-        hz_p = array(hz_p_array, typecode=Float)
-        hz_p = reshape(hz_p, (len(lon), len(lat)))
-        hz_p = transpose(hz_p)  # mux has lat varying fastest, nc has long v.f.
+        hz_p = num.array(hz_p_array, typecode=num.Float)
+        hz_p = num.reshape(hz_p, (len(lon), len(lat)))
+        hz_p = num.transpose(hz_p)  # mux has lat varying fastest, nc has long v.f.
 
         #write time slice to nc file
@@ -4711 +4685 @@
     outfile.variables[lat_name][:] = ensure_numeric(lat)
 
-    depth = reshape(depth_vector, (len(lat), len(lon)))
+    depth = num.reshape(depth_vector, (len(lat), len(lon)))
     outfile.variables[zname][:] = depth
 
@@ -4772 +4746 @@
     QUANTITY = 2
 
-    long_lat_dep = ensure_numeric(long_lat_dep, Float)
+    long_lat_dep = ensure_numeric(long_lat_dep, num.Float)
 
     num_points = long_lat_dep.shape[0]
@@ -4799 +4773 @@
         last = first + lenlat
 
-        assert allclose(long_lat_dep[first:last,LAT], lat), msg
-        assert allclose(long_lat_dep[first:last,LONG], long[i]), msg
+        assert num.allclose(long_lat_dep[first:last,LAT], lat), msg
+        assert num.allclose(long_lat_dep[first:last,LONG], long[i]), msg
 
     msg = 'Out of range latitudes/longitudes'
@@ -4810 +4784 @@
     # FIXME - make this faster/do this a better way
     # use numeric transpose, after reshaping the quantity vector
-    quantity = zeros(num_points, Float)
+    quantity = num.zeros(num_points, num.Float)
 
     for lat_i, _ in enumerate(lat):
@@ -5282 +5256 @@
     sww = Write_sww()
     sww.store_header(outfile, times, len(volumes), len(points_utm),
-                     verbose=verbose, sww_precision=Float)
+                     verbose=verbose, sww_precision=num.Float)
     outfile.mean_stage = mean_stage
     outfile.zscale = zscale
@@ -5306 +5280 @@
                                  xmomentum=xmomentum,
                                  ymomentum=ymomentum,
-                                 sww_precision=Float)
+                                 sww_precision=num.Float)
         j += 1
 
@@ -5349 +5323 @@
     """
 
-    from Numeric import ones,Float,compress,zeros,arange
     from urs_ext import read_mux2
 
     numSrc = len(filenames)
 
-    file_params = -1 * ones(3, Float)                    # [nsta,dt,nt]
+    file_params = -1 * num.ones(3, num.Float)                    # [nsta,dt,nt]
 
     # Convert verbose to int C flag
@@ -5363 +5336 @@
 
     if weights is None:
-        weights = ones(numSrc)
+        weights = num.ones(numSrc)
 
     if permutation is None:
-        permutation = ensure_numeric([], Float)
+        permutation = ensure_numeric([], num.Float)
 
     # Call underlying C implementation urs2sts_ext.c
@@ -5373 +5346 @@
 
     msg = 'File parameter values were not read in correctly from c file'
-    assert len(compress(file_params > 0, file_params)) != 0, msg
+    assert len(num.compress(file_params > 0, file_params)) != 0, msg
 
     msg = 'The number of stations specifed in the c array and in the file ' \
@@ -5413 +5386 @@
     parameters_index = data.shape[1] - OFFSET
 
-    times = dt * arange(parameters_index)
-    latitudes = zeros(number_of_selected_stations, Float)
-    longitudes = zeros(number_of_selected_stations, Float)
-    elevation = zeros(number_of_selected_stations, Float)
-    quantity = zeros((number_of_selected_stations, parameters_index), Float)
+    times = dt * num.arange(parameters_index)
+    latitudes = num.zeros(number_of_selected_stations, num.Float)
+    longitudes = num.zeros(number_of_selected_stations, num.Float)
+    elevation = num.zeros(number_of_selected_stations, num.Float)
+    quantity = num.zeros((number_of_selected_stations, parameters_index), num.Float)
 
     starttime = 1e16
@@ -5444 +5417 @@
         mux_times_start_i = 0
     else:
-        mux_times_start_i = searchsorted(mux_times, mint)
+        mux_times_start_i = num.searchsorted(mux_times, mint)
 
     if maxt == None:
@@ -5451 +5424 @@
         maxt += 0.5 # so if you specify a time where there is
                     # data that time will be included
-        mux_times_fin_i = searchsorted(mux_times, maxt)
+        mux_times_fin_i = num.searchsorted(mux_times, maxt)
 
     return mux_times_start_i, mux_times_fin_i
@@ -5518 +5491 @@
     import os
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import Float, Int, Int32, searchsorted, zeros, array
-    from Numeric import allclose, around,ones,Float
     from types import ListType,StringType
     from operator import __and__
@@ -5543 +5514 @@
     if weights is None:
         # Default is equal weighting
-        weights = ones(numSrc, Float) / numSrc
+        weights = num.ones(numSrc, num.Float) / numSrc
     else:
         weights = ensure_numeric(weights)
@@ -5622 +5593 @@
             msg = '%s, %s and %s have inconsistent gauge data' \
                   % (files_in[0], files_in[1], files_in[2])
-            assert allclose(times, times_old), msg
-            assert allclose(latitudes, latitudes_old), msg
-            assert allclose(longitudes, longitudes_old), msg
-            assert allclose(elevation, elevation_old), msg
-            assert allclose(starttime, starttime_old), msg
+            assert num.allclose(times, times_old), msg
+            assert num.allclose(latitudes, latitudes_old), msg
+            assert num.allclose(longitudes, longitudes_old), msg
+            assert num.allclose(elevation, elevation_old), msg
+            assert num.allclose(starttime, starttime_old), msg
         times_old = times
         latitudes_old = latitudes
@@ -5665 +5636 @@
     # 0 to number_of_points-1
     if permutation is None:
-        permutation = arange(number_of_points, typecode=Int)
+        permutation = num.arange(number_of_points, typecode=num.Int)
 
     # NetCDF file definition
@@ -5679 +5650 @@
                      description=description,
                      verbose=verbose,
-                     sts_precision=Float)
+                     sts_precision=num.Float)
 
     # Store
     from anuga.coordinate_transforms.redfearn import redfearn
 
-    x = zeros(number_of_points, Float)  # Easting
-    y = zeros(number_of_points, Float)  # Northing
+    x = num.zeros(number_of_points, num.Float)  # Easting
+    y = num.zeros(number_of_points, num.Float)  # Northing
 
     # Check zone boundaries
@@ -5715 +5686 @@
     geo_ref = write_NetCDF_georeference(origin, outfile)
 
-    elevation = resize(elevation, outfile.variables['elevation'][:].shape)
-    outfile.variables['permutation'][:] = permutation.astype(Int32) # Opteron 64
+    elevation = num.resize(elevation, outfile.variables['elevation'][:].shape)
+    outfile.variables['permutation'][:] = permutation.astype(num.Int32) # Opteron 64
     outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
     outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
@@ -5782 +5753 @@
     y = fid.variables['y'][:] + yllcorner
 
-    x = reshape(x, (len(x),1))
-    y = reshape(y, (len(y),1))
-    sts_points = concatenate((x,y), axis=1)
+    x = num.reshape(x, (len(x),1))
+    y = num.reshape(y, (len(y),1))
+    sts_points = num.concatenate((x,y), axis=1)
 
     return sts_points.tolist()
@@ -5832 +5803 @@
                      smoothing=True,
                      order=1,
-                     sww_precision=Float32,
+                     sww_precision=num.Float32,
                      verbose=False):
         """Write an SWW file header.
@@ -5865 +5836 @@
         # This is being used to seperate one number from a list.
         # what it is actually doing is sorting lists from numeric arrays.
-        if type(times) is list or type(times) is ArrayType:
+        if type(times) is list or type(times) is num.ArrayType:
             number_of_times = len(times)
             times = ensure_numeric(times)
@@ -5914 +5885 @@
         outfile.createVariable('z', sww_precision, ('number_of_points',))
 
-        outfile.createVariable('volumes', Int, ('number_of_volumes',
-                                                'number_of_vertices'))
+        outfile.createVariable('volumes', num.Int, ('number_of_volumes',
+                                                    'number_of_vertices'))
 
         # Doing sww_precision instead of Float gives cast errors.
-        outfile.createVariable('time', Float,
+        outfile.createVariable('time', num.Float,
                                ('number_of_timesteps',))
 
@@ -5932 +5903 @@
             outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
 
-        if type(times) is list or type(times) is ArrayType:
+        if type(times) is list or type(times) is num.ArrayType:
             outfile.variables['time'][:] = times    #Store time relative
 
@@ -5987 +5958 @@
 
         number_of_points = len(points_utm)
-        volumes = array(volumes)
-        points_utm = array(points_utm)
+        volumes = num.array(volumes)
+        points_utm = num.array(points_utm)
 
         # given the two geo_refs and the points, do the stuff
@@ -6035 +6006 @@
         outfile.variables['z'][:] = elevation
         outfile.variables['elevation'][:] = elevation  #FIXME HACK
-        outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
+        outfile.variables['volumes'][:] = volumes.astype(num.Int32) #On Opteron 64
 
         q = 'elevation'
@@ -6051 +6022 @@
     # @param verbose True if this function is to be verbose.
     # @param **quant
-    def store_quantities(self, outfile, sww_precision=Float32,
+    def store_quantities(self, outfile, sww_precision=num.Float32,
                          slice_index=None, time=None,
                          verbose=False, **quant):
@@ -6244 +6215 @@
                      number_of_points,
                      description='Converted from URS mux2 format',
-                     sts_precision=Float32,
+                     sts_precision=num.Float32,
                      verbose=False):
         """
@@ -6267 +6238 @@
         # This is being used to seperate one number from a list.
         # what it is actually doing is sorting lists from numeric arrays.
-        if type(times) is list or type(times) is ArrayType:
+        if type(times) is list or type(times) is num.ArrayType:
             number_of_times = len(times)
             times = ensure_numeric(times)
@@ -6287 +6258 @@
 
         # Variable definitions
-        outfile.createVariable('permutation', Int, ('number_of_points',))
+        outfile.createVariable('permutation', num.Int, ('number_of_points',))
         outfile.createVariable('x', sts_precision, ('number_of_points',))
         outfile.createVariable('y', sts_precision, ('number_of_points',))
@@ -6302 +6273 @@
 
         # Doing sts_precision instead of Float gives cast errors.
-        outfile.createVariable('time', Float, ('number_of_timesteps',))
+        outfile.createVariable('time', num.Float, ('number_of_timesteps',))
 
         for q in Write_sts.sts_quantities:
@@ -6314 +6285 @@
             outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
 
-        if type(times) is list or type(times) is ArrayType:
+        if type(times) is list or type(times) is num.ArrayType:
             outfile.variables['time'][:] = times    #Store time relative
 
@@ -6365 +6336 @@
 
         number_of_points = len(points_utm)
-        points_utm = array(points_utm)
+        points_utm = num.array(points_utm)
 
         # given the two geo_refs and the points, do the stuff
@@ -6423 +6394 @@
     # @param verboseTrue if this function is to be verbose.
     # @param **quant Extra keyword args.
-    def store_quantities(self, outfile, sts_precision=Float32,
+    def store_quantities(self, outfile, sts_precision=num.Float32,
                          slice_index=None, time=None,
                          verbose=False, **quant):
@@ -6514 +6485 @@
         lonlatdep = p_array.array('f')
         lonlatdep.read(mux_file, columns * self.points_num)
-        lonlatdep = array(lonlatdep, typecode=Float)
-        lonlatdep = reshape(lonlatdep, (self.points_num, columns))
+        lonlatdep = num.array(lonlatdep, typecode=num.Float)
+        lonlatdep = num.reshape(lonlatdep, (self.points_num, columns))
         self.lonlatdep = lonlatdep
 
@@ -6550 +6521 @@
         hz_p_array = p_array.array('f')
         hz_p_array.read(self.mux_file, self.points_num)
-        hz_p = array(hz_p_array, typecode=Float)
+        hz_p = num.array(hz_p_array, typecode=num.Float)
         self.iter_time_step += 1
 
@@ -6695 +6666 @@
     # array to store data, number in there is to allow float...
     # i'm sure there is a better way!
-    data = array([], typecode=Float)
-    data = resize(data, ((len(lines)-1), len(header_fields)))
+    data = num.array([], typecode=num.Float)
+    data = num.resize(data, ((len(lines)-1), len(header_fields)))
 
     array_number = 0
@@ -6853 +6824 @@
     from Scientific.IO.NetCDF import NetCDFFile
     from shallow_water import Domain
-    from Numeric import asarray, transpose, resize
     from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
 
@@ -6871 +6841 @@
 
     # Mesh (nodes (Mx2), triangles (Nx3))
-    nodes = concatenate((x[:,NewAxis], y[:,NewAxis]), axis=1)
+    nodes = num.concatenate((x[:,num.NewAxis], y[:,num.NewAxis]), axis=1)
     triangles = fid.variables['volumes'][:]
 
@@ -7294 +7264 @@
 
         # Get the relevant quantities (Convert from single precison)
-        elevation = array(fid.variables['elevation'][:], Float)
-        stage = array(fid.variables['stage'][:], Float)
+        elevation = num.array(fid.variables['elevation'][:], num.Float)
+        stage = num.array(fid.variables['stage'][:], num.Float)
 
         # Here's where one could convert nodal information to centroid
@@ -7312 +7282 @@
             # and call it here
 
-            points = concatenate((x[:,NewAxis], y[:,NewAxis]), axis=1)
+            points = num.concatenate((x[:,num.NewAxis], y[:,num.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)
+            elevation = num.take(elevation, point_indices)
+            stage = num.take(stage, point_indices, axis=1)
 
             # Get info for location of maximal runup
-            points_in_polygon = take(points, point_indices)
+            points_in_polygon = num.take(points, point_indices)
             x = points_in_polygon[:,0]
             y = points_in_polygon[:,1]
         else:
             # Take all points
-            point_indices = arange(len(x))
+            point_indices = num.arange(len(x))
 
         # Temporal restriction
         time = fid.variables['time'][:]
-        all_timeindices = arange(len(time))
+        all_timeindices = num.arange(len(time))
         if time_interval is not None:
             msg = 'time_interval must be a sequence of length 2.'
@@ -7343 +7313 @@
 
             # Take time indices corresponding to interval (& is bitwise AND)
-            timesteps = compress((time_interval[0] <= time) \
+            timesteps = num.compress((time_interval[0] <= time) \
                                      & (time <= time_interval[1]),
-                                 all_timeindices)
+                                     all_timeindices)
 
             msg = 'time_interval %s did not include any model timesteps.' \
                   % time_interval
-            assert not alltrue(timesteps == 0), msg
+            assert not num.alltrue(timesteps == 0), msg
         else:
             # Take them all
@@ -7372 +7342 @@
             # Get wet nodes i.e. nodes with depth>0 within given region
             # and timesteps
-            wet_nodes = compress(depth > minimum_allowed_height,
-                                 arange(len(depth)))
-
-            if alltrue(wet_nodes == 0):
+            wet_nodes = num.compress(depth > minimum_allowed_height,
+                                     num.arange(len(depth)))
+
+            if num.alltrue(wet_nodes == 0):
                 runup = None
             else:
                 # Find maximum elevation among wet nodes
-                wet_elevation = take(elevation, wet_nodes)
-                runup_index = argmax(wet_elevation)
+                wet_elevation = num.take(elevation, wet_nodes)
+                runup_index = num.argmax(wet_elevation)
                 runup = max(wet_elevation)
                 assert wet_elevation[runup_index] == runup       # Must be True
@@ -7388 +7358 @@
 
                 # Record location
-                wet_x = take(x, wet_nodes)
-                wet_y = take(y, wet_nodes)
+                wet_x = num.take(x, wet_nodes)
+                wet_y = num.take(y, wet_nodes)
                 maximal_runup_location = [wet_x[runup_index],wet_y[runup_index]]