====== GIS ukázka ======

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í [[http://postgis.refractions.net/|PostGIS]], které databázi //spatially enables//, tedy přidá prostorovou podporu.



==== Založení databáze pro GIS ====

**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.
<code sql>
CREATE LANGUAGE plpgsql;
</code>

Dále je potřeba načíst definice potřebných funkci a objektů, které jsou v souboru
<code>
lwpostgis.sql
</code>
případně (v novejších verzích PostgreSQL)
<code>
postgis.sql
</code>
který se nachází např. v 
<code>
/usr/share/postgresql-8.3-postgis/
</code>

Dále je šikovné si nahrát definici mnoha souřadnicových systémů do tabulky ''spatial_ref_sys'', které jsou v souboru
<code>
spatial_ref_sys.sql
</code>

A ještě se hodí podpora pro //integer arrays// v souboru:
<code>
/usr/share/postgresql/8.3/contrib/_int.sql
</code>

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 ====

[[http://www.openstreetmap.org/|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 [[http://wiki.openstreetmap.org/wiki/Main_Page|OpenStreetMap wiki]]. Zkuste si prohlédnout vygenerované rastrové [[http://www.openstreetmap.org/|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 [[http://planet.openstreetmap.org/|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 [[http://download.geofabrik.de/osm/europe/|Geofabrik.de]]. Ta má zabalená cca 300MB, po importu do databáze řádově 4GB.

==== Projekce ====

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)
<code>
50°4'36.289"N, 14°25'5.136"E
</code>
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):
<code>
14.4180933, 50.0767469
</code>

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.
<code sql>
SELECT * FROM spatial_ref_sys WHERE srid=4326;
</code>
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:
<code sql>
SELECT * FROM spatial_ref_sys WHERE srid=900913;
</code>
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á".

==== GIS objekty ====

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:
<code>
bytea WKB = ST_AsBinary(geometry);
text WKT = ST_AsText(geometry);
geometry = ST_GeomFromWKB(bytea WKB, SRID);
geometry = ST_GeometryFromText(text WKT, SRID);
</code>

==== Dotazování ====

=== Výpočet vzdálenosti ===

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ě
<code>
50°6'9.105"N, 14°23'34.029"E
</code>
a v čistě číselné podobě:
<code>
14.3927858, 50.10252917
</code>

Nejprve učiníme naivní pokus: zavoláme funkci na výpočet vzdálenosti na tyto souřadnice:
<code sql>
SELECT ST_Distance(
  ST_GeometryFromText('POINT(14.4180933 50.0767469)', 4326),
  ST_GeometryFromText('POINT(14.3927858 50.10252917)', 4326)
);
</code>
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.
<code sql>
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)
);
</code>
Vypadá výsledek reálně ? Je přesný ? Proč (ne) ? Porovnejte s měřením vzdálenosti např. na [[http://www.mapy.cz/|mapy.cz]].

Zkuste totéž, ale v promítání vhodném pro ČR (S-JTSK):
<code sql>
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)
);
</code>

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.

=== Hledání bodů ===

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:
<code sql>
SELECT *
FROM planet_osm_point
WHERE osm_id=290036012
</code>

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:

<code sql>
SELECT ST_AsText(way)
FROM planet_osm_point
WHERE osm_id=290036012
</code>

Vidíme, že se skutečne jedná o bod. Můžeme si ho převést do WGS 84:
<code sql>
SELECT ST_AsText(ST_Transform(way, 4326))
FROM planet_osm_point
WHERE osm_id=290036012
</code>
Zobrazte si souřadnice na mapě, např. [[http://maps.google.com|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á ?

<code sql>
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));
</code>

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:

<code sql>
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;
</code>

Jak dlouho trvá tento dotaz ? Porovnejte s předchozím časem, a najděte příčinu. Lze vzdálenosti indexovat ?

==== Prostorové indexy ====

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.

==== Vizualizace ====

Pro ty, kterým zbyde čas: můžete si zkusit vektorová data graficky prohlédnout např. programem [[http://udig.refractions.net/download/|uDig]]. Zde zvolíte novou vrstvu, PostGIS, a vyplníte údaje pro připojení k databázi.

===== Odkazy na mapové aplikace =====
  * nedokonalé normování velikosti pixelu (pohled pod úhlem) u fotomap 
    * srovnej trať a elektrické vedení [[https://mapy.cz/letecka?x=14.3849637&y=49.9387105&z=20&l=0| mapy.cz 2015 ]] a [[https://mapy.cz/letecka-2006?x=14.3849637&y=49.9387105&z=19&l=0|mapy.cz 2006]]
  * příklad práce s vrstvami - [[http://www.geodata.cz/geoportal/kraslice/| Geoportál města Kraslice ]]
  * historické mapy z 50tých let - [[http://kontaminace.cenia.cz/| Kontaminace @ cenia ]]
  * "nepovedený" postprocessing dat - srovnej [[https://mapy.cz/letecka?x=19.4834282&y=49.0977578&z=18&l=0|mapy.cz ]] a [[http://mapy.atlas.sk/mapa/vlachy#{%22loc%22:%22vlachy%22,%22z%22:11,%22lat%22:49.09761460783054,%22lon%22:19.483924303476744,%22src%22:%22ortoGroup%22}|atlas.sk]]

====== Datový sklad ======

Cílem cvičení je seznámení se s prací v [[http://krizik.felk.cvut.cz:9980/mondrian/testpage.jsp|datovém skladu pomocí aplikace JPivot]].

Dokumentace k [[http://mondrian.pentaho.com/documentation/schema.php|Mondrian OLAP]]

Opět se pracuje s daty FoodMart, jejichž schéma je [[https://cw.felk.cvut.cz/lib/exe/fetch.php/courses/a4b33ds/foodmarter.png|zde]].