Komplexní příklad přístupu k datům RÚIAN

Tato kapitola prezentuje příklad použití knihovny GDAL pro čtení online dat z Registru územní identifikace, adres a nemovitostí (RÚIAN) dostupných v rámci Veřejného dálkové přístupu.

Pro zadanou obec (řádek 8) skript vypíše informace o parcelách (parcelní číslo, výměru uloženou v atributech a vypočtenou z geometrie). Parcely, u kterých je rozdíl výměr větší než tolerance (řádek 5), uloží do nové vrstvy ve formátu Esri Shapefile do aktuálního adresáře (řádek 27).

import datetime
from osgeo import ogr, osr

# tolerance ve vymere
tol = 100

# kod zajmove obce
obec = 505528

# posledni datum v mesici (k tomuto dni jsou publikovana stavova data)
today = datetime.date.today()
if today.month == 12:
    day = today.replace(day=31)
day = (today.replace(month=today.month, day=1) - datetime.timedelta(days=1))
datum = day.strftime("%Y%m%d")

# URL souboru VDP
url='https://vdp.cuzk.cz/vymenny_format/soucasna/{}_OB_{}_UKSH.xml.zip'.format(datum, obec)

# otevrit vstupni datasource (RUIAN)
ds = ogr.Open('/vsizip/vsicurl/' + url)
# nacist vrstvu parcel
layer = ds.GetLayerByName('Parcely')

# vytvorit vystupni datasource (Shapefile - parcely)
driver = ogr.GetDriverByName("ESRI Shapefile")
dso = driver.CreateDataSource("parcely.shp")
srs = osr.SpatialReference()
srs.ImportFromEPSG(5514)
olayer = dso.CreateLayer("parcely", srs, ogr.wkbPolygon)

# definovat atributovou tabulku vystupniho Shapefile
field_name = ogr.FieldDefn("cislo", ogr.OFTString)
field_name.SetWidth(255)
olayer.CreateField(field_name)
olayer.CreateField(ogr.FieldDefn("vymera", ogr.OFTReal))
olayer.CreateField(ogr.FieldDefn("rozdil", ogr.OFTReal))

# pocet parcel
count = layer.GetFeatureCount()

# pocet parcel nad toleranci
count_dif = 0

# sekvencne cist parcely
layer.ResetReading()
while True:
    feat = layer.GetNextFeature()
    if feat is None:
        break
		
    # nacist zajmove atributy
    kc = feat.GetField('KmenoveCislo')
    pod = feat.GetField('PododdeleniCisla')
    vym = feat.GetField('VymeraParcely')
	
    # prvek ma vice geometrickych reprezentaci, zajima nas hranice parcel
    geom = feat.GetGeomFieldRef('OriginalniHranice')
	
    # na zakladne geometrie spocitat vymeru
    vym2 = geom.GetArea()
	
    # vypocitat rozdil mezi vymerou z atributove tabulky a geometrie
    dif = abs(vym - vym2)

    # oznaceni parcely
    pc = str(kc)
    if pod:
        pc += '/' + str(pod)
	
    # vytisknout report pro kazdou parcelu
    print ('{0}: {1} {2:.1f} {3:.1f}'.format(pc, vym, vym2, dif))
	
    if dif < tol:
        continue
	
    # parcely, u kterych je rozdil vymer vetsi nez 100m2, ulozime
    # do vystupniho shapefile
    ofeature = ogr.Feature(olayer.GetLayerDefn())
    ofeature.SetField("cislo", pc)
    ofeature.SetField("vymera", vym)
    ofeature.SetField("rozdil", dif)
    ofeature.SetGeometry(geom)
    olayer.CreateFeature(ofeature)
    
    count_dif += 1
    
    # dealokovat pamet
    ofeature = None
    
# vytisknout zaverecnou statistiku
print ('-' * 80)
print ('Pocet parcel: {}'.format(count))
print ('Pocet parcel nad toleranci: {}'.format(count_dif))

# zavrit datove zdroje
ds.Destroy()
dso.Destroy()

Poznámka

Skript ke stažení zde.

Tip

Pro efektivní práci s daty RÚIAN existuje speciální knihovna GDAL-VFR a její konzolové nastroje.

Pro srovnání ještě stejný skript za použití knihoven Fiona a Shapely

import datetime
import fiona
from shapely.geometry import shape

# tolerance ve vymere
tol = 100

# kod zajmove obce
obec = 505528

# posledni datum v mesici (k tomuto dni jsou publikovana stavova data)
today = datetime.date.today()
if today.month == 12:
    day = today.replace(day=31)
day = (today.replace(month=today.month, day=1) - datetime.timedelta(days=1))
datum = day.strftime("%Y%m%d")

# URL souboru VDP
url = 'https://vdp.cuzk.cz/vymenny_format/soucasna/{}_OB_{}_UKSH.xml.zip'.format(datum, obec)

# Fiona nepodporuje zdroje s vice geometriemi, proto data nejprve predzpracujeme
# pomoci knihovny GDAL
from osgeo import gdal
parcely_poly = "/tmp/parcely.gml"
gdal.VectorTranslate(parcely_poly, '/vsizip/vsicurl/' + url, format="GML",
                     SQLStatement="select OriginalniHranice,* from parcely")
# with fiona.open('/vsizip/vsicurl/' + url, layer='Parcely') as ds:
with fiona.open(parcely_poly, layer='Parcely') as ds:
    schema = {
            "geometry": "Polygon",
            "properties": {
                "cislo":  "str:25",
                "vymera": "float",
                "rozdil": "float"
                }
            }
    with fiona.open("parcely.shp", 'w', driver="ESRI Shapefile", crs="epsg:5514", schema=schema) as dso:
        count = 0
        count_dif = 0

        for f in ds:
            count += 1
            
            cislo = f["properties"]["KmenoveCislo"]
            if f["properties"]["PododdeleniCisla"]:
                cislo  = "{}/{}".format(cislo, f["properties"]["PododdeleniCisla"])

            vym = float(f["properties"]["VymeraParcely"])
            vym2 = shape(f["geometry"]).area
            dif = abs(vym - vym2)
            if dif < tol:
                continue

            dso.write({
                "type": "Feature",
                "properties": {
                    "cislo": cislo,
                    "vymera": f["properties"]["VymeraParcely"],
                    "rozdil": vym - vym2
                },
                "geometry": f["geometry"]
            })

            count_dif += 1

        print("-"*80)
        print("Počet parcel: {}".format(count))
        print("Počet parcel nad tolerancí: {}".format(count_dif))