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 _______________________________________________ Talk-cz mailing list Talk-cz@openstreetmap.org http://lists.openstreetmap.org/listinfo/talk-cz