Search
Ukážeme si použití hotového GIS (geografického informačního systému) uloženého v databázi PostgreSQL. K těmto účelům slouží rozšíření PostGIS, které databázi spatially enables, tedy přidá prostorovou podporu.
Upozornění: vše v této sekci je už v naší databázi zařízeno, uvedené kroky proto neprovádějte. Kromě toho některé operace vyžadují administrátorská práva na databázi. Tato sekce je pro ilustraci, jak byste si takovou databázi sami zakládali.
Zbytek této sekce předpokládá, že je rozšíření PostGIS již nainstalováno v systému (ukázky jsou na linuxu), ať už ze zdrojových kódu, nebo z repository.
Pro založení databáze s GIS daty je potřeba mít v této databázi nainstalovaný SQL procedurální jazyk PL/pgSQL, který PostGIS používá pro funkce a triggery.
CREATE LANGUAGE plpgsql;
Dále je potřeba načíst definice potřebných funkci a objektů, které jsou v souboru
lwpostgis.sql
postgis.sql
/usr/share/postgresql-8.3-postgis/
Dále je šikovné si nahrát definici mnoha souřadnicových systémů do tabulky spatial_ref_sys, které jsou v souboru
spatial_ref_sys
spatial_ref_sys.sql
A ještě se hodí podpora pro integer arrays v souboru:
/usr/share/postgresql/8.3/contrib/_int.sql
Při přidávání GIS tabulek je potřeba jejich geometrické sloupce zaregistrovat v tabulce geometry_columns, buď pomocí funkce AddGeometryColumn, nebo ručně. Zde si PostGIS udržuje metadata o těchto sloupcích (typ objektů, SRID, …).
geometry_columns
AddGeometryColumn
OpenStreetMap je veřejně přístupná otevřená editovatelná mapa světa. Do této mapy může každý přispívat, ale velkým zdrojem jsou dostupná data od různých organizací (vládních, apod.), která ale musí splňovat licenci OpenStreetMap. Detaily viz OpenStreetMap wiki. Zkuste si prohlédnout vygenerované rastrové mapy. Uvidíte, že jsou na velmi dobré úrovni i v ČR.
Rastrová podoba mapy je generovaná z vektorových dat uložených právě v PostGISu na PostgreSQL. Tím, že je OpenStreetMap otevřená, můžeme si vektorová data map stáhnout. To jiné mapy (např. GoogleMaps, mapy.cz, apod.) často nenabízejí. Data jsou dostupná ve vlastním formátu – XML nebo binárním, ke kterým OpenStreetMap nabízí utilitu osm2pgsql, kterou se provede import do PostGIS databáze. Data jsou k dispozici např. na této adrese. Pro představu, vektorová data zabalená do BZIP2 mají cca 27GB. Po importu do databáze mají stovky GB. Pro ukázku jsme si do společné databáze naimportovali jen výňatek – Českou republiku, dostupnou na Geofabrik.de. Ta má zabalená cca 300MB, po importu do databáze řádově 4GB.
osm2pgsql
V GIS pracujeme s body na povrchu planety, která má tvar přibližně jako elipsoid. Abychom mohli pracovat s dvourozměrnými souřadnicemi bodů, je třeba body z povrchu planety popsat sférickými souřadnicemi (pak použijeme jen úhlové souřadnice – předpokládáme, že body leží na povrchu), nebo promítnout na nějakou plochu popsatelnou 2D souřadnicemi.
PostGIS tak umožňuje dva režimy při práci s body:
OpenStreetMap používá geometrický režim, a také ve zbytku tutoriálu se budeme zabývat výhradně jím.
Obvyklým souřadnicovým systémem je WGS 84, se kterým se můžete často setkávat např. v systémech GPS, nebo na mapách a mapových serverer. Souřadnice WGS 84 budovy E FELu na Karlově námestí je např. (zobrazte si ji na oblíbeném mapovém serveru)
50°4'36.289"N, 14°25'5.136"E
Jinak lze tytéž souřadnice vyjádřit prostou dvojicí desetinných čísel x a y (budeme dál používat v ukázkách):
14.4180933, 50.0767469
Pozor: souřadnice “šířka a délka” (“latitude a longitude”) mají opačné pořadí než “x a y” – šířka = y, délka = x.
Pro geometrický režim musíme vždy říci, v jakém promítání pracujeme – jak body z povrchu planety promítáme do roviny. Způsobů je mnoho, každý je vhodný na něco jiného. Většina používaných je standardizována pod označením EPSG, např. zmíněné WGS 84 má označení EPSG:4326. Toto a mnoho dalších najdete v tabulce spatial_ref_sys. WGS 84 tedy ale vlastně není projekce do roviny, pouze se vezmou úhlové sférické souřadnice jako souřadnice 2D prostoru. Takto s nimi tedy nelze např. počítat vzdálenosti.
PostGIS pro promítání používá projekční knihovnu PROJ.4. Zobrazte si parametry zobrazení pro WGS 84 následujícím dotazem. Uvidíte zde inicializační řetězec pro nastavení knihovny PROJ.4.
SELECT * FROM spatial_ref_sys WHERE srid=4326;
longlat
OpenStreetMap ale používá projekci Google Mercator projection, která defaultně v PostGISu není – je třeba ji tam doplnit (skript na toto je součástí utility osm2pgsql). Ta nese označení 900913, lze ji tedy vypsat následujícím dotazem:
SELECT * FROM spatial_ref_sys WHERE srid=900913;
Na okraj pro zájemce: pro ČR se často používá tzv. Křovákovo zobrazení, neboli S-JTSK, nesoucí označení EPSG:2065. Jedná se o kuželové zobrazení, které vzniklo pro českou a slovenskou kartografii, protože na území bývalého Československa vykazuje malé zkreslení – kužel je nastaven tak, že na tomto území “dobře doléhá”.
Základní objekty v GIS systémech jsou obvykle
Abychom je mohli standardizovaně zapisovat, existují definice
Příklady textové reprezentace objektů ve WKT jsou:
Ještě může být před slovem POINT/LINESTRING/POLYGON předpona MULTI, čímž se objekt daného typu může vícekrát v “zapouzdřujícím” objektu opakovat.
Pro převod pak slouží funkce:
bytea WKB = ST_AsBinary(geometry); text WKT = ST_AsText(geometry); geometry = ST_GeomFromWKB(bytea WKB, SRID); geometry = ST_GeometryFromText(text WKT, SRID);
Nejprve si zkusíme vypočítat vzdálenost mezi dvěma body. Jako jeden použijeme výše uvedený, reprezentující polohu bodovy E na Karlově náměstí, jako druhý souřadnice vchodu do bodovy FELu v Dejvicích, které jsou v textové podobě
50°6'9.105"N, 14°23'34.029"E
14.3927858, 50.10252917
Nejprve učiníme naivní pokus: zavoláme funkci na výpočet vzdálenosti na tyto souřadnice:
SELECT ST_Distance( ST_GeometryFromText('POINT(14.4180933 50.0767469)', 4326), ST_GeometryFromText('POINT(14.3927858 50.10252917)', 4326) );
Zkusíme tedy body nejprve promítnout zvoleným zobrazením, např. tím v OpenStreetMap – 900913, a pak teprve určit vzdálenost.
SELECT ST_Distance( ST_Transform(ST_GeometryFromText('POINT(14.4180933 50.0767469)', 4326), 900913), ST_Transform(ST_GeometryFromText('POINT(14.3927858 50.10252917)', 4326), 900913) );
Zkuste totéž, ale v promítání vhodném pro ČR (S-JTSK):
SELECT ST_Distance( ST_Transform(ST_GeometryFromText('POINT(14.4180933 50.0767469)', 4326), 2065), ST_Transform(ST_GeometryFromText('POINT(14.3927858 50.10252917)', 4326), 2065) );
Porovnejte s hodnotou, kterou vypočítá použití geografického režimu – výsledek je 3391.7738m. Geografický režim bohužel není ve starší verzi PostGISu, která je nainstalována na našem serveru, ještě implementován.
Nejprve si zkusíme prohlédnout záznam o nějakém bodu (v OpenStreetMap jsou samozřejmě i liniové a plošné objekty), které jsou mimo jiné v tabulce planet_osm_point, např. zvolený obchod na Karlově náměstí, který najdeme podle identifikátoru:
planet_osm_point
SELECT * FROM planet_osm_point WHERE osm_id=290036012
Podívejte se, co všechno lze v tabulce o objektu nalézt. Podívejte se na binární reprezentaci geometrie ve sloupci way. Tu si můžeme zobrazit ve WKT formátu:
way
SELECT ST_AsText(way) FROM planet_osm_point WHERE osm_id=290036012
Vidíme, že se skutečne jedná o bod. Můžeme si ho převést do WGS 84:
SELECT ST_AsText(ST_Transform(way, 4326)) FROM planet_osm_point WHERE osm_id=290036012
Nyní zkusíme hledání bodů poblíž referenčního bodu podle vzdálenosti. Např. chceme najít body nacházející se do 200m (na to se hodí predikát ST_DWithin) od budovy E FELu na Karlově náměstí. Následující příkaz najde tyto body, seřadí je podle nalezené vzdálenosti, a tyto vzdálenosti s dalšími detaily o bodech zobrazí. Rozeberte si následující příkaz, najděte parametr maximální vzdálenosti, a příkaz si vyzkoušejte. Co všechno se v tomto okolí nalézá ?
ST_DWithin
SELECT osm_id,"addr:housenumber",amenity,building,name,railway,shop,tourism,way, ST_Distance(way, ST_Transform(ST_GeometryFromText('POINT(14.4180933 50.0767469)', 4326), 900913)) FROM planet_osm_point WHERE ST_DWithin(way, ST_Transform(ST_GeometryFromText('POINT(14.4180933 50.0767469)', 4326), 900913), 200) ORDER BY ST_Distance(way, ST_Transform(ST_GeometryFromText('POINT(14.4180933 50.0767469)', 4326), 900913));
V tabulce je více jak milion bodů. Jaktože je hledání tak rychlé ? Odpovědí jsou opět indexy, tentokrát prostorové – v tomto případě je prostorový index na sloupci way. Která část příkazu index využije ?
Další příklad ukazuje, jak najdeme např. 10 nejbližších bodů. Rozdíl oproti předchozímu je, že nepoužijeme omezení na vzdálenost predikátem ST_DWithin, ale jen si je opět seřadíme podle vzdálenosti, a omezíme počet řádků ve výsledku:
SELECT osm_id,"addr:housenumber",amenity,building,name,railway,shop,tourism,way, ST_Distance(way, ST_Transform(ST_GeometryFromText('POINT(14.4180933 50.0767469)', 4326), 900913)) FROM planet_osm_point ORDER BY ST_Distance(way, ST_Transform(ST_GeometryFromText('POINT(14.4180933 50.0767469)', 4326), 900913)) LIMIT 10;
Jak dlouho trvá tento dotaz ? Porovnejte s předchozím časem, a najděte příčinu. Lze vzdálenosti indexovat ?
PostGIS používá R-stromy (ještě s pomocí GiST – Generalized Search Trees). R-stromy jsou rozšířením B+ stromů do 2D. Dělí prostor na MBR – minimal bounding rectangles, které nemusí být disjunktní. Nové body se umísťují do listu, který vyžaduje nejmenší zvětšení MBR. Ukázka na cvičení na tabuli.
Pro ty, kterým zbyde čas: můžete si zkusit vektorová data graficky prohlédnout např. programem uDig. Zde zvolíte novou vrstvu, PostGIS, a vyplníte údaje pro připojení k databázi.
Cílem cvičení je seznámení se s prací v datovém skladu pomocí aplikace JPivot.
Dokumentace k Mondrian OLAP
Opět se pracuje s daty FoodMart, jejichž schéma je zde.