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).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
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='http://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, viz dokumentace.

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
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='http://vdp.cuzk.cz/vymenny_format/soucasna/{}_OB_{}_UKSH.xml.zip'.format(datum, obec)

count_dif = 0

with fiona.open('/vsizip/vsicurl/' + url, 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:
        for f in ds:
            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
                }
            })

            count_dif += 1

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