[Talk-cz] import dibavod C02 -- koupaci oblasti
Pavel Machek
pavel na ucw.cz
Sobota Říjen 23 20:50:07 UTC 2010
Ahoj!
Uvazuju jak importovat "koupaci oblasti".
Tenhle import je generovan ponekud jednodussim hackem, ale jsou to
jenom body -> melo by to byt ok.
Aktualni verse je:
<?xml version='1.0' encoding='UTF-8'?>
<osm version='0.5' generator='JOSM'>
<node id='-264' action='modify' timestamp='2010-10-23T20:41:21Z'
visible='true' lat='49.650985' lon='16.604522'>
<tag k='note' v='koupaci oblast ve volne prirode dle DIBAVODu' />
<tag k='source' v='dibavod C02' />
<tag k='sport' v='swimming' />
</node>
Je nejake vhodnejsi tagovani?
Hmm.... v databazi to vypada takhle...
Rekr id : KO530901
Kraj : PARDUBICKÝ
Koup naz : písník Březhrad (u nádra<9E>í)
Tok naz : Plačický potok
Voda typ : P
Hlgp vuv : 103010170
Tok id : 104550000100
Pozn : není přímé spojení s vodním tokem
Idvt : 10100651
Nadr gid : 103010170003
Orp : Pardubice
Orp id : 5309
Obec : Opatovice nad Labem
Obec id : 575429
...a uvazuju ze ten nazev koupaliste by se asi hodil... Zatim to
vyrabim *strasnym* hackem (prilozen), ktery byl puvodne urcen na
plneni databaze, ne vyrobu xmlka... Bohuzel se v nem moc nevyznam :-(.
Napady?
Pavel
#!/usr/bin/python
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 "</osm>"
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
elif typ == POINT_SHAPE:
( x, y ) = obj
( lon, lat ) = jtsk2wgs84(-y, -x)
print "<node id='%d' lon='%f' lat='%f' />" % ( -total, lon, lat )
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()
if nameColumn != '':
inazev = fieldnames.index(nameColumn)
else:
inazev = None
iid = fieldnames.index(idColumn)
i = 0
for row in db:
n = None
if inazev:
n = unicode(row[inazev], 'cp1250').strip()
if len(n) == 0:
n = None
did = long(row[iid])
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 = {}
if 1:
print "<?xml version='1.0' encoding='UTF-8'?>"
print "<osm version='0.5' generator='shapefile'>"
#importVody('A02_Vodni_tok_JU', 'NAZ_TOK', 'UTOKJ_ID')
#importVody('A05_Vodni_nadrze', 'NAZ_NA', 'NADR_GID')
#importVody('A01_Vodni_tok_CEVT', 'NAZ_TOK', 'TOK_ID')
#importVody('A04zvm_Melioracni_kanaly', '', 'LCRM_ID')
#importVody('A06_Bazina_mocal', '', 'ET_ID')
#importVody('I01zvm_Jezy', 'TOK_ID', 'OT06_ID')
importVody('C02_Koupaci_oblasti', '', '')
--
(english) http://www.livejournal.com/~pavelmachek
(cesky, pictures) http://atrey.karlin.mff.cuni.cz/~pavel/picture/horses/blog.html
Další informace o konferenci talk-cz