Topologie

PostGIS vznikl jako projekt implementující standard OGC Simple Features pro práci s vektorovými daty ve formě jednoduchých prvků. Od verze 2.0 nicméně umožnuje práci s vektorovými daty také v topologické formě.

Rozšíření pro topologická vektorová data nahrajeme příkazem:

CREATE EXTENSION postgis_topology;

Topologický datový model

Rozdíl mezi datovým modelem jednoduchých geoprvků (simple features) a topologickým modelem si ukážeme na případě dvou sousedících polygonů.

Jednoduché prvky

../_images/postgis-polygon-example.png

Obrázek 1: Dva sousedící polygony jako jednoduché geoprvky.

 polygon |              geometrie (WKT)
---------+----------------------------------------------
    A    | POLYGON((100 0,0 0,0 100,100 100,100 0))
    B    | POLYGON((100 0,100 100,200 100,200 0,100 0))

Z výše uvedeného je evidentní, že jsou oba prvky uloženy odděleně. To vede k tomu, že hranice sousedících polygonů je uložena dvakrát. Jednou jako součást polygonu A podruhé jako součást polygonu B.

Topologický datový model

PostGIS Topology používá datový model Topo-Geo, který vychází z technické normy SQL/MM (ISO 13249-3:2006).

Model pracuje s třemi základními topologickými primitivy:

  • uzly (nodes)
  • hranami (edges)
  • stěnami (faces)

Kompozice znázorněná na obrázku výše bude v topologickém modelu PostGISu popsána:

  • dvěma uzly N1 a N2
  • třemi hranami E1, E2 a E3
  • dvěma stěnami F1 a F2
../_images/postgis-topo-polygon-example.png

Obrázek 2: Dva sousedící polygony z obr. 1 v topologickém modelu.

Ve výsledku bude tedy společná hranice polygonů A a B uložena pouze jednou a to jako hrana E1.

Tip

Podrobné informace k tomuto tématu zde.

Příklad

-- vytvoříme pracovní schéma a nastavíme vyhledávací cestu
CREATE schema topo_test;
-- schéma topology a public musí být v cestě uvedeno vždy
SET search_path TO topo_test,topology,public;

-- nahrání dat ve formě simple features
CREATE TABLE p2 (fid serial PRIMARY KEY, geom geometry(Polygon));
INSERT INTO p2 (geom) VALUES (ST_GeomFromText('Polygon(
 (0 0, 100 0, 100 100, 0 100, 0 0))'));
INSERT INTO p2 (geom) VALUES (ST_GeomFromText('Polygon(
 (100 0, 200 0, 200 100, 100 100, 100 0))'));

Každá datová vrstva s topologii je uložena ve zvláštním schématu, nové schéma vytvoříme pomocí funkce CreateTopology.

SELECT CreateTopology('topo_p2');

Tip

Topologická schéma jsou uložena v tabulce topology (schéma topology).

Do tohoto schématu vložíme nový atribut, do kterého posléze sestavíme topologii prvků. K tomu použijeme funkce AddTopoGeometryColumn.

SELECT AddTopoGeometryColumn('topo_p2', 'topo_test', 'p2', 'topo', 'POLYGON');

Ve výsledku se v tabulce p2 vytvoří nový sloupce s názvem topo a datovým typem TopoGeometry.

Tip

Atributy s topologií jsou uloženy v tabulce layer (schéma topology).

Topologická primitiva sestavíme z jednoduchým prvků pomocí funkce toTopoGeom.

UPDATE p2 SET topo = toTopoGeom(geom, 'topo_p2', 1);

Poznámka

Poslední argument určuje toleranci se kterou budeme topologii sestavovat. Zde jsme zvolili toleranci 1 metr.

Datový typ TopoGeometry

Datový typ TopoGeometry reprezentuje geometrii definovanou topologickými primitivy. Je složen ze čtyř složek:

  • topology_id (id topologického schématu v tabulce topology)
  • layer_id (id topologického atributu v tabulce layer)
  • id (id topologického primitiva)
  • type (geometrický typ jednoduchého prvku)
  • 1 bod (point)
  • 2 lomená čára (linestring)
  • 3 polygon
SELECT fid,ST_AsText(geom),topo FROM p2;

V našem případě:

 fid |                  st_astext                   |   topo
-----+----------------------------------------------+-----------
   1 | POLYGON((0 0,100 0,100 100,0 100,0 0))       | (1,1,1,3)
   2 | POLYGON((100 0,200 0,200 100,100 100,100 0)) | (1,1,2,3)

Tabulky s topologickými primitivy

Topologická primitiva jsou uloženy v tabulkách topologického schématu node, edge a face.

-- seznam uzlů
SELECT node_id,containing_face,st_astext(geom) from topo_p2.node;

-- seznam hran
SELECT edge_id,start_node,end_node,next_left_edge,next_right_edge,
 left_face,right_face,st_astext(geom) from topo_p2.edge;

-- seznam stěn
SELECT face_id,ST_AsText(mbr) from topo_p2.face;

Kontrola konzistence dat

Pro kontrolu topologické konzistence můžete použít dvě funkce TopologySummary a ValidateTopology. První z nich vypisuje souhrné informace o topologii, druhá provádí validaci topologických primitiv.

SELECT TopologySummary('topo_p2');
SELECT ValidateTopology('topo_p2');

Praktická ukázka

Z důvodu časové náročnosti si topologii sestavíme pouze na vzorku parcel na uzemí Hlavního města Prahy.

-- vybereme část parcel na území Hl. města Prahy
CREATE TABLE parcely_732583 AS
 SELECT * FROM ruian_praha.parcely WHERE katastralniuzemikod = 732583;

-- přídáme primární klíč
 ALTER TABLE parcely_732583 ADD PRIMARY KEY (ogc_fid);

-- a prostorové indexy
CREATE INDEX ON parcely_732583 USING gist (geom);

Vytvoříme nové schéma a atribut pro topologii.

-- topologické schéma
SELECT CreateTopology('topo_parcely_732583', 5514);

-- topologický atribut
SELECT AddTopoGeometryColumn('topo_parcely_732583', 'topo_test',
 'parcely_732583', 'topo', 'POLYGON');

Tip

Souřadnicový systém pro topologické schéma můžeme odvodit dynamicky pomocí funkce find_srid, např. find_srid('topo_test', 'parcely_732583', 'geom').

Nakonec se pokusíme topologii sestavit z naimportovaných jednoduchých prvků.

UPDATE parcely_732583 SET topo = toTopoGeom(geom, 'topo_parcely_732583', 1);

Poznámka

Sestavení topologie z jednoduchých geoprvků je poměrně časově náročná činnost. Na výše uvedeném katastrálním území může trvat až několik minut. Funkce toTopoGeom je navíc velmi náchylná na topologické chyby na vstupu a často skončí chybou.

Poznámka pro pokročilé

Pro sestavení topologie můžete použít jako externí nástroj GRASS GIS. Následuje zkracený návod. Detaily tohoto řešení jsou nad rámec tohoto školení a spadají spíše do školení GRASS GIS pro pokročilé.

v.in.ogr in=PG:dbname=pokusnik layer=ukol_1.parcely out=parcely
v.out.postgis -l in=parcely out=PG:dbname=pokusnik out_layer=parcely_topo

Zadání

Najděte parcely, které sousedí s parcelou, ve které se nachází adresní bod s označením kod=22560840.

Řešení

WITH original_face_id AS(
   SELECT
   topology.getFaceByPoint('topo_parcely_732583', geom, 0) face_id
   FROM ruian_praha.adresnimista WHERE kod = 22560840
)
, sousedni AS (
   SELECT
   CASE right_face
      WHEN original_face_id.face_id THEN left_face
      ELSE right_face
   END face
   FROM original_face_id
   , topology.ST_GetFaceEdges('topo_parcely_732583', face_id)
   JOIN topo_parcely_732583.edge ON edge_id = @edge
)
SELECT *
FROM topo_test.parcely_732583 p
JOIN topo_parcely_732583.relation r
ON (p.topo).id = r.topogeo_id
AND (p.topo).layer_id = r.layer_id
AND (p.topo).type = 3
WHERE r.element_id IN (
   SELECT face
   FROM sousedni
);

Praktická ukázka Brloh u Drhovle

Topologickou funkcionalitu PostGISu můžeme s výhodou využít pro zaplochování a tvorbu polygonů u dat, kde máme síť hranic a centroidy. Typickým příkladem dat, která takto dostáváme jsou data ve formátech VFK a VKM vydávané Českým úřadem zeměměřickým a katastrálním.

Další velice užitečnou možností využití funkcionalit topologie v PostGISu je topologická generalizace.

Zadání

Stáhneme vstupní soubor VKM Brlohu u Drhovle a pomocí skriptu v Bashi vybereme linie hranic parcel a definiční body. Data jsou připravena k nahrání do schématu brloh_data.

Sestavíme z dat polygony parcel. Vytvoříme jednoduchou geometrii a topologickou geometrii. Dále provedeme generalizaci hranic parcel a porovnáme výsledek generalizace jednoduché a topologické geometrie.

Postup

Nahrání linií do topologického schématu

Nejdříve vytvoříme nové topologické schéma. Nazveme ho brloh.

SELECT topology.CreateTopology('brloh', 5514, 0.01, false);

Zvolíme SRID "5514" a centimetrovou přesnost. Poslední parametr funkce CreateTopology udává, že ve schématu nepočítáme se souřadnicí Z.

Poznámka

Všiměte si, že funkce nemá prefix "ST", není tedy součástí standardu ISO SQL M/M.

Přidáme všechny linie z tabulky brloh_data.hranice do topologie.

SELECT topology.TopoGeo_AddLineString('brloh', geom, 0) from brloh_data.hranice;

Zaplochujeme

SELECT topology.Polygonize('brloh');

Poznámka

V novějších verzích PostGIS se stěny vytvoří automaticky.

Nyní můžeme vytvořit pohled parcely, který přidá k definičním bodům polygony.

CREATE VIEW parcely AS
SELECT debo.id, debo.label
, topology.ST_GetFaceGeometry(
   'brloh', topology.GetFaceByPoint('brloh', geom, 0)
) geom FROM debo;

Případně můžeme geometrii k bodům doplnit přímo do tabulky.

ALTER TABLE debo ADD geom_polygon geometry(POLYGON, 5514);
UPDATE debo
SET geom_polygon = topology.ST_GetFaceGeometry(
   'brloh'
   , topology.GetFaceByPoint('brloh', geom, 0)
);

Z již vytvořené topologie můžeme vytvořit zpětně topogeometrie bez toho, abychom museli topografický sloupec vytvářet z již získané "jednoduché" geometrie.

Nejdřív přidáme nový atribut pro datový typ TopoGeometry.

SELECT topology.AddTopoGeometryColumn(
   'brloh'
   , 'brloh_data'
   , 'debo'
   , 'topo'
   , 'POLYGON'
);

Číslo vrstvy zjistíme dotazem v tabulce topology.layer.

SELECT layer_id
FROM topology.layer
WHERE schema_name = 'brloh_data' AND table_name = 'debo';

K vygenerování topogeometrie použijeme funkci CreateTopoGeom.

UPDATE debo
SET topo = topology.CreateTopoGeom(
   'brloh'
   , 3 --polygonová vrstva
   , 1 --id vrstvy
   , ARRAY[
      ARRAY[
         topology.GetFaceByPoint('brloh', geom, 0) --face_id
         , 3 --značí, že jde o face
      ]
   ]::topology.TopoElementArray
);

Nyní můžeme zobrazit geometrii sestavenou z topologických primitiv.

SELECT id, label, topo::geometry FROM brloh_data.debo;

Případně se přesvědčit, zda se geometrie sestavená z topologie shoduje s dříve uloženou geometrií vycházející z jednoduchých prvků.

SELECT count(*)
FROM debo
WHERE NOT ST_Equals(geom_polygon, topo::geometry);

Generalizace

Nejdříve vyzkoušíme generalizaci jednoduché geometrie. Použijeme funkci ST_SimplifyPreserveTopology, která na rozdíl od ST_Simplify dohlíží i na validitu vrácených prvků.

SELECT id, st_simplifypreservetopology(geom_polygon, 5) geom
FROM brloh_data.debo;

Výsledek si mohu zobrazit v QGIS.

../_images/simplify.png

Obrázek 3: Výsledek generalizace jednoduchých prvků zobrazený v QGIS.

Z obrázku je celkem jasně patrné, že v důsledku toho, že každý prvek byl generalizován samostatně vznikla řada překryvů a mezer.

Vyzkoušíme to samé v topologii pomocí funkce TP_ST_Simplify.

SELECT id, topology.st_simplify(topo, 15) geom FROM brloh_data.debo;
../_images/simplify_topo.png

Obrázek 4: Výsledek generalizace s ohledem na topologii zobrazený v QGIS.

Vrácený výsledek naprosto přesně odpovídá nezjednodušeným geometriím. Je to proto, že se generalizují jednotlivé hrany a zde použité hrany mají všechny právě dva body a nelze je tedy zjednodušit. Je to důsledek importu z VKM.

Smažeme tedy celé pracně vytvořené topologické schéma. A vytvoříme ho znova.

SELECT topology.DropTopology('brloh');
SELECT topology.CreateTopology('brloh', 5514, 0.01, false);

Tentokrát linie spojíme a rozpustíme pomocí funkce ST_Dump. Navíc je třeba použít funkci ST_LineMerge, která propojí jednotlivé úseky v rámci multilinie.

SELECT topology.TopoGeo_AddLineString('brloh', geom, 0)
FROM (
   SELECT (ST_Dump(u_geom)).geom
   FROM (
      SELECT ST_LineMerge(ST_Union(geom)) u_geom
      FROM brloh_data.hranice
   ) uni
) geoms;

Dále budeme postupovat stejně.

../_images/simplify_topo_2.png

Obrázek 5: Opravený výsledek generalizace na základě topologie zobrazený v QGIS.

Původní hrany jsou vyznačeny světlou linkou.