Ú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]:
Copied!
import fiona
parks = fiona.open('ParkA.shp')
parks
import fiona
parks = fiona.open('ParkA.shp')
parks
Out[1]:
<open Collection '/tmp/uzpr/ParkA.shp:ParkA', mode 'r' at 0x7f71077f65f0>
Vypíšeme informace o datovém zdroji.
In [2]:
Copied!
parks.driver, parks.path, parks.bounds
parks.driver, parks.path, parks.bounds
Out[2]:
('ESRI Shapefile', '/tmp/uzpr/ParkA.shp', (-891816.2666999996, -1209930.3480999991, -440147.6123999993, -943089.5646000016))
Zaměříme se na souřadnicový systém:
In [3]:
Copied!
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]:
Copied!
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', '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]:
Copied!
len(parks)
len(parks)
Out[5]:
38
K prvkům lze přistupovat náhodně na základě indexu:
In [6]:
Copied!
feat = parks[0]
feat
feat = parks[0]
feat
Out[6]:
<fiona.model.Feature at 0x7f71077f71f0>
Prkvek má dvě části:
- geometrii (geometry)
- atributy (properties)
In [7]:
Copied!
feat.geometry, feat.properties
#### nebo
# feat["geometry"], feat["properties"]
feat.geometry, feat.properties
#### nebo
# feat["geometry"], feat["properties"]
Out[7]:
(<fiona.model.Geometry at 0x7f71077f69e0>, <fiona.model.Properties at 0x7f71077f6d70>)
In [8]:
Copied!
feat.geometry.type, len(feat.geometry.coordinates[0])
feat.geometry.type, len(feat.geometry.coordinates[0])
Out[8]:
('Polygon', 1289)
In [9]:
Copied!
feat.properties["NAMN1"]
feat.properties["NAMN1"]
Out[9]:
'Broumovsko'
In [10]:
Copied!
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 SHAPE_Leng => 111219.544184 SHAPE_Area => 432206596.814 SND => 5220000 ANND => 7120200
Prvky můžeme procházet sekvenčně:
In [11]:
Copied!
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]:
Copied!
parks.close()
parks.close()
In [13]:
Copied!
fn_gpkg = "Geonames.gpkg"
fiona.listlayers(fn_gpkg)[:5]
fn_gpkg = "Geonames.gpkg"
fiona.listlayers(fn_gpkg)[:5]
Out[13]:
['MestaObce', 'CastiMestObci', 'MistniCasti', 'ChranenaUzemiPrirody', 'Regiony']
In [14]:
Copied!
vm = fiona.open(fn_gpkg, layer="Regiony")
vm.meta
vm = fiona.open(fn_gpkg, layer="Regiony")
vm.meta
Out[14]:
{'driver': 'GPKG', 'schema': {'properties': {'id': 'float', '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': 'int', 'umisteni': 'int', 'typbodu': 'int', 'vyslovnost': 'str:255', 'datum_stan': 'str:12', 'zkratka': 'str:20', 'velpis': 'int', 'lokalizace': 'str:255', 'jmeno_druh_lokalizace': 'str:255', 'fid_jmeno': 'str:20', 'datum_aktu': 'str:20'}, '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 [21]:
Copied!
rlw = fiona.open('/tmp/uzpr/RailrdL.shp')
rlw.schema
rlw = fiona.open('/tmp/uzpr/RailrdL.shp')
rlw.schema
Out[21]:
{'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', 'SHAPE_Leng': 'float:19.11', 'SND': 'int:10'}, 'geometry': 'LineString'}
Zjistíme seznam unikátních hodnot ve sloupci RRC
.
In [22]:
Copied!
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[22]:
{16, 17, 999}
In [23]:
Copied!
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[23]:
3044
In [18]:
Copied!
# 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]:
3044
Zápis dat¶
In [24]:
Copied!
out_fn = 'main_rlw.gpkg'
out = fiona.open(out_fn, "w", driver="GPKG", crs=rlw.crs, schema=rlw.schema)
out.writerecords(main_rlw)
out_fn = '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]:
Copied!
import geopandas as gpd
gpd.read_file(out_fn).explore()
import geopandas as gpd
gpd.read_file(out_fn).explore()
/usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.0 warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
Out[20]:
Make this Notebook Trusted to load map: File -> Trust Notebook