Ahoj, cvičně jsem si taky zkusil tři KÚ, první dojmy:
* Skript by měl slučovat duplicitní uzly (triviální). * Skript by měl eliminovat uzel v cestě, pokud je totožný jako předcházející uzel (triviální). * Musí se ručně řešit multipolygony. * Spousta sousedících polygonů má překryvy, i v rámci jednoho KU. * LPIS obsahuje jen zemědělskou půdu, takže v mapě vznikají "ostrůvky" polygonů místo souvislé plochy. * Podezřele hodně luk je označeno jako oplocených (kul. 71) -- skript předpokládá, že každá pastvina má ohradník/plot? * Jak už zaznělo, budou potíže na rozhraních s existujícími lesy. * LPIS na řadě míst nekopíruje hranice parcel, takže budou nesoulady s daty RUIANu, plochami kreslenými podle KM atd. * Od pohledu mi přijde, že parcely RUIANu by byly kvalitnější zdroj pro plošné mapování. Na druhou stranu, už teď skript ušetří mraky ručního obkreslování, stačí jen v JOSM doopravit problémy. :-) (V příloze je skript s přidaným primitivním odstraňováním duplicitních uzlů. Nicméně, Python vůbec neznám a z jeho syntaxe dostávám osypky, takže zbytek nechám pythonistům :-)) Martin On 22.7.2014 15:28, Pavel Machek wrote:
Hi! I'd like to start import of LPIS farmland database, as we have very good coverage of houses, forests and water, but farmland is very good at places and completely missing at different places. Import page is at https://wiki.openstreetmap.org/wiki/LPIS , import script is attached to this email. Best regards, Pavel _______________________________________________ Talk-cz mailing list [email protected] https://lists.openstreetmap.org/listinfo/talk-cz
#!/usr/bin/python # http://www.lpis.cz/index.html # Mapa: # http://eagri.cz/public/app/lpisext/lpis/verejny/ # Vyzadat data na: # http://eagri.cz/public/app/lpisext/lpis/verejny/exportDat.html # Vyzadat: Pudni bloky # Souradny systrem: wgs-84 # Aha, a ty cisla katastralnich uzemi jde dokonce vykoukat z osm, jako ref u boundary=administrative, admin_level=10 # A... nefunguje jim dobre captcha, takze jde stahnout vic uzemi jednim formularem. # doubek: 631035 # babice: 600601 # for A in *.zip; do unzip $A; done import sys import shapefile from pyproj import * print "<?xml version='1.0' encoding='UTF-8'?>" print "<osm version='0.6' generator='pyshp'>" #p1 = Proj(init='esri:102067') p1 = Proj(init='EPSG:4326') p2 = Proj(init='EPSG:4326') sf = shapefile.Reader(sys.argv[1]) global field field = {} def a_field(num, name): global field i1, i2, i3, i4 = sf.fields[num] if name != i1: print 'Unexpected field name <', name field[name] = num a_field(1, 'ID_FB') a_field(3, 'NKODFB') a_field(6, 'PLATNYOD') a_field(7, 'VYMERAM') a_field(8, 'KULTURA') a_field(9, 'KULTURA_KL') a_field(10, 'KULTURAOD') a_field(11, 'EKO') a_field(21, 'VYSKA') a_field(22, 'SVAZITOST') a_field(26, 'KU_KOD') global nodeid nodeid = 0 global nodemap nodemap = {} def convert(point): lon, lat = point lon, lat = transform(p1, p2, lon, lat) #lat += 50.013082018919185-50.013834 #lon += 14.748617781670555-14.749757 return lon, lat def write_point(point): global nodeid global nodemap lon, lat = convert(point) nodekey = '%f/%f' % (lon, lat) if nodekey in nodemap: return nodemap[nodekey] tags = '<tag k="created_by" v="shpupload"/>' tags += '<tag k="source" v="lpis"/>' nodeid -= 1 print '<node id="%d" lon="%f" lat="%f\">%s</node>' % ( nodeid, lon, lat, tags ) nodemap[nodekey] = nodeid; return nodeid def attr(shrec, name): return shrec.record[field[name]] for shrec in sf.shapeRecords(): shape = shrec.shape prevpt = 6666 pts = [] for point in shape.points: curpt = write_point(point) if prevpt != curpt: pts += [ curpt ] prevpt = curpt nodeid -= 1 print '<way id="%d">' % nodeid print ' <tag k="created_by" v="pyshp"/>' print ' <tag k="id_fb" v="%s"/>' % attr(shrec, 'ID_FB') print ' <tag k="source" v="lpis"/>' print ' <tag k="lpis:kultura" v="%d"/>' % attr(shrec, 'KULTURA') kul = attr(shrec, 'KULTURA') if kul == 2: print ' <tag k="landuse" v="farmland"/>' elif kul == 3: print ' <tag k="landuse" v="farmland"/><tag k="crop" v="hop"/>' elif kul == 30: print ' <tag k="landuse" v="farmland"/><tag k="crop" v="hop"/>' elif kul == 31: print ' <tag k="landuse" v="farmland"/><tag k="crop" v="hop"/>' elif kul == 41: print ' <tag k="landuse" v="vineyard"/><tag k="barrier" v="fence"/>' elif kul == 61: print ' <tag k="landuse" v="orchard"/><tag k="barrier" v="fence"/>' elif kul == 62: print ' <tag k="landuse" v="orchard"/>' elif kul == 7: print ' <tag k="landuse" v="meadow"/><tag k="meadow" v="agricultural"/>' elif kul == 71: print ' <tag k="landuse" v="meadow"/><tag k="meadow" v="agricultural"/><tag k="barrier" v="fence"/>' elif kul == 72: print ' <tag k="landuse" v="meadow"/><tag k="meadow" v="agricultural"/>' elif kul == 91: print ' <tag k="landuse" v="forest"/><tag k="barrier" v="fence"/>' elif kul == 92: print ' <tag k="landuse" v="farm"/><tag k="crop" v="vegetables"/>' elif kul == 97: print ' <tag k="landuse" v="reservoir"/>' elif kul == 98: print ' <tag k="landuse" v="scrub"/><tag k="note" v="rychle rostouci dreviny"/>' elif kul == 99: print ' <tag k="landuse" v="forest"/>' else: print ' <tag k="landuse" v="unknown_farmland_%d"/>' % kul print ' <tag k="ele" v="%d"/>' % attr(shrec, 'VYSKA') print ' <tag k="ref" v="%s"/>' % attr(shrec, 'NKODFB') for pt in pts: print ' <nd ref="%d"/>' % pt print '</way>' print '</osm>'
_______________________________________________ Talk-cz mailing list [email protected] https://lists.openstreetmap.org/listinfo/talk-cz
