4.2 Úvod do knihovny Shapely¶

Navážeme na předchozí cvičení a otevřeme vrstvu s hlavními železničními tratěmi.

In [1]:

import fiona
main_rlw = fiona.open('/tmp/out.gpkg')
len(main_rlw)

import fiona main_rlw = fiona.open('/tmp/out.gpkg') len(main_rlw)

Out[1]:

4554

Pokud chceme pracovat s geometrickou sloužkou popisu, tak potřebujeme knihovnu Shapely. Prvním úkolem bude spočitát celkovou délku hlavních tratí v ČR. Ukázka na prvním segmentu:

In [2]:

from shapely.geometry import shape
# pozor: GPKG čísluje prvky od 1
shape(main_rlw[1].geometry).length

from shapely.geometry import shape # pozor: GPKG čísluje prvky od 1 shape(main_rlw[1].geometry).length

Out[2]:

233.59242604638075

Geometrii ve formě shape lze i jednoduše vykreslit:

In [3]:

shape(main_rlw[1].geometry)

shape(main_rlw[1].geometry)

Out[3]:

In [4]:

sum_length = 0
for f in main_rlw:
    sum_length += shape(f.geometry).length

print(f"Délka v km: {round(sum_length / 1000)}")

sum_length = 0 for f in main_rlw: sum_length += shape(f.geometry).length print(f"Délka v km: {round(sum_length / 1000)}")

Délka v km: 3802

Transformace mezi souřadnicovými systémy¶

In [5]:

from fiona.transform import transform_geom
from fiona.crs import from_epsg

target_crs = from_epsg(4326)

geom_4326 = transform_geom(
    main_rlw.crs, target_crs, main_rlw[1].geometry
)
geom_4326.coordinates

from fiona.transform import transform_geom from fiona.crs import from_epsg target_crs = from_epsg(4326) geom_4326 = transform_geom( main_rlw.crs, target_crs, main_rlw[1].geometry ) geom_4326.coordinates

Out[5]:

[(17.393853580516584, 49.61166705199342),
 (17.395362460117184, 49.61352454611004)]

Spočítejme prostorový rozsah dat ve WGS-84:

In [6]:

shape(geom_4326).envelope

shape(geom_4326).envelope

Out[6]:

In [7]:

main_rlw.bounds

main_rlw.bounds

Out[7]:

(-896176.085999999, -1217278.3142, -438599.808200002, -948355.5218)

In [8]:

from shapely.ops import unary_union

all_rlw = unary_union(
    [shape(transform_geom(main_rlw.crs, target_crs, f.geometry)) for f in main_rlw]
)
all_rlw.bounds

from shapely.ops import unary_union all_rlw = unary_union( [shape(transform_geom(main_rlw.crs, target_crs, f.geometry)) for f in main_rlw] ) all_rlw.bounds

Out[8]:

(12.253141281385112, 48.593071341515476, 18.763925446381837, 50.95739586957717)

Tento způsob není zcela optimální. Vylepšená verze:

In [9]:

from shapely.geometry import mapping
from shapely.ops import unary_union

mall_rlw = transform_geom(
    main_rlw.crs,
    target_crs,
    mapping(unary_union([shape(f.geometry) for f in main_rlw]))
)
shape(all_rlw).bounds

from shapely.geometry import mapping from shapely.ops import unary_union mall_rlw = transform_geom( main_rlw.crs, target_crs, mapping(unary_union([shape(f.geometry) for f in main_rlw])) ) shape(all_rlw).bounds

Out[9]:

(12.253141281385112, 48.593071341515476, 18.763925446381837, 50.95739586957717)

Určete, jak ovlivňují hlavní železniční tratě největší CHKO v ČR¶

Nejprve najdeme největší CHKO (NA3 == 'CAT V') v ČR.

Výběr chráněného území¶

In [10]:

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

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

Out[10]:

{'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'}

In [11]:

parks_area = {}
for f in parks:
    if f.properties["NA3"] != "CAT V":
        continue
    name = f.properties['NAMN1']
    if name not in parks_area:
        parks_area[name] = 0
    parks_area[name] += shape(f.geometry).area
        
parks_area

parks_area = {} for f in parks: if f.properties["NA3"] != "CAT V": continue name = f.properties['NAMN1'] if name not in parks_area: parks_area[name] = 0 parks_area[name] += shape(f.geometry).area parks_area

Out[11]:

{'Broumovsko': 432206596.8137291,
 'Krkonošský národní park': 363570935.9759183,
 'Šumava': 995845458.4680469,
 'Labské pískovce': 242104236.44435042,
 'Litovelské Pomoraví': 93419634.2008796,
 'Lužické hory': 270103882.5698882,
 'Jeseníky': 741332206.9236288,
 'Pálava': 85639555.81207283,
 'Bílé Karpaty': 746694026.0388199,
 'Poodří': 83688261.05104716,
 'Český ráj': 182198623.25407416,
 'Blaník': 40456954.83506415,
 'Český kras': 132283622.52276358,
 'Blanský les': 219631393.77284405,
 'Beskydy': 1203345814.8093934,
 'Český les': 466001615.58559895,
 'Orlické hory': 233859151.06964803,
 'České středohoří': 1068605271.770449,
 'Jizerské hory': 374453354.19508487,
 'Třeboňsko': 687380581.1151866,
 'Moravský kras': 97391405.49233718,
 'Železné hory': 284996668.22751456,
 'Křivoklátsko': 625378137.957978,
 'Kokořínsko-Máchův kraj': 410767218.5685966,
 'Žďárské vrchy': 709008240.7698,
 'Brdy': 344998430.0609124,
 'Slavkovský les': 610279695.4088577}

In [12]:

parkname = max(parks_area, key=parks_area.get)

# pokud bychom si chteli procvicit pandas:

# import pandas
#
# parkname = pandas.DataFrame.from_dict(parks_area, orient="index").sort_values(
#     0, ascending=False
# ).iloc[0].name

parkname = max(parks_area, key=parks_area.get) # pokud bychom si chteli procvicit pandas: # import pandas # # parkname = pandas.DataFrame.from_dict(parks_area, orient="index").sort_values( # 0, ascending=False # ).iloc[0].name

In [13]:

sel_parks = []
for f in parks:
    if f.properties["NAMN1"] == parkname:
        sel_parks.append(shape(f.geometry))
sel_park = unary_union(sel_parks)
sel_park

sel_parks = [] for f in parks: if f.properties["NAMN1"] == parkname: sel_parks.append(shape(f.geometry)) sel_park = unary_union(sel_parks) sel_park

Out[13]:

Výběr hlavních tratí¶

In [14]:

sel_rlws = []
for f in main_rlw:
    geom = shape(f.geometry)
    if geom.intersects(sel_park):
        sel_rlws.append(geom)
        
sel_rlw = unary_union(sel_rlws)
sel_rlw

sel_rlws = [] for f in main_rlw: geom = shape(f.geometry) if geom.intersects(sel_park): sel_rlws.append(geom) sel_rlw = unary_union(sel_rlws) sel_rlw

Out[14]:

Vytvoření obalové zóny¶

In [15]:

sel_rlw_200 = sel_rlw.buffer(200)
sel_rlw_200

sel_rlw_200 = sel_rlw.buffer(200) sel_rlw_200

Out[15]:

Průnik obalové zóny s chraněným územím¶

In [16]:

sel_park_rlw_200 = sel_rlw_200.intersection(sel_park)
sel_park_rlw_200

sel_park_rlw_200 = sel_rlw_200.intersection(sel_park) sel_park_rlw_200

Out[16]:

In [17]:

sel_rlw_200.area, sel_park_rlw_200.area

sel_rlw_200.area, sel_park_rlw_200.area

Out[17]:

(9409617.171661964, 3846705.3217628486)

In [18]:

print(f"Celková výměra CHKO (ha): {round(sel_park.area / 1e4)}")
print(f"Zasažená oblast (ha): {round(sel_park_rlw_200.area / 1e4)}")
print(f"Zasažená oblast (%): {round((sel_rlw_200.area / sel_park.area) * 100, 2)}")

print(f"Celková výměra CHKO (ha): {round(sel_park.area / 1e4)}") print(f"Zasažená oblast (ha): {round(sel_park_rlw_200.area / 1e4)}") print(f"Zasažená oblast (%): {round((sel_rlw_200.area / sel_park.area) * 100, 2)}")

Celková výměra CHKO (ha): 120335
Zasažená oblast (ha): 385
Zasažená oblast (%): 0.78

Vše v jednom¶

In [19]:

def report_chko_rlw(parkname, distance=200):
    sel_parks = []
    for f in parks:
        if f.properties["NAMN1"] == parkname:
            sel_parks.append(shape(f.geometry))
    sel_park = unary_union(sel_parks)
    
    sel_rlws = []
    for f in main_rlw:
        geom = shape(f.geometry)
        if geom.intersects(sel_park):
            sel_rlws.append(shape(f.geometry))
        
    sel_rlw = unary_union(sel_rlws)

    sel_rlw_200 = sel_rlw.buffer(distance)

    sel_park_rlw_200 = sel_rlw_200.intersection(sel_park)

    print(f"Celková výměra CHKO (ha): {round(sel_park.area / 1e4)}")
    print(f"Zasažená oblast (ha): {round(sel_park_rlw_200.area / 1e4)}")
    print(f"Zasažená oblast (%): {round((sel_park_rlw_200.area / sel_park.area) * 100, 2)}")
    
report_chko_rlw("Beskydy")

def report_chko_rlw(parkname, distance=200): sel_parks = [] for f in parks: if f.properties["NAMN1"] == parkname: sel_parks.append(shape(f.geometry)) sel_park = unary_union(sel_parks) sel_rlws = [] for f in main_rlw: geom = shape(f.geometry) if geom.intersects(sel_park): sel_rlws.append(shape(f.geometry)) sel_rlw = unary_union(sel_rlws) sel_rlw_200 = sel_rlw.buffer(distance) sel_park_rlw_200 = sel_rlw_200.intersection(sel_park) print(f"Celková výměra CHKO (ha): {round(sel_park.area / 1e4)}") print(f"Zasažená oblast (ha): {round(sel_park_rlw_200.area / 1e4)}") print(f"Zasažená oblast (%): {round((sel_park_rlw_200.area / sel_park.area) * 100, 2)}") report_chko_rlw("Beskydy")

Celková výměra CHKO (ha): 120335
Zasažená oblast (ha): 385
Zasažená oblast (%): 0.32

In [20]:

report_chko_rlw("Beskydy", 500)

report_chko_rlw("Beskydy", 500)

Celková výměra CHKO (ha): 120335
Zasažená oblast (ha): 1003
Zasažená oblast (%): 0.83

In [21]:

report_chko_rlw("České středohoří")

report_chko_rlw("České středohoří")

Celková výměra CHKO (ha): 106861
Zasažená oblast (ha): 4466
Zasažená oblast (%): 4.18