###############################################
#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 anuga.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()




    def NOT_test_dem2pts(self):
        """Test conversion from dem in ascii format to native NetCDF xya format
        """

        import time, os
        from Numeric import array, zeros, allclose, Float, concatenate
        from Scientific.IO.NetCDF import NetCDFFile

        #Write test asc file
        root = 'demtest'

        filename = root+'.asc'
        fid = open(filename, 'w')
        fid.write("""ncols         5
nrows         6
xllcorner     2000.5
yllcorner     3000.5
cellsize      25
NODATA_value  -9999
""")
        #Create linear function

        ref_points = []
        ref_elevation = []
        for i in range(6):
            y = (6-i)*25.0
            for j in range(5):
                x = j*25.0
                z = x+2*y

                ref_points.append( [x,y] )
                ref_elevation.append(z)
                fid.write('%f ' %z)
            fid.write('\n')

        fid.close()

        #Write prj file with metadata
        metafilename = root+'.prj'
        fid = open(metafilename, 'w')


        fid.write("""Projection UTM
Zone 56
Datum WGS84
Zunits NO
Units METERS
Spheroid WGS84
Xshift 0.0000000000
Yshift 10000000.0000000000
Parameters
""")
        fid.close()

        #Convert to NetCDF pts
        convert_dem_from_ascii2netcdf(root)
        dem2pts(root)

        #Check contents
        #Get NetCDF
        fid = NetCDFFile(root+'.pts', 'r')

        # Get the variables
        #print fid.variables.keys()
        points = fid.variables['points']
        elevation = fid.variables['elevation']

        #Check values

        #print points[:]
        #print ref_points
        assert allclose(points, ref_points)

        #print attributes[:]
        #print ref_elevation
        assert allclose(elevation, ref_elevation)

        #Cleanup
        fid.close()


        os.remove(root + '.pts')
        os.remove(root + '.dem')
        os.remove(root + '.asc')
        os.remove(root + '.prj')



    def NOT_test_dem2pts_bounding_box(self):
        """Test conversion from dem in ascii format to native NetCDF xya format
        """

        import time, os
        from Numeric import array, zeros, allclose, Float, concatenate
        from Scientific.IO.NetCDF import NetCDFFile

        #Write test asc file
        root = 'demtest'

        filename = root+'.asc'
        fid = open(filename, 'w')
        fid.write("""ncols         5
nrows         6
xllcorner     2000.5
yllcorner     3000.5
cellsize      25
NODATA_value  -9999
""")
        #Create linear function

        ref_points = []
        ref_elevation = []
        for i in range(6):
            y = (6-i)*25.0
            for j in range(5):
                x = j*25.0
                z = x+2*y

                ref_points.append( [x,y] )
                ref_elevation.append(z)
                fid.write('%f ' %z)
            fid.write('\n')

        fid.close()

        #Write prj file with metadata
        metafilename = root+'.prj'
        fid = open(metafilename, 'w')


        fid.write("""Projection UTM
Zone 56
Datum WGS84
Zunits NO
Units METERS
Spheroid WGS84
Xshift 0.0000000000
Yshift 10000000.0000000000
Parameters
""")
        fid.close()

        #Convert to NetCDF pts
        convert_dem_from_ascii2netcdf(root)
        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
                northing_min=3035.0, northing_max=3125.5)

        #Check contents
        #Get NetCDF
        fid = NetCDFFile(root+'.pts', 'r')

        # Get the variables
        #print fid.variables.keys()
        points = fid.variables['points']
        elevation = fid.variables['elevation']

        #Check values
        assert fid.xllcorner[0] == 2010.0
        assert fid.yllcorner[0] == 3035.0

        #create new reference points
        ref_points = []
        ref_elevation = []
        for i in range(4):
            y = (4-i)*25.0 + 25.0
            y_new = y + 3000.5 - 3035.0
            for j in range(4):
                x = j*25.0 + 25.0
                x_new = x + 2000.5 - 2010.0
                z = x+2*y

                ref_points.append( [x_new,y_new] )
                ref_elevation.append(z)

        #print points[:]
        #print ref_points
        assert allclose(points, ref_points)

        #print attributes[:]
        #print ref_elevation
        assert allclose(elevation, ref_elevation)

        #Cleanup
        fid.close()


        os.remove(root + '.pts')
        os.remove(root + '.dem')
        os.remove(root + '.asc')
        os.remove(root + '.prj')



    def NOT_test_dem2pts_remove_Nullvalues(self):
        """Test conversion from dem in ascii format to native NetCDF xya format
        """

        import time, os
        from Numeric import array, zeros, allclose, Float, concatenate
        from Scientific.IO.NetCDF import NetCDFFile

        #Write test asc file
        root = 'demtest'

        filename = root+'.asc'
        fid = open(filename, 'w')
        fid.write("""ncols         5
nrows         6
xllcorner     2000.5
yllcorner     3000.5
cellsize      25
NODATA_value  -9999
""")
        #Create linear function
        # ref_ will write all the values
        # new_ref_ will write the values except for NODATA_values
        ref_points = []
        ref_elevation = []
        new_ref_pts = []
        new_ref_elev = []
        NODATA_value = -9999
        for i in range(6):
            y = (6-i)*25.0
            for j in range(5):
                x = j*25.0
                z = x+2*y
                if j == 4: z = NODATA_value # column
                if i == 2 and j == 2: z = NODATA_value # random
                if i == 5 and j == 1: z = NODATA_value
                if i == 1: z = NODATA_value # row
                if i == 3 and j == 1: z = NODATA_value # two pts/row
                if i == 3 and j == 3: z = NODATA_value


                if z <> NODATA_value:
                    new_ref_elev.append(z)
                    new_ref_pts.append( [x,y] )

                ref_points.append( [x,y] )
                ref_elevation.append(z)

                fid.write('%f ' %z)
            fid.write('\n')

        fid.close()


        #Write prj file with metadata
        metafilename = root+'.prj'
        fid = open(metafilename, 'w')


        fid.write("""Projection UTM
Zone 56
Datum WGS84
Zunits NO
Units METERS
Spheroid WGS84
Xshift 0.0000000000
Yshift 10000000.0000000000
Parameters
""")
        fid.close()

        #Convert to NetCDF pts
        convert_dem_from_ascii2netcdf(root)
        dem2pts(root)

        #Check contents
        #Get NetCDF
        fid = NetCDFFile(root+'.pts', 'r')

        # Get the variables
        #print fid.variables.keys()
        points = fid.variables['points']
        elevation = fid.variables['elevation']

        #Check values
        #print 'points', points[:]
        assert len(points) == len(new_ref_pts), 'length of returned points not correct'
        assert allclose(points, new_ref_pts), 'points do not align'

        #print 'elevation', elevation[:]
        assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct'
        assert allclose(elevation, new_ref_elev), 'elevations do not align'

        #Cleanup
        fid.close()


        os.remove(root + '.pts')
        os.remove(root + '.dem')
        os.remove(root + '.asc')
        os.remove(root + '.prj')

    def NOT_test_dem2pts_bounding_box_Nullvalues(self):
        """Test conversion from dem in ascii format to native NetCDF xya format
        """

        import time, os
        from Numeric import array, zeros, allclose, Float, concatenate
        from Scientific.IO.NetCDF import NetCDFFile

        #Write test asc file
        root = 'demtest'

        filename = root+'.asc'
        fid = open(filename, 'w')
        fid.write("""ncols         5
nrows         6
xllcorner     2000.5
yllcorner     3000.5
cellsize      25
NODATA_value  -9999
""")
        #Create linear function

        ref_points = []
        ref_elevation = []
        new_ref_pts1 = []
        new_ref_elev1 = []
        NODATA_value = -9999
        for i in range(6):
            y = (6-i)*25.0
            for j in range(5):
                x = j*25.0
                z = x+2*y
                if j == 4: z = NODATA_value # column
                if i == 2 and j == 2: z = NODATA_value # random
                if i == 5 and j == 1: z = NODATA_value
                if i == 1: z = NODATA_value # row
                if i == 3 and j == 1: z = NODATA_value # two pts/row
                if i == 3 and j == 3: z = NODATA_value

                if z <> NODATA_value:
                    new_ref_elev1.append(z)
                    new_ref_pts1.append( [x,y] )

                ref_points.append( [x,y] )
                ref_elevation.append(z)
                fid.write('%f ' %z)
            fid.write('\n')

        fid.close()

        #Write prj file with metadata
        metafilename = root+'.prj'
        fid = open(metafilename, 'w')


        fid.write("""Projection UTM
Zone 56
Datum WGS84
Zunits NO
Units METERS
Spheroid WGS84
Xshift 0.0000000000
Yshift 10000000.0000000000
Parameters
""")
        fid.close()

        #Convert to NetCDF pts
        convert_dem_from_ascii2netcdf(root)
        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
                northing_min=3035.0, northing_max=3125.5)

        #Check contents
        #Get NetCDF
        fid = NetCDFFile(root+'.pts', 'r')

        # Get the variables
        #print fid.variables.keys()
        points = fid.variables['points']
        elevation = fid.variables['elevation']

        #Check values
        assert fid.xllcorner[0] == 2010.0
        assert fid.yllcorner[0] == 3035.0

        #create new reference points
        ref_points = []
        ref_elevation = []
        new_ref_pts2 = []
        new_ref_elev2 = []
        for i in range(4):
            y = (4-i)*25.0 + 25.0
            y_new = y + 3000.5 - 3035.0
            for j in range(4):
                x = j*25.0 + 25.0
                x_new = x + 2000.5 - 2010.0
                z = x+2*y

                if j == 3: z = NODATA_value # column
                if i == 1 and j == 1: z = NODATA_value # random
                if i == 4 and j == 0: z = NODATA_value
                if i == 0: z = NODATA_value # row
                if i == 2 and j == 0: z = NODATA_value # two pts/row
                if i == 2 and j == 2: z = NODATA_value

                if z <> NODATA_value:
                    new_ref_elev2.append(z)
                    new_ref_pts2.append( [x_new,y_new] )


                ref_points.append( [x_new,y_new] )
                ref_elevation.append(z)

        #print points[:]
        #print ref_points
        #assert allclose(points, ref_points)

        #print attributes[:]
        #print ref_elevation
        #assert allclose(elevation, ref_elevation)


        assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
        assert allclose(points, new_ref_pts2), 'points do not align'

        #print 'elevation', elevation[:]
        assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct'
        assert allclose(elevation, new_ref_elev2), 'elevations do not align'
        #Cleanup
        fid.close()


        os.remove(root + '.pts')
        os.remove(root + '.dem')
        os.remove(root + '.asc')
        os.remove(root + '.prj')

#********************
#*** END OF OBSOLETE FUNCTIONS
#***************