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

případně (v novejších verzích PostgreSQL)

postgis.sql

který se nachází např. v

/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.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, …).

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.

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:

• geografický – popisuje body sférickými souřadnicemi (v úhlových jednotkách) – šířka a délka (latitude a longitude)
  • nejkratší spojnice dvou bodů je oblouk na kouli
  • složité výpočty
  • pro přesnější výpočty (planeta není koule) je třeba uvažovat sféroid, což vše ještě komplikuje
• geometrický – pracuje s rovinou (kartézskými souřadnicemi)
  • nejkratší spojnice dvou bodů je rovná čára
  • hodně zjednodušuje výpočty (vzdáleností, ploch, průsečíků…)

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

vyjádřeno v obvyklém textové formátu, kde šířka a délka je udána ve stupních, minutách a desetinném čísle sekund. Znaménko jednotlivých souřadnic označující polokoule je udáno anglickou zkratkou světové strany (N = severní / S = jižní, E = východní / W = západní polokoule).

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;

Zde vidíte, že projekce je longlat, tedy prosté úhlové souřadnice, a elipsoid je daný standardem WGS 84.

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;

Zkuste odhadnout, co jednotlivé parametry projekce znamenají. Které popisují elipsoid ? Čím je ten elipsoid specifický ? Zkratka projekce merc odpovídá Mercator projection, což je cylindrické zobrazení. Posuďte, kde zobrazení vykazuje největší zkreslení, a ověřte na rastrové mapě na webu OpenStreetMap.

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

• bod
• linie, která se skládá z rovných úseček
• polygon – plocha ohraničená rovnými úsečkami

Abychom je mohli standardizovaně zapisovat, existují definice

• WKT – well-known text, textová reprezentace, vhodná pro zápis dotazů
• WKB – well-known binary, binární zápis, použitý pro ukládání do databáze

Příklady textové reprezentace objektů ve WKT jsou:

• POINT(0 0)
• LINESTRING(0 0,1 1,1 2)
• POLYGON( (0 0,4 0,4 4,0 4,0 0), (1 1, 2 1, 2 2, 1 2,1 1) )

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

a v čistě číselné podobě:

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)
);

Co vám dotaz vrátí za hodnotu ? V jakých je jednotkách ? Má hodnota nějaký smysl ? Jak ji PostGIS spočítal ? Zkuste dosáhnout stejného výsledku sami.

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)
);

Vypadá výsledek reálně ? Je přesný ? Proč (ne) ? Porovnejte s měřením vzdálenosti např. na mapy.cz.

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:

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:

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

Zobrazte si souřadnice na mapě, např. GoogleMaps akceptují souřadnice WGS 84 v desetinných číslech. Pozor na pořadí souřadnic (my máme x,y, ale GoogleMaps očekávají šířka,délka).

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á ?

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.

• nedokonalé normování velikosti pixelu (pohled pod úhlem) u fotomap
  • srovnej trať a elektrické vedení mapy.cz 2015 a mapy.cz 2006
• příklad práce s vrstvami - Geoportál města Kraslice
• historické mapy z 50tých let - Kontaminace @ cenia
• “nepovedený” postprocessing dat - srovnej mapy.cz a atlas.sk

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.