Geometrie geoprvků a knihovna Shapely

V této části se blíže seznámíme jak pracovat s geometrií vektorových prvků.

Knihovna Shapely (stejného autora jako Fiony) nám umožňuje pracovat s geometrickou složkou popisu geoprvků opět ve stylu jazyka Python. Stejně jako Fiona, převádí Shapely geometrické vlastnosti na objekty typu JSON.

Jak jsme si řekli, Fiona převádí vstupní data z různých formátů na strukturu GeoJSON:

from shapely.geometry import shape
cr = chko[54]

Shapely

Shapely obsahuje i některé funkce pro modifikaci geometrií, například generalizaci, obalovou zónu (buffer) nebo porovnání dvou geometrií. Shapely využívá pro manipulaci s geometriemi knihovu GEOS, která je zase re-implementací nástroje JTS (Java Topoology Suite) do jazyka C, a používá se například v databázové nadstavbě PostGIS. GEOS a JTS odpovídají specifikaci OGC Simple Feature. Knihovna Shapely nám umožní řešit například úlohy typu:

• Jakou výměru má plošný prvek?
• Překrývají se dva prvky?
• Jak vypadá společný průnik více prvků?
• Vytvořit obalovou zónu okolo prvku.

Pro práci s geometriemi proto potřebujeme převést data na struktury Shapely:

poly = shape(cr['geometry'])
print (poly.bounds)

simple = poly.simplify(10)
print (simple.intersects(poly))

(-683329.1875, -993228.75, -681265.625, -991528.0)

Poznámka

V prostředí Jupyter notebook můžeme geometrii vizualizovat prostě tak, že obsah proměnné s geometrií necháme vypsat, např.:

geom = shape(chko[54]["geometry"])
geom

Todo

Příklad vykreslení geometrie pomocí matplotlib https://deparkes.co.uk/2015/03/11/how-to-plot-polygons-in-python/

Nyní se můžeme vypsat celou řadu metadat o daném geometrickém objektu

print(geom.type)
print(geom.area)
...

Úkol

Zjistěte, jaký má naše geometrie obvod a najděte souřadnice jejího centroidu. Jaká metoda se použije pro vytvoření obalové zóny a jak funguje?

Poznámka

Kompletní ukázkový příklad na využití knihoven Fiona a Shapely zde.

Převod geometrie zpět na formát JSON

Na převod zpět do formátu JSON použijeme funkci mapping z modulu shapely.geometry

Jako příklad si vytvoříme „ručně“ celý nový vektorový objekt

new_feature = {
    "type": "Feature",
    "properties": {
        "NAZEV": "prvek s obalovou zonou 100m",
    },
    "geometry": mapping(geom.buffer(100))
}