import struct, dbf, cPickle, time
import sqlite3, os.path, math

NULL_SHAPE = 0
POINT_SHAPE = 1
POLYLINE_SHAPE = 3
POLYGON_SHAPE = 5

def pnInPoly(pts, pt):
    c = False
    j = len(pts) - 1
    for i in xrange(len(pts)):
        if ((pts[i][1] <= pt[1]) and (pt[1] < pts[j][1])) or ((pts[j][1] <= pt[1]) and (pt[1] < pts[i][1])):
            if pt[0] < (float(pts[j][0] - pts[i][0]) * (pt[1] - pts[i][1]) / (pts[j][1] - pts[i][1]) + pts[i][0]):
                c = not c
        j = i
    return c 

def reader(filename, records = -1):
    f = open(filename, 'rb')
    f.seek(100)

    while 1:
        try:
            (number, length) = struct.unpack('>ii', f.read(8))
        except:
            print "end of file..."
            break

        record = f.read(length * 2)
        if ord(record[0]) == NULL_SHAPE:
            # Null shape
            assert (len(record) == 4)
            yield (number, 0, None)
        elif ord(record[0]) == POINT_SHAPE:
            # Point shape
            assert (len(record) == 20)
            (typ, x, y) = struct.unpack('<idd', record)
            yield (number, typ, (x, y))
        elif ord(record[0]) in (POLYLINE_SHAPE, POLYGON_SHAPE):
            # Polygon or polyline shape
            (typ, bbleft, bbtop, bbright, bbbottom, numParts, numPoints) = struct.unpack('<iddddii', record[:44])
            record = record[44:]
            parts = []
            for i in xrange(numParts):
                parts.append(struct.unpack('<i', record[i*4:(i+1)*4])[0])
            
            record = record[numParts * 4:]
            points = []
            for i in xrange(numPoints):
                points.append(struct.unpack('<dd', record[i * 16:(i + 1) * 16]))

            polygon = []
            for i in xrange(len(parts)):
                if i + 1 >= len(parts):
                    stop = -1
                else:
                    stop = parts[i + 1]
                current = parts[i]
                polygonpart = []
                while current != stop:
                    if current >= len(points):
                        break
                    polygonpart.append(points[current])
                    current += 1
                polygon.append(polygonpart)
            yield (number, typ, polygon)
        else:
            raise Exception('Unknown shape')

        records -= 1
        if records == 0:
            break
    f.close()
    
def isIn(index, pt):
    x = pt[0] / 1000000
    y = pt[1] / 1000000

    for poly in index.get((x, y), []):
        for p in poly[1]:
            if pnInPoly(p, pt):
                return poly[0]
    
    return None

def jtsk2wgs84(X, Y):
    # Prepocet vstupnich udaju
    H = 245

    # Vypocet zemepisnych souradnic z rovinnych souradnic
    a = 6377397.15508
    e = 0.081696831215303
    n = 0.97992470462083
    konst_u_ro = 12310230.12797036
    sinUQ = 0.863499969506341
    cosUQ = 0.504348889819882
    sinVQ = 0.420215144586493
    cosVQ = 0.907424504992097
    alfa = 1.000597498371542
    k = 1.003419163966575
    ro = math.sqrt(X * X + Y * Y)
    epsilon = 2 * math.atan(Y / (ro + X))
    D = epsilon / n
    S = 2 * math.atan(math.exp(1 / n * math.log(konst_u_ro / ro))) - math.pi / 2
    sinS = math.sin(S)
    cosS = math.cos(S)
    sinU = sinUQ * sinS - cosUQ * cosS * math.cos(D)
    cosU = math.sqrt(1 - sinU * sinU)
    sinDV = math.sin(D) * cosS / cosU
    cosDV = math.sqrt(1 - sinDV * sinDV)
    sinV = sinVQ * cosDV - cosVQ * sinDV
    cosV = cosVQ * cosDV + sinVQ * sinDV
    Ljtsk = 2 * math.atan(sinV / (1 + cosV)) / alfa
    t = math.exp(2 / alfa * math.log((1 + sinU) / cosU / k))
    pom = (t - 1) / (t + 1)
    
    while True:
        sinB = pom
        pom = t * math.exp(e * math.log((1 + e * sinB) / (1 - e * sinB)))
        pom = (pom - 1) / (pom + 1)
        if abs(pom - sinB) < 1e-15:
            break
    Bjtsk = math.atan(pom / math.sqrt(1 - pom * pom))

    # Pravouhle souradnice ve S-JTSK
    a = 6377397.15508
    f_1 = 299.152812853
    e2 = 1 - (1 - 1 / f_1) * (1 - 1/f_1)
    ro = a / math.sqrt(1 - e2 * math.sin(Bjtsk) * math.sin(Bjtsk))
    x = (ro + H) * math.cos(Bjtsk) * math.cos(Ljtsk)
    y = (ro + H) * math.cos(Bjtsk) * math.sin(Ljtsk)
    z = ((1 - e2) * ro + H) * math.sin(Bjtsk)

    # Pravouhle souradnice v WGS-84
    dx = 570.69
    dy = 85.69
    dz = 462.84
    wz = -5.2611 / 3600 * math.pi / 180
    wy = -1.58676 / 3600 * math.pi / 180
    wx = -4.99821 / 3600 * math.pi / 180
    m = 3.543e-6
    xn = dx + (1 + m) * (x + wz * y - wy * z)
    yn = dy + (1 + m) * (-wz * x + y + wx * z)
    zn = dz + (1 + m) * (wy * x - wx * y + z)

    # Geodeticke souradnice v systemu WGS-84
    a = 6378137.0
    f_1 = 298.257223563
    a_b = f_1 / (f_1 - 1)
    p = math.sqrt(xn * xn + yn * yn)
    e2 = 1 - (1 - 1 / f_1) * (1 - 1/f_1)
    theta = math.atan(zn * a_b / p)
    st = math.sin(theta)
    ct = math.cos(theta)
    t = (zn + e2 * a_b * a * st * st * st) / (p - e2 * a * ct * ct * ct)
    B = math.atan(t)
    L = 2 * math.atan(yn / (p + xn))
    H = math.sqrt(1 + t * t) * (p - a / math.sqrt(1 + (1 - e2) * t * t))

    # Format vystupnich hodnot
    return (180 * L / math.pi , 180 * B / math.pi) 

def createTables(cur):
    cur.execute("""CREATE TABLE IF NOT EXISTS Nodes(ID_NODE INTEGER PRIMARY KEY AUTOINCREMENT,
    X INTEGER NOT NULL, Y INTEGER NOT NULL)""")
    cur.execute("""CREATE TABLE IF NOT EXISTS Lines(ID_LINE INTEGER PRIMARY KEY AUTOINCREMENT,
    ID_DIBAVOD INTEGER)""")
    cur.execute("""CREATE TABLE IF NOT EXISTS Areas(ID_AREA INTEGER PRIMARY KEY AUTOINCREMENT,
    ID_DIBAVOD INTEGER)""")
    cur.execute("""CREATE TABLE IF NOT EXISTS Relations(ID_RELATION INTEGER PRIMARY KEY AUTOINCREMENT)""")
    cur.execute("""CREATE TABLE IF NOT EXISTS Dibavod(ID_DIBAVOD INTEGER PRIMARY KEY AUTOINCREMENT,
    NAME TEXT, OBJ_ID INTEGER NOT NULL)""")

    cur.execute("""CREATE TABLE IF NOT EXISTS LineNodes(ID_LINE INTEGER,
    ID_NODE INTEGER, SEQ INTEGER)""")
    cur.execute("""CREATE TABLE IF NOT EXISTS LineRelations(ID_RELATION INTEGER,
    ID_LINE INTEGER, SEQ INTEGER)""")

    cur.execute("""CREATE TABLE IF NOT EXISTS AreaNodes(ID_AREA INTEGER,
    ID_NODE INTEGER, SEQ INTEGER)""")
    cur.execute("""CREATE TABLE IF NOT EXISTS AreaRelations(ID_RELATION INTEGER,
    ID_AREA INTEGER, SEQ INTEGER)""")

def importVody(prefix, nameColumn, idColumn):
    global allNodes
    con = sqlite3.connect("vody.sqlite")
    con.isolation_level = None
    cur = con.cursor()

    createTables(cur)

    cur.execute("BEGIN")
    total = 0
    allLines = []
    lastTyp = None
    for (number, typ, obj) in reader(os.path.join('data', prefix + '.shp')):
        total += 1
        if typ in (POLYLINE_SHAPE, POLYGON_SHAPE):
            if typ == POLYLINE_SHAPE:
                tablePrefix = "Line"
            elif typ == POLYGON_SHAPE:
                tablePrefix = "Area"
                
            if lastTyp is None:
                lastTyp = typ
                
            assert lastTyp == typ
            
            lineIDs = []
            for line in obj:
                nodeIDs = []
                if typ == POLYGON_SHAPE:
                    line = line[:-1]
                for pt in line:
                    pt2 = jtsk2wgs84(-pt[1], -pt[0])
                    pt2 = (int(pt2[0] * 10000000 + 0.5), int(pt2[1] * 10000000 + 0.5))
                    if pt2 in allNodes:
                        nodeID = allNodes[pt2]
                    else:
                        cur.execute("insert into Nodes(X, Y) values (?, ?)", pt2)
                        nodeID = cur.lastrowid
                        allNodes[pt2] = nodeID
                    if len(nodeIDs) > 0 and nodeID == nodeIDs[-1]:
                        pass
                    else:
                        nodeIDs.append(nodeID)
                        
                cur.execute("insert into %ss(ID_DIBAVOD) values (NULL)" % tablePrefix)
                lineID = cur.lastrowid
                lineIDs.append(lineID)
                
                for i, nodeID in enumerate(nodeIDs):
                    cur.execute("insert into %sNodes(ID_%s, ID_NODE, SEQ) values (?, ?, ?)" % (tablePrefix, tablePrefix),
                        (lineID, nodeID, i + 1))
            allLines.append(lineIDs)

            if len(obj) > 1:
                print "Relations: ", len(obj)
                cur.execute("insert into Relations VALUES(NULL)")
                relationID = cur.lastrowid
                for i, lineID in enumerate(lineIDs):
                    cur.execute("insert into %sRelations(ID_RELATION, ID_%s, SEQ) values (?, ?, ?)" % (tablePrefix, tablePrefix),
                        (relationID, lineID, i + 1))
#                break                
        else:
            raise Exception('Unknown shape')
        #if total >= 10:
        #    break
        #break
        if total % 200 == 0:
            print total
    print total

    f = open(os.path.join('data', prefix + '.dbf'), "rb")
    db = dbf.dbfreader(f)

    fieldnames = db.next()
    fieldspecs = db.next()

    inazev = fieldnames.index(nameColumn)
    iid = fieldnames.index(idColumn)

    i = 0
    for row in db:
        n = unicode(row[inazev], 'cp1250').strip()
        did = long(row[iid])
        if len(n) == 0:
            n = None
        cur.execute("INSERT INTO Dibavod(NAME, OBJ_ID) VALUES (?, ?)",
            (n, did))
        idDibavod = cur.lastrowid
        for line in allLines[i]:
            cur.execute("UPDATE %ss SET ID_DIBAVOD=? WHERE ID_%s=?" % (tablePrefix, tablePrefix),
                (idDibavod, line))
        i += 1
        if i % 200 == 0:
            print i
        if i == len(allLines):
            break
    f.close()
    cur.execute("COMMIT")
    con.close()
       
if __name__ == '__main__':
    allNodes = {}
    importVody('A02_Vodni_tok_JU', 'NAZ_TOK', 'UTOKJ_ID')
    #importVody('A05_Vodni_nadrze', 'NAZ_NA', 'NADR_GID')    
    #70402
    #71922 - 2341 + 821
    