source: inundation/pyvolution/data_manager_obsolete_stuff.py @ 2501

###############################################
#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 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
#***************