Changeset 2344

Timestamp:

Feb 7, 2006, 1:20:39 PM (19 years ago)

Author:

ole

Message:

Moved obsolete stuff out from data_manager

Location:

inundation/pyvolution

Files:

• data_manager.py (modified) (2 diffs)
• data_manager_obsolete_stuff.py (added)

inundation/pyvolution/data_manager.py

@@ -2834 +2834 @@
 
 
-
-
-
-
-
-
-
-
-###############################################
-#OBSOLETE STUFF
-#Native checkpoint format.
-#Information needed to recreate a state is preserved
-#FIXME: Rethink and maybe use netcdf format
-def cpt_variable_writer(filename, t, v0, v1, v2):
-    """Store all conserved quantities to file
-    """
-
-    M, N = v0.shape
-
-    FN = create_filename(filename, 'cpt', M, t)
-    #print 'Writing to %s' %FN
-
-    fid = open(FN, 'w')
-    for i in range(M):
-        for j in range(N):
-            fid.write('%.16e ' %v0[i,j])
-        for j in range(N):
-            fid.write('%.16e ' %v1[i,j])
-        for j in range(N):
-            fid.write('%.16e ' %v2[i,j])
-
-        fid.write('\n')
-    fid.close()
-
-
-def cpt_variable_reader(filename, t, v0, v1, v2):
-    """Store all conserved quantities to file
-    """
-
-    M, N = v0.shape
-
-    FN = create_filename(filename, 'cpt', M, t)
-    #print 'Reading from %s' %FN
-
-    fid = open(FN)
-
-
-    for i in range(M):
-        values = fid.readline().split() #Get one line
-
-        for j in range(N):
-            v0[i,j] = float(values[j])
-            v1[i,j] = float(values[3+j])
-            v2[i,j] = float(values[6+j])
-
-    fid.close()
-
-def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2):
-    """Writes x,y,z,z,z coordinates of triangles constituting the bed
-    elevation.
-    FIXME: Not in use pt
-    """
-
-    M, N = v0.shape
-
-
-    print X0
-    import sys; sys.exit()
-    FN = create_filename(filename, 'cpt', M)
-    print 'Writing to %s' %FN
-
-    fid = open(FN, 'w')
-    for i in range(M):
-        for j in range(2):
-            fid.write('%.16e ' %X0[i,j])   #x, y
-        for j in range(N):
-            fid.write('%.16e ' %v0[i,j])       #z,z,z,
-
-        for j in range(2):
-            fid.write('%.16e ' %X1[i,j])   #x, y
-        for j in range(N):
-            fid.write('%.16e ' %v1[i,j])
-
-        for j in range(2):
-            fid.write('%.16e ' %X2[i,j])   #x, y
-        for j in range(N):
-            fid.write('%.16e ' %v2[i,j])
-
-        fid.write('\n')
-    fid.close()
-
-
-
-#Function for storing out to e.g. visualisation
-#FIXME: Do we want this?
-#FIXME: Not done yet for this version
-def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2):
-    """Writes x,y,z coordinates of triangles constituting the bed elevation.
-    """
-
-    M, N = v0.shape
-
-    FN = create_filename(filename, 'dat', M)
-    #print 'Writing to %s' %FN
-
-    fid = open(FN, 'w')
-    for i in range(M):
-        for j in range(2):
-            fid.write('%f ' %X0[i,j])   #x, y
-        fid.write('%f ' %v0[i,0])       #z
-
-        for j in range(2):
-            fid.write('%f ' %X1[i,j])   #x, y
-        fid.write('%f ' %v1[i,0])       #z
-
-        for j in range(2):
-            fid.write('%f ' %X2[i,j])   #x, y
-        fid.write('%f ' %v2[i,0])       #z
-
-        fid.write('\n')
-    fid.close()
-
-
-
-def dat_variable_writer(filename, t, v0, v1, v2):
-    """Store water height to file
-    """
-
-    M, N = v0.shape
-
-    FN = create_filename(filename, 'dat', M, t)
-    #print 'Writing to %s' %FN
-
-    fid = open(FN, 'w')
-    for i in range(M):
-        fid.write('%.4f ' %v0[i,0])
-        fid.write('%.4f ' %v1[i,0])
-        fid.write('%.4f ' %v2[i,0])
-
-        fid.write('\n')
-    fid.close()
-
-
-def read_sww(filename):
-    """Read sww Net CDF file containing Shallow Water Wave simulation
-
-    The integer array volumes is of shape Nx3 where N is the number of
-    triangles in the mesh.
-
-    Each entry in volumes is an index into the x,y arrays (the location).
-
-    Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions
-    number_of_timesteps, number_of_points.
-
-    The momentum is not always stored.
-
-    """
-    from Scientific.IO.NetCDF import NetCDFFile
-    print 'Reading from ', filename
-    fid = NetCDFFile(filename, 'r')    #Open existing file for read
-#latitude, longitude
-    # Get the variables as Numeric arrays
-    x = fid.variables['x']             #x-coordinates of vertices
-    y = fid.variables['y']             #y-coordinates of vertices
-    z = fid.variables['elevation']     #Elevation
-    time = fid.variables['time']       #Timesteps
-    stage = fid.variables['stage']     #Water level
-    #xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
-    #ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
-
-    volumes = fid.variables['volumes'] #Connectivity
-
-    #FIXME (Ole): What is this?
-    #             Why isn't anything returned?
-    #             Where's the unit test?
-
 def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None,
                  verbose=False):
@@ -3138 +2962 @@
 
 
-
-def sww2asc_obsolete(basename_in, basename_out = None,
-            quantity = None,
-            timestep = None,
-            reduction = None,
-            cellsize = 10,
-            verbose = False,
-            origin = None):
-    """Read SWW file and convert to Digitial Elevation model format (.asc)
-
-    Example:
-
-    ncols         3121
-    nrows         1800
-    xllcorner     722000
-    yllcorner     5893000
-    cellsize      25
-    NODATA_value  -9999
-    138.3698 137.4194 136.5062 135.5558 ..........
-
-    Also write accompanying file with same basename_in but extension .prj
-    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
-
-
-    if quantity is given, out values from quantity otherwise default to
-    elevation
-
-    if timestep (an index) is given, output quantity at that timestep
-
-    if reduction is given use that to reduce quantity over all timesteps.
-
-    """
-    from Numeric import array, Float, concatenate, NewAxis, zeros,\
-         sometrue
-    from utilities.polygon import inside_polygon
-
-    #FIXME: Should be variable
-    datum = 'WGS84'
-    false_easting = 500000
-    false_northing = 10000000
-
-    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'
-    ascfile = basename_out + '.asc'
-    prjfile = basename_out + '.prj'
-
-
-    if verbose: print 'Reading from %s' %swwfile
-    #Read sww file
-    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'][:]
-
-    ymin = min(y); ymax = max(y)
-    xmin = min(x); xmax = max(x)
-
-    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]
-
-
-    #Get quantity and reduce if applicable
-    if verbose: print 'Reading quantity %s' %quantity
-
-    if quantity.lower() == 'depth':
-        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
-    else:
-        q = fid.variables[quantity][:]
-
-
-    if len(q.shape) == 2:
-        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
-
-    #Now q has dimension: number_of_points
-
-    #Create grid and update xll/yll corner and x,y
-    if verbose: print 'Creating grid'
-    ncols = int((xmax-xmin)/cellsize)+1
-    nrows = int((ymax-ymin)/cellsize)+1
-
-    newxllcorner = xmin+xllcorner
-    newyllcorner = ymin+yllcorner
-
-    x = x+xllcorner-newxllcorner
-    y = y+yllcorner-newyllcorner
-
-    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
-    assert len(vertex_points.shape) == 2
-
-
-    from Numeric import zeros, Float
-    grid_points = zeros ( (ncols*nrows, 2), Float )
-
-
-    for i in xrange(nrows):
-        yg = i*cellsize
-        for j in xrange(ncols):
-            xg = j*cellsize
-            k = i*ncols + j
-
-            grid_points[k,0] = xg
-            grid_points[k,1] = yg
-
-    #Interpolate
-    from least_squares import Interpolation
-
-
-    #FIXME: This should be done with precrop = True, otherwise it'll
-    #take forever. With expand_search  set to False, some grid points might
-    #miss out....
-
-    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
-                           precrop = False, expand_search = True,
-                           verbose = verbose)
-
-    #Interpolate using quantity values
-    if verbose: print 'Interpolating'
-    grid_values = interp.interpolate(q).flat
-
-    #Write
-    #Write prj file
-    if verbose: print 'Writing %s' %prjfile
-    prjid = open(prjfile, 'w')
-    prjid.write('Projection    %s\n' %'UTM')
-    prjid.write('Zone          %d\n' %zone)
-    prjid.write('Datum         %s\n' %datum)
-    prjid.write('Zunits        NO\n')
-    prjid.write('Units         METERS\n')
-    prjid.write('Spheroid      %s\n' %datum)
-    prjid.write('Xshift        %d\n' %false_easting)
-    prjid.write('Yshift        %d\n' %false_northing)
-    prjid.write('Parameters\n')
-    prjid.close()
-
-
-
-    if verbose: print 'Writing %s' %ascfile
-    NODATA_value = -9999
-
-    ascid = open(ascfile, 'w')
-
-    ascid.write('ncols         %d\n' %ncols)
-    ascid.write('nrows         %d\n' %nrows)
-    ascid.write('xllcorner     %d\n' %newxllcorner)
-    ascid.write('yllcorner     %d\n' %newyllcorner)
-    ascid.write('cellsize      %f\n' %cellsize)
-    ascid.write('NODATA_value  %d\n' %NODATA_value)
-
-
-    #Get bounding polygon from mesh
-    P = interp.mesh.get_boundary_polygon()
-    inside_indices = inside_polygon(grid_points, P)
-
-    for i in range(nrows):
-        if verbose and i%((nrows+10)/10)==0:
-            print 'Doing row %d of %d' %(i, nrows)
-
-        for j in range(ncols):
-            index = (nrows-i-1)*ncols+j
-
-            if sometrue(inside_indices == index):
-                ascid.write('%f ' %grid_values[index])
-            else:
-                ascid.write('%d ' %NODATA_value)
-
-        ascid.write('\n')
-
-    #Close
-    ascid.close()
-    fid.close()
-
-#********************
-#*** END OF OBSOLETE FUNCTIONS
-#***************