Datafile.py

'''
Created on Feb 2, 2011

@author: john
'''
import sqlite3
import scipy
import random
#Input: (time, lat, lon, a, b, c,)

def getValue(timestep, p, a, b, c):
    """getValue(args) -> tuple(direction, velocity, water surface elevation)
    
    Given a timestep in seconds, a (latitude, longitude) tuple of a point inside a triangle,
    and three numbers describing the points of that triangle, returns the interpolated
    values for the three values in the database"""
    conn = sqlite3.connect('dahv.sqlite')
    cursor = conn.cursor()
    
    #Pull a, b, and c points from the database
    # fetchall()[0] returns just the first item (i.e. only item)
    # fetchall()[0][1:3] returns just items 1 and 2 (i.e. lat and lon) from that first item
    # scipy.array() just wraps the whole thing in an array for matrix calculations
    res = cursor.execute('''select * from coordinates where num=%s''' % a).fetchall()[0][1:3]
    A = scipy.array(res)
    res = cursor.execute('''select * from coordinates where num=%s''' % b).fetchall()[0][1:3]
    B = scipy.array(res)
    res = cursor.execute('''select * from coordinates where num=%s''' % c).fetchall()[0][1:3]
    C = scipy.array(res)

    P = scipy.array(p)
    
    # Pull the values for a, b, and c.
    # Note: there should only be one valueset per number per timestep, hence fetchall()[0]
    # We only grab the last 3 values, because we're only interested in the direction, velocity, and water surface elevation
    a = cursor.execute('''select * from data where num=%s and time=%s''' % (a, timestep)).fetchall()[0][3:]
    b = cursor.execute('''select * from data where num=%s and time=%s''' % (b, timestep)).fetchall()[0][3:]
    c = cursor.execute('''select * from data where num=%s and time=%s''' % (c, timestep)).fetchall()[0][3:]

    #Interpolate for all three values
    vals = []
    for i in xrange(3):
        vals.append(interpN(P, A, B, C, a[i], b[i], c[i]))
    return(vals)
    
def interpN(P, A, B, C, a, b, c):
    """interpSurface(*args) -> interpolated value
    
    Returns an n-dimensional interpolation for a value on a surface
    
    Expects 4 tuples and 3 values:
        P -> Positional scipy.array describing the point in question
        A,B,C -> Positional scipy.arrays describing points of surrounding triangle ABC
        a,b,c -> values of interest for points ABC which, interpolated, describe point P's value
    """
    # Culled from http://www.blackpawn.com/texts/pointinpoly/default.html
    
    # Compute vectors        
    v0 = C - A
    v1 = B - A
    v2 = P - A
    
    # Compute dot products
    dot = lambda x,y : scipy.dot(x,y)
    dot00 = dot(v0, v0)
    dot01 = dot(v0, v1)
    dot02 = dot(v0, v2)
    dot11 = dot(v1, v1)
    dot12 = dot(v1, v2)
    
    # Compute barycentric coordinates
    invDenom = 1 / (dot00 * dot11 - dot01 * dot01)
    u = (dot11 * dot02 - dot01 * dot12) * invDenom
    v = (dot00 * dot12 - dot01 * dot02) * invDenom
    if (u < 0) or (v < 0) or (u+v > 1):
        raise ValueError('u and v value place point p outside triangle abc')

    value = a + v*(b - a) + u * (c - a)
    return value

if __name__ == "__main__":
    # Test data:
    a,b,c =  (5337, 5528, 5527)
    p = (-123.93706, 45.8849480)

    #print the resulting interpolation
    print(getValue(4800, p, a, b, c))

sql_input.py

'''
Created on Feb 1, 2011

@author: john
'''
import sqlite3
import Datapoint
from math import sqrt


###########################################################
### This is only used to populate the database initially
def insert():
    """Use the information provided to populate a sqlite database with
    the data found in an object that is returned by PySELFE.Dataset.read_time_series()"""
    
    # Hardcoded datapaths on my system
    datapath = "/Users/john/Dropbox/Cannon Beach/data/"
    #elev_pickle = "elev.timeseries.out"
    dahv_pickle = "dahv.timeseries.out"
    
    # Water surface elevation is actually in dahv, so we don't need them both
    #elev = cPickle.load(open(datapath + elev_pickle,'r'))
    
    #Unpickle the dahv.timeseries pickle
    import cPickle
    dahv = cPickle.load(open(datapath + dahv_pickle,'r'))
    del cPickle #Big file, save some space
    
    #Open a DB connection to the given file and create a cursor
    conn = sqlite3.connect('dahv.sqlite')
    c = conn.cursor()
    
    # Create tables
    
    # time = time in seconds
    # iter = iteration number
    # num    = number of datapoint in the list. 
    #        Note: pulling data for all times for a given x and ordering by
    #        time will result in a timeseries of data for that point.
    # direction = What I believe is direction
    # velocity = What I believe is velocity
    # elevation = Water surface elevation
    c.execute('''drop table data''')
    c.execute('''drop table coordinates''')
    c.execute('''create table data
                (time real, iter real, num real, direction real, velocity real, elevation real)''')
    c.execute('''create table coordinates 
                (num real, lon real, lat real, bath real)''')

    # Cycle through the data and fill the database.
    # The "i" value here is the element's number in the lists described
    # in the note above.
    for i in xrange(len(dahv[0])):
        t = dahv[0][i]          # Zeroth element at each timestep is time in seconds
        t_iter = dahv[1][i]     # First element at each timestep is iteration number
        eta = dahv[2][i]        # Second element at each timestep is water surface elevation
        data = dahv[4][i]       # fourth element at each timestep is the data value tuple
    
        # Data list for this timestep at this position
        for j in xrange(len(data)):
            val = data[j] #We're assuming these equate!
            c.execute("""insert into data
                   values(%i, %i, %i, %f, %f, %f)""" % (t, t_iter, j, val[0][0], val[0][1], eta[j]))
    # Cycle through the data in the hgrid.LL file and pull the lat, lon, and bathymetry
    lines = [j.split() for j in [i[:-1] for i in open(datapath+"hgrid.LL", 'r').readlines()[2:742450]]]

    n = 1
    for row in lines:
        print("line %i" % n)
        c.execute('''insert into coordinates values(%i, %f, %f, %f)''' % (int(row[0]), float(row[1]), float(row[2]), float(row[3])))
        n += 1

    #Index num
    c.execute('''create index number on data (num)''')
    # Save (commit) the changes
    conn.commit()
    
    # We can also close the cursor if we are done with it
    c.close()

if __name__ == "__main__":
    insert()