Ú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]:
Copied!
import fiona
main_rlw = fiona.open('main_rlw.gpkg')
len(main_rlw)
import fiona
main_rlw = fiona.open('main_rlw.gpkg')
len(main_rlw)
Out[1]:
3044
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]:
Copied!
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]:
238.743298120115
Geometrii ve formě shape
lze i jednoduše vykreslit:
In [3]:
Copied!
shape(main_rlw[1].geometry)
shape(main_rlw[1].geometry)
Out[3]:
In [4]:
Copied!
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: 2402
Transformace mezi souřadnicovými systémy¶
In [5]:
Copied!
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]:
[(16.868457358462393, 49.19663693512262), (16.871299618777353, 49.19770455779107)]
Spočítejme prostorový rozsah dat ve WGS-84:
In [6]:
Copied!
shape(geom_4326).envelope
main_rlw.bounds
shape(geom_4326).envelope
main_rlw.bounds
Out[6]:
(-896176.085999999, -1217278.3142, -438599.808200002, -955983.221799999)
In [7]:
Copied!
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.envelope.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.envelope.bounds
Out[7]:
(12.253141281385112, 48.593071341515476, 18.76391449120931, 50.85966753785939)
Tento způsob není zcela optimální. Vylepšená verze:
In [8]:
Copied!
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).envelope.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).envelope.bounds
Out[8]:
(12.253141281385112, 48.593071341515476, 18.76391449120931, 50.85966753785939)
In [9]:
Copied!
parks = fiona.open('ParkA.shp')
parks.schema
parks = fiona.open('ParkA.shp')
parks.schema
Out[9]:
{'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'}
In [10]:
Copied!
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[10]:
{'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': 741321209.5754136, 'Pálava': 85639555.81207283, 'Bílé Karpaty': 746694026.0388199, 'Poodří': 83733448.9926418, 'Český ráj': 182198623.25407416, 'Blaník': 40456954.83506415, 'Český kras': 132283622.52276358, 'Blanský les': 219650117.06156796, 'Beskydy': 1203345814.8093934, 'Český les': 466001615.58559895, 'Orlické hory': 233859151.06964803, 'České středohoří': 1068605271.770449, 'Jizerské hory': 374453030.2032833, 'Třeboňsko': 687380581.1151866, 'Moravský kras': 97450000.49400458, 'Ž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 [11]:
Copied!
import pandas
pandas.DataFrame.from_dict(parks_area, orient="index").sort_values(0, ascending=False)
# TODO: variable
parkname = "Beskydy"
import pandas
pandas.DataFrame.from_dict(parks_area, orient="index").sort_values(0, ascending=False)
# TODO: variable
parkname = "Beskydy"
In [12]:
Copied!
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[12]:
Výběr hlavních tratí¶
In [13]:
Copied!
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
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
Out[13]:
Vytvoření obalové zóny¶
In [14]:
Copied!
sel_rlw_200 = sel_rlw.buffer(200)
sel_rlw_200
sel_rlw_200 = sel_rlw.buffer(200)
sel_rlw_200
Out[14]:
Průnik obalové zóny s chraněným územím¶
In [15]:
Copied!
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[15]:
In [16]:
Copied!
sel_rlw_200.area, sel_park_rlw_200.area
sel_rlw_200.area, sel_park_rlw_200.area
Out[16]:
(9409617.171661964, 3846705.3217628486)
In [17]:
Copied!
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 [18]:
Copied!
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_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_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.78
In [19]:
Copied!
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 (%): 2.03
In [20]:
Copied!
report_chko_rlw("České středohoří")
report_chko_rlw("České středohoří")
Celková výměra CHKO (ha): 106861 Zasažená oblast (ha): 3430 Zasažená oblast (%): 3.42