[Talk-cz] Import of farmland from LPIS

[Martin Švec - OSM]

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
> Talk-cz na openstreetmap.org
> https://lists.openstreetmap.org/listinfo/talk-cz

------------- další část ---------------
#!/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>'