Materiál k cvičení DB2 týden #12, sezona 2024/2025

Geospatial data (PostgreSQL)

Geospaciální data

Geospaciální data (geodata, prostorová data) je typ dat, který se používá k reprezentaci geografických nebo prostorových informací. Typicky tedy obsahují informaci o tom, kde se nějaký záznam nebo událost nachází. Lze si představit:

V aplikacích se s geospaciálními daty lze setkat v různých oblastech nebo use-case, jako jsou:

Specifikem geospaciálních dat je, že se jedná o data, která mají prostorovou složku. To znamená, že kromě klasických atributů (např. jméno, věk) obsahují také informace o poloze (např. zeměpisná šířka a délka). Tato data jsou reprezentována complexními geometrickými objekty jako:

Zároveň pro práci s takovými daty je třeba mít k dispozici i specializované nástroje a databáze, které umožňují efektivní zpracování geometrií a dotazů např. na vzdálenosti, překryvy, sousedství atd.

Příprava prostředí

V rámci cvičení budeme používat databázi PostgreSQL (v14.17) s rozšířením PostGIS (v3.2.0), která momentálně a dočasně běží na mém VPS pro zjednodušení a urychlení práce.

V public schematu jsou připravené dvě tabulky cities a countries, které obsahují geospaciální data o městech a zemích světa.

Vlastní instalace

Pokud chcete rozjet vlastní instalaci a nechcete používat docker, já jsem postupoval následovně.

Nainstalovat PostgreSQL a PostGIS (v tomto případě jednak jako rozšíření databáze pro běžnou práci a jednak jako samostatný balíček pro předzpracování dat při načítání):


  

sudo apt install postgresql postgresql-contrib
sudo apt install postgresql-postgis
sudo apt install postgis

Zkontrolovat, že postgreSQL běží:


  

sudo systemctl status postgresql
sudo systemctl enable postgresql

Přihlásit se do databáze jako postgres uživatel, přepnout na existující databázi postgres, přidat rozšíření PostGIS do současně používané databáze:


  

sudo -u postgres psql
\c postgres
CREATE EXTENSION postgis;

Pro ověření, že extenze byla úspěšně nainstalována, můžeme skrze psql spustit příkaz \dx, který vypíše všechny nainstalované extenze v aktuální databázi.

Je vhodné nastavit heslo pro uživatele postgres, aby bylo možné se k databázi připojit i z jiných aplikací (např. pgAdmin, DBeaver, ...):


  

ALTER USER postgres WITH PASSWORD 'foobar';

Pokud s PostgreSQL chcete pracovat vzdáleně a ne skrze localhost, je potřeba ještě editovat konfigurační soubory.


  

sudo nano /etc/postgresql/*/main/postgresql.conf

Najít řádek listen_addresses a změnit na:listen_addresses = '*'

Pak je potřeba povolit připojení z jiných IP adres než localhost. To se dělá v souboru pg_hba.conf(viz řádky dole):


  

sudo nano /etc/postgresql/*/main/pg_hba.conf

Po změně bude potřeba restartovat PostgreSQL server:


  

sudo systemctl restart postgresql

Vlastní nahrání dat

Během cvičení využívám data z Natural Earth. Konkrétně Admin 0 - Countries a Populated Places. Z nějakého důvodu nelze dát přímý link ke stažení, tak se musíte proklikat skrze web.

Stažené archivy dočasně rozbalte do libovolného adresáře a spusťte načítací příkazy:


  

shp2pgsql -I -s 4326 ne_110m_admin_0_countries.shp public.countries | psql -h localhost -U postgres -d postgres
shp2pgsql -I -s 4326 ne_110m_populated_places.shp public.cities | psql -h localhost -U postgres -d postgres

Příkaz shp2pgsql slouží k převodu shapefile do SQL příkazu, který se následně spustí v databázi. Parametr -I znamená, že se vytvoří i index pro geometrii a -s 4326 určuje SRID (Spatial Reference Identifier), což je systém pro identifikaci prostorových referencí. V našem případě je to WGS 84, což je standardní systém pro GPS souřadnice.

Dotazování

V rámci dotazování budeme používat standardní SQL příkazy, které jsou rozšířeny o geospaciální funkce. Tyto funkce umožňují provádět různé prostorové operace, jako je výpočet vzdálenosti, překryv, sousedství atd.

Prozkoumejte nejprve strukturu tabulek cities a countries a najděte, jak je v nich zapsaná geometrie.

Výpis států a jejich geometrie


  

SELECT name, iso_a3, geom
FROM countries
LIMIT 10;

Při vypsání dat z tabulky countries uvidíte, že geometrie je uložena jako typ geometry, což je typ dat pro prostorové objekty v PostgreSQL. Při použití DBeaveru se automaticky použije funkce ST_AsText, která převede geometrii do textového formátu WKT (Well-Known Text) a zobrazí tedy čitelnější výsledek.

Nicméně z tohoto výpisu stejně není zřejmé, jak polygony vypadají a kde leží. Využijme tedy funkci ST_AsGeoJSON, která převede geometrii do formátu GeoJSON, který následně můžeme vizualizovat v mapovém nástroji (např. geojson.io).


  

SELECT name, iso_a3, ST_AsGeoJSON(geom)
FROM countries
LIMIT 10;

Q1 Odpovídají hranice rastrovému obrázku při výkresu pomocí geojson.io? Jaké jsou rozdíly a jaká je případně jejich příčina?

※ V reálných situacích s prostorovými daty může docházet k problému v situacích kdy se bod nachází na hranici polygonu, případně když má geometrie nízkou přesnost.

Ve kterém státě leží Praha?


  

SELECT c.name AS city, d.name AS country
FROM cities c
JOIN countries d
  ON ST_Within(c.geom, d.geom)
WHERE c.name = 'Prague';

Nově je použitý příkaz ST_Within, který vrací TRUE, pokud je geometrie prvního argumentu uvnitř geometrie druhého argumentu. Všimněte si, že tato funkce je použita jako podmínka pro spojení (technicky je v datech cities i textový sloupec s názvem státu, ale pro účely cvičení to ignoruji).

ST_Within je jednou z mnoha prostorových funkcí, které PostGIS poskytuje. Lze ji využít třeba i k upravné úloze ve smyslu lokalize libolného bodu popsaného souřadnicemi (např. v WGS 84) do konkrétního státu.


  

SELECT d.name AS country
FROM countries d
WHERE ST_Within(
          ST_SetSRID(ST_MakePoint(13.3525253, 49.7267239), 4326),
          d.geom
);

Zde dochází k definici konkrétního bodu pomocí funkce ST_MakePoint, která vytvoří bod na základě zadaných souřadnic. Funkce ST_SetSRID nastaví, že se jedná o souřadnice v systému WGS 84 (SRID 4326). Poté se tento bod porovná s geometrií států pomocí funkce ST_Within, která vrátí název státu, ve kterém se bod nachází.

※ Jedná se o souřadnice budovy FAV

Vzdálenost mezi městy


  

WITH prague AS (
  SELECT geom FROM cities WHERE name = 'Prague'
)
SELECT c.name, ST_Distance(c.geom::geography, p.geom::geography)/1000 AS dist_km
FROM cities c, prague p
WHERE c.name <> 'Prague'
  AND ST_DWithin(c.geom::geography, p.geom::geography, 500000)
ORDER BY dist_km
LIMIT 10;

V tomto dotazu je použita funkce ST_Distance, která vrací vzdálenost mezi dvěma geometrickými objekty. Všimněte si, že geometrie je převedena na typ geography, což umožňuje výpočet vzdálenosti na zakřiveném povrchu Země. Výsledek je poté převeden na kilometry dělením 1000.

ST_DWithin je funkce, která vrací TRUE, pokud je vzdálenost mezi dvěma geometrickými objekty menší než zadaná hodnota (v tomto případě 500 km).

Q2 Porovnejte výsledek s dotazem, který by nepoužil převod na typ geography.

※ Typ geometry je vhodný pro výpočty na menší ploše (např. v rámci města). Typ geography je přesnější při výpočtech vzdáleností a plochy ve skutečné zakřivené geometrii Země.

V předchozím příkladě srovnáváme umístění 2 bodů. Nicméně lze vypočítat třeba i vzálenost mezi bodem a polygonem. Lehkou úpravou tak vytvoříme dotaz na nalezení vzdálenosti místa od hranic států:


  

WITH place AS (
    SELECT ST_SetSRID(ST_MakePoint(13.3525253, 49.7267239), 4326)::geography AS geom
)
SELECT c.name, ST_Distance(c.geom::geography, p.geom)/1000 AS dist_km
FROM countries c, place p
ORDER BY dist_km
LIMIT 10;

Výběr států v regionu


  

SELECT name
FROM countries
WHERE ST_Intersects(
  geom,
  ST_MakeEnvelope(13, 48, 19, 52, 4326)  -- Roughly Central Europe
);

V tomto dotazu je použita funkce ST_Intersects, která vrací TRUE, pokud se geometrie prvního argumentu překrývá s geometrií druhého argumentu. Druhý argument je definován pomocí funkce ST_MakeEnvelope, která vytvoří obdélník (bounding box) na základě zadaných souřadnic. Koordináty (13, 48, 19, 52, 4326) představují levý dolní a pravý horní roh obdélníku v systému WGS 84 (SRID 4326).

Q3 Jaké státy jsou v tomto regionu? Jaký je rozdíl mezi ST_Intersects a ST_Within?

Výpočet rozlohy států


  

SELECT name, ST_Area(geom::geography)/1000000 AS area_km2
FROM countries
ORDER BY area_km2 DESC
LIMIT 10;

Funkce ST_Area vrací plochu geometrického objektu. Opět je převedena na typ geography, aby byla plocha vypočítána na zakřiveném povrchu Země. Výsledek je poté převeden na km2 z m2.

Q4 Jak se změní pořadí států při použití geometry místo geography?

Použití convex hull

Konvexní obálka (convex hull) je nejmenší konvexní polygon, který obklopuje danou množinu bodů. V PostGIS je možné ji vypočítat pomocí funkce ST_ConvexHull. Tato funkce je užitečná pro zjednodušení geometrie a pro analýzu tvaru objektu.


  

SELECT ST_AsGeoJSON(ST_ConvexHull(ST_Collect(c.geom)))
FROM cities c
JOIN countries d ON ST_Within(c.geom, d.geom)
WHERE d.iso_a3 = 'USA';

Tento dotaz spojuje tabulky cities a countries na základě toho, že města leží uvnitř států. Poté se použije funkce ST_Collect, která shromáždí všechny geometrie do jednoho objektu, a následně se použije ST_ConvexHull, aby se vypočítala konvexní obálka pro všechna města v USA. Výsledek je opět převeden na GeoJSON pro vizualizaci.

Q5 Jaký tvar má konvexní obálka pro města v USA?

Rozdíl mezi konvexní obálkou a spojením polygonů

Zkusme porovnat výstupy funkcí ST_ConvexHull a ST_Union. Funkce ST_Union slouží ke sloučení více geometrických objektů do jednoho. Na rozdíl od konvexní obálky, která vytváří nejmenší konvexní polygon obklopující objekty, ST_Union vytvoří složenou geometrii, která může mít libovolný tvar.

Pro ilustraci vyřešme úlohu v nalezení spojené geometrie států v Evropě (bez Ruska).


  

SELECT ST_AsGeoJSON(ST_Union(geom)) AS europe_union
FROM countries
WHERE region_un = 'Europe'
AND iso_a3 != 'RUS';


  

SELECT ST_AsGeoJSON(ST_ConvexHull(ST_Collect(geom))) AS europe_convex
FROM countries
WHERE region_un = 'Europe'
AND iso_a3 != 'RUS';

Q6 Co by se stalo při tvorbě konvexní obálky pro Evropu včetně Ruska?

Jaké státy zasahují do konvexní obálky Evropy?


  

WITH europe_hull AS (
    SELECT ST_ConvexHull(ST_Collect(geom)) AS geom
    FROM countries
    WHERE region_un = 'Europe'
    AND iso_a3 != 'RUS'
)
SELECT c.name, c.region_un, ST_Area(ST_Intersection(c.geom, h.geom)) AS overlap_area
FROM countries c, europe_hull h
WHERE c.region_un <> 'Europe'
  AND ST_Intersects(c.geom, h.geom)
ORDER BY overlap_area DESC
;

Závěrečný dotaz vrací státy, které zasahují do obálky Evropy. Používá se funkce ST_Intersection, která vrací plochu překryvu mezi geometrií státu a konvexní obálkou Evropy. Filtrace je pomocí ST_Intersects, která vrací TRUE, pokud se geometrie překrývají.

Shrnutí vybraných prostorových funkcí

Funkce Popis Příklad použití
ST_Within(a, b) Vrátí TRUE, pokud geometrie a leží uvnitř geometrie b Najdi, ve kterém státě leží město
ST_Intersects(a, b) Vrátí TRUE, pokud se geometrie a a b překrývají Vyber státy zasahující do bounding boxu
ST_Distance(a, b) Vrací nejkratší vzdálenost mezi dvěma geometriemi Vzdálenost mezi městy
ST_DWithin(a, b, d) Zda je vzdálenost mezi a a b menší než d Města do 100 km od bodu
ST_Buffer(geom, r) Vytvoří zónu r okolo geometrie Oblast 20 km kolem města
ST_Union(geoms) Spojí více geometrií do jedné složené Spojení hranic států Evropy
ST_Collect(geoms) Spojí více geometrií do jednoho multigeometrického objektu Příprava pro konvexní obálku
ST_ConvexHull(geom) Nejmenší konvexní obálka obklopující geometrii Konvexní tvar kolem měst v zemi
ST_Area(geom) Vrací plochu geometrie (např. stát) Výpočet rozlohy států
ST_MakePoint(x, y) Vytvoří bod z longitude, latitude Lokace dle GPS
ST_SetSRID(geom, srid) Nastaví prostorový referenční systém Zajištění, že bod je ve WGS 84
ST_AsText(geom) Převede geometrii do čitelného WKT formátu Ladění nebo výpis
ST_AsGeoJSON(geom) Export do GeoJSON pro mapovou vizualizaci Vizualizace přes geojson.io

📚 Více informací o funkcích a jejich použití najdete v dokumentaci PostGIS.