Dálnice přes České středohoří

Úkol

Nová dálnice přes CHKO České středohoří

Najděte oblasti, kterých se bezprostředně dotkla výstavba dálnice D8 z Prahy do Drážďan. Předpokládejme, že se jedná o oblast 200m od dálnice. O jak velké území se jedná?

Obr. 1 Dálniční síť České republiky spolu s mapou chráněných krajinných oblastí.

Potřebné kroky pro vyřešení úlohy

1. Otevření obou datových zdrojů (highways.geojson a chko.shp)
2. Převod na společný souřadnicový systém (první datový zdroj je v EPSG:4326, druhý v EPSG:5514)
3. Vytvoření obalové zóny okolo dálnic
4. Zjištění průniku všech geometrií

Načtení datových zdrojů

Nejprve oba soubory otevřeme pomocí klasického fiona.open(). Použijeme buď klauzuli with nebo nesmíme nakonec zapomenout soubory uzavřít.

import fiona
import json
from shapely.geometry import mapping, shape
from fiona.transform import transform_geom
from shapely.ops import cascaded_union
from shapely import speedups

speedups.enable()

with fiona.open("data/chko.shp", "r") as pas:
    with fiona.open("data/highways.geojson", "r") as hws:

Společný souřadnicový systém

Data, jak můžeme ověřit, jsou v různých souřadnicových systémech.

Úkol

Zjistěte, v jakých souřadnicových systémech jsou dálnice a chráněné krajinné oblasti.

Pro další práci musíme data převést do společného souřadnicového systému.

Pro transformaci vektorových prvků (resp. jejich geometrií) můžeme použít funkci trasnform_geom().

from fiona.transform import transform_geom

wgs84 = "epsg:4326"
jtsk = {"init": "epsg:5514", "towgs84": "570.8,85.7,462.8,4.998,1.587,5.261,3.56"}

geom_transformed = transform_geom(wgs84, jtsk, feature["geometry"])

Vytvoření obalové zóny okolo dálnic

Pro každý prvek v dálniční síti, nejprve vyfiltrujeme pouze prvky se jménem dálnice D8, transformujeme geometrii na S-JTSK a vytvoříme obalovou zónu:

import fiona
import json
from shapely.geometry import mapping, shape
from fiona.transform import transform_geom
from shapely.ops import cascaded_union
from shapely import speedups

speedups.enable()

with fiona.open("data/chko.shp", "r") as pas:
    with fiona.open("data/highways.geojson", "r") as hws:

        buffered_highways = []

        wgs84 = "EPSG:4326"
        jtsk = pas.crs

        d8 = list(filter(lambda hw: hw["properties"]["ref"] == "D8", hws))

        for hw in d8:
            geom_transformed = transform_geom(wgs84, jtsk, hw["geometry"])

            g = shape(geom_transformed)
            buffered_highways.append(g.buffer(100))

Obr. 2 Obalová zóna 200m okolo dálnice

Poznámka

Zrychlení topologických operací

Modul shapely.speedups obsahuje zrychlení některých operací optimalizovaných v jazyce C. Pokud během instalace byla dostupná knihovan GEOS a pokud byl dostupný kompilátor jazyka C, jsou tyto optimalizace automaticky zavedeny.

Ověřit jejich přítomnost, případně je explicitně povolit, můžeme modulem speedups:

>>> from shapely import speedups
>>> print(speedups.available)
True
>>> speedups.enable()

Průnik geometrií

Pro každou chráněnou krajinnou oblast musíme nejprve zjistit, jestli se geometrie CHKO a obalové zóny okolo dálnice protínají - teprve potom má smysl pokračovat s interkací.

import fiona
import json
from shapely.geometry import mapping, shape
from fiona.transform import transform_geom
from shapely.ops import cascaded_union
from shapely import speedups

speedups.enable()

with fiona.open("data/chko.shp", "r") as pas:
    with fiona.open("data/highways.geojson", "r") as hws:

        buffered_highways = []

        wgs84 = "EPSG:4326"
        jtsk = pas.crs

        d8 = list(filter(lambda hw: hw["properties"]["ref"] == "D8", hws))

        for hw in d8:
            geom_transformed = transform_geom(wgs84, jtsk, hw["geometry"])

            g = shape(geom_transformed)
            buffered_highways.append(g.buffer(100))

        intersections = []

        for hw in buffered_highways:

            for pa in pas:

                pa_geom = shape(pa["geometry"])
                if hw.intersects(pa_geom):
                    out_geom = hw.intersection(pa_geom)
                    intersections.append(out_geom)

Spojení do jedné geometrie

V tuto chvíli máme data rozkouskovaná, bylo by ale vhodné je spojit do jednoho geometrického objektu.

Úkol

Spojení geometrií do jednoho objektu

Navrhněte postup, jak spojit více geometrií do jedné.

Zápis do souboru - OGC GeoPackage

Pokud potřebujeme uložit některý z binárních typů vektorových souborů (ESRI Shapefile, OGC GeoPackage, …) musíme nejprve nadefinovat datové schéma - atributy a typ geometrií.

schema = {
        'properties': {
                'highway': 'str'
        },
        'geometry': 'MultiPolygon'
}

Následně můžeme otevřít nový soubor pro zápis a uložit do něj námi vytvořenou geometrii.

with fiona.open("chko_x_highway.gpkg", "w", driver="GPKG", crs="EPSG:5514", schema=schema) as out:
        feature = {
                'type': 'Feature',
                'properties': {
                        'highway': 'D8'
                },
                'geometry': mapping(one_geometry)
        }
        out.write(feature)

Zápis do souboru - GeoJSON

Zapsat soubor ve formátu GeoJSON je od dost jednodušší, protože se nemusíme spoléhat na (relativně) komplikovanou funkci fiona.write(), ale můžeme data zapsat napřímo pomocí modulu json.

To nám samozřejmě umožňuje nedodržet stritkní schéma atributů, což může a nemusí být výhoda.

Mějte na paměti, že podle GeoJSON specifikace se předpokládá, že souřadnicovým systémem je WGS84.

Poznámka

Existuje i způsob, jak do GeoJSONu uložit geometrie v jiných souřadnicový systémech a i když s ním některé softwary počítají, do finální specifikace se tento způsob nedostal.

Úkol

Zapište výstup do formátu GeoJSON

Vytvořte část skriptu, která převede geometrie na WGS84 a zapíše výsledek do souboru GeoJSON.

Tip

Pro převod použijte funkci transform_geom(), povinný „obal“ seznamu vektorových prvků (features) vypadá pro formát GeoJSON následovně.

{
    "type": "FeatureCollection",
    "features": [
            {
                "type":"Feature",
                "properties": {...},
                "geometry": {...}
            },
            {
                ...
            },
                ...

    ]
}

Pokud je pro vás úloha příliš komplikovaná, můžete se na celý skript samozřejmě podívat.

Další atributy

Někdy je výhodné pro další práci s vektorovými daty, uložit vlastnosti geometrie do atributové tabulky, např. rozlohu plochy nebo délku linie.

Úkol

Uložení dalších atributů

Upravíte skript tak, aby výstup obsahoval také celkovou zasaženou plochu (m2) CHKO?

Obr. 3 Obalová zóna okolo dálnice protínající CHKO České středohoří.

Optimalizace (diskuse)

Skript je ve své podstatě celkem neefektivní a stálo by za to popřemýšlet o jeho optimalizaci, aby proběhl rychleji:

1. Dálniční těleso je reprezentováno dvěmi liniemi, stálo by za zvážení použít pouze jednu z nich
2. Buffer by se mohl spojit pomocí cascaded_union() do jedné geometrie, tím by se měl následný průnik zrychlit
3. …