4.1 Úvod do zpracování vektorových dat pomocí knihovny Fiona¶

V této lekci budeme pracovat s otevřenými daty Data250 z ČÚZK (dokumentace).

Čtení dat¶

In [1]:

import fiona
parks = fiona.open('/tmp/DATA250_2025_SHP/ParkA.shp')
parks

import fiona parks = fiona.open('/tmp/DATA250_2025_SHP/ParkA.shp') parks

Out[1]:

<open Collection '/tmp/DATA250_2025_SHP/ParkA.shp:ParkA', mode 'r' at 0x7f6fc036d400>

Vypíšeme informace o datovém zdroji.

In [2]:

parks.driver, parks.path, parks.bounds

parks.driver, parks.path, parks.bounds

Out[2]:

('ESRI Shapefile',
 '/tmp/DATA250_2025_SHP/ParkA.shp',
 (-891816.2666999996,
  -1209930.3480999991,
  -440147.6123999993,
  -943089.5646000016))

Zaměříme se na souřadnicový systém:

In [3]:

parks.crs, parks.crs_wkt

parks.crs, parks.crs_wkt

Out[3]:

(CRS.from_epsg(5514),
 'PROJCS["S-JTSK / Krovak East North",GEOGCS["S-JTSK",DATUM["System_of_the_Unified_Trigonometrical_Cadastral_Network",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6156"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4156"]],PROJECTION["Krovak"],PARAMETER["latitude_of_center",49.5],PARAMETER["longitude_of_center",24.8333333333333],PARAMETER["azimuth",30.2881397527778],PARAMETER["pseudo_standard_parallel_1",78.5],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","5514"]]')

Kompletní informace:

In [4]:

parks.meta

parks.meta

Out[4]:

{'driver': 'ESRI Shapefile',
 'schema': {'properties': {'FCSUBTYPE': 'int:10',
   'F_CODE': 'str:5',
   'ICC': 'str:5',
   'NAMN1': 'str:50',
   'NA3': 'str:7',
   'NAMN2': 'str:50',
   'NLN1': 'str:3',
   'NLN2': 'str:3',
   'GLOBALID': 'str:38',
   'SHAPE_Leng': 'float:19.11',
   'SHAPE_Area': 'float:19.11',
   'SND': 'int:10',
   'ANND': 'int:10'},
  'geometry': 'Polygon'},
 'crs': CRS.from_epsg(5514),
 'crs_wkt': 'PROJCS["S-JTSK / Krovak East North",GEOGCS["S-JTSK",DATUM["System_of_the_Unified_Trigonometrical_Cadastral_Network",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6156"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4156"]],PROJECTION["Krovak"],PARAMETER["latitude_of_center",49.5],PARAMETER["longitude_of_center",24.8333333333333],PARAMETER["azimuth",30.2881397527778],PARAMETER["pseudo_standard_parallel_1",78.5],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","5514"]]'}

Počet prvků:

In [5]:

len(parks)

len(parks)

Out[5]:

38

K prvkům lze přistupovat náhodně na základě indexu:

In [6]:

feat = parks[0]
feat

feat = parks[0] feat

Out[6]:

fiona.Feature(geometry=fiona.Geometry(coordinates=[[(-615978.6818999983, -994555.0846999995), ...]], type='Polygon'), id='0', properties=fiona.Properties(FCSUBTYPE=2, F_CODE='FA081', ICC='CZ', NAMN1='Broumovsko', NA3='CAT V', NAMN2='UNK', NLN1='cze', NLN2='UNK', GLOBALID='{E2275BAB-7CF6-47A2-B5C5-3F89204BDF6E}', SHAPE_Leng=111219.544184, SHAPE_Area=432206596.814, SND=5220000, ANND=7120200))

Prkvek má dvě části:

• geometrii (geometry)
• atributy (properties)

In [7]:

feat.geometry, feat.properties
#### nebo
# feat["geometry"], feat["properties"]

feat.geometry, feat.properties #### nebo # feat["geometry"], feat["properties"]

Out[7]:

(fiona.Geometry(coordinates=[[(-615978.6818999983, -994555.0846999995), ...]], type='Polygon'),
 fiona.Properties(FCSUBTYPE=2, F_CODE='FA081', ICC='CZ', NAMN1='Broumovsko', NA3='CAT V', NAMN2='UNK', NLN1='cze', NLN2='UNK', GLOBALID='{E2275BAB-7CF6-47A2-B5C5-3F89204BDF6E}', SHAPE_Leng=111219.544184, SHAPE_Area=432206596.814, SND=5220000, ANND=7120200))

In [8]:

feat.geometry.type, len(feat.geometry.coordinates[0])

feat.geometry.type, len(feat.geometry.coordinates[0])

Out[8]:

('Polygon', 1289)

In [9]:

feat.properties["NAMN1"]

feat.properties["NAMN1"]

Out[9]:

'Broumovsko'

In [10]:

for c, v in feat.properties.items():
    print(c, "=>", v)

for c, v in feat.properties.items(): print(c, "=>", v)

FCSUBTYPE => 2
F_CODE => FA081
ICC => CZ
NAMN1 => Broumovsko
NA3 => CAT V
NAMN2 => UNK
NLN1 => cze
NLN2 => UNK
GLOBALID => {E2275BAB-7CF6-47A2-B5C5-3F89204BDF6E}
SHAPE_Leng => 111219.544184
SHAPE_Area => 432206596.814
SND => 5220000
ANND => 7120200

Prvky můžeme procházet sekvenčně:

In [11]:

for f in parks[:3]:
    print(f.geometry.type, f.properties["NAMN1"], f.properties["NA3"])

for f in parks[:3]: print(f.geometry.type, f.properties["NAMN1"], f.properties["NA3"])

Polygon Broumovsko CAT V
Polygon Krkonošský národní park CAT V
Polygon Šumava CAT V

Datový zdroj nezapomeneme zavřít.

In [12]:

parks.close()

parks.close()

Čtení datového zdroje obsahujícího více vrstev¶

V tomto případě použijeme Geonames v režimu otevřených dat šířený ve formátu GeoPackage (otevřený OGC formát založený na souborové databázi SQLite).

In [13]:

fn_gpkg = "/tmp/GEONAMES_RESULTS.gpkg"

fiona.listlayers(fn_gpkg)[:5]

fn_gpkg = "/tmp/GEONAMES_RESULTS.gpkg" fiona.listlayers(fn_gpkg)[:5]

Out[13]:

['MestaObce',
 'CastiMestObci',
 'MistniCasti',
 'ChranenaUzemiPrirody',
 'Regiony']

In [14]:

vm = fiona.open(fn_gpkg, layer="Regiony")
vm.meta

vm = fiona.open(fn_gpkg, layer="Regiony") vm.meta

Out[14]:

{'driver': 'GPKG',
 'schema': {'properties': {'jmeno': 'str:100',
   'hranicni_jmeno': 'str:100',
   'jazyk': 'str:30',
   'zdroj_hr': 'str:20',
   'druh_kod': 'str:20',
   'druh': 'str:100',
   'kat_id': 'str:4',
   'kat_nazev': 'str:50',
   'rotace': 'int16',
   'umisteni': 'int16',
   'typbodu': 'int16',
   'vyslovnost': 'str:255',
   'datum_stan': 'str:12',
   'zkratka': 'str:20',
   'velpis': 'int16',
   'lokalizace': 'str:255',
   'jmeno_druh_lokalizace': 'str:255',
   'fid_jmeno': 'str:20',
   'datum_aktu': 'str:21'},
  'geometry': 'Point'},
 'crs': CRS.from_epsg(5514),
 'crs_wkt': 'PROJCS["S-JTSK / Krovak East North",GEOGCS["S-JTSK",DATUM["System_of_the_Unified_Trigonometrical_Cadastral_Network",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6156"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4156"]],PROJECTION["Krovak"],PARAMETER["latitude_of_center",49.5],PARAMETER["longitude_of_center",24.8333333333333],PARAMETER["azimuth",30.2881397527778],PARAMETER["pseudo_standard_parallel_1",78.5],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","5514"]]'}

Filtrování (atributový dotaz)¶

In [15]:

rlw = fiona.open('/tmp/DATA250_2025_SHP/RailrdL.shp')
rlw.schema

rlw = fiona.open('/tmp/DATA250_2025_SHP/RailrdL.shp') rlw.schema

Out[15]:

{'properties': {'FCSUBTYPE': 'int:10',
  'F_CODE': 'str:5',
  'ICC': 'str:5',
  'EXS': 'int:10',
  'FCO': 'int:10',
  'GAW': 'int:10',
  'LLE': 'int:10',
  'NAMN1': 'str:50',
  'NLN1': 'str:3',
  'RGC': 'int:10',
  'RRA': 'int:10',
  'RRC': 'int:10',
  'RSD': 'int:10',
  'RSU': 'int:10',
  'TUC': 'int:10',
  'GLOBALID': 'str:38',
  'SHAPE_Leng': 'float:19.11',
  'SND': 'int:10'},
 'geometry': 'LineString'}

Zjistíme seznam unikátních hodnot ve sloupci RRC.

In [16]:

rrc = set()
for f in rlw:
    rrc.add(f.properties['RRC'])
rrc

rrc = set() for f in rlw: rrc.add(f.properties['RRC']) rrc

Out[16]:

{16, 17, 999}

In [17]:

main_rlw = []
for f in rlw:
    if f.properties['RRC'] == 16:
        main_rlw.append(f)
len(main_rlw)

main_rlw = [] for f in rlw: if f.properties['RRC'] == 16: main_rlw.append(f) len(main_rlw)

Out[17]:

4554

In [18]:

# zkracený zápis
main_rlw = list(filter(lambda hw: hw.properties['RRC'] == 16, rlw))
len(main_rlw)

# zkracený zápis main_rlw = list(filter(lambda hw: hw.properties['RRC'] == 16, rlw)) len(main_rlw)

Out[18]:

4554

Zápis dat¶

In [19]:

out_fn = '/tmp/main_rlw.gpkg'
out = fiona.open(out_fn, "w", driver="GPKG", crs=rlw.crs, schema=rlw.schema)
out.writerecords(main_rlw)

out_fn = '/tmp/main_rlw.gpkg' out = fiona.open(out_fn, "w", driver="GPKG", crs=rlw.crs, schema=rlw.schema) out.writerecords(main_rlw)

Vysledek si zkontrolujeme pomocí GeoPandas:

In [20]:

import geopandas as gpd

gpd.read_file(out_fn).explore()

import geopandas as gpd gpd.read_file(out_fn).explore()

Out[20]:

Make this Notebook Trusted to load map: File -> Trust Notebook