FI MUNI, Milan Drášil 2019 Úvod do GIS 1 Fakulta Informatiky, Masarykova Universita v Brně Učební text k přednášce PV019 – Úvod do geoinformačních systémů FI MUNI, Milan Drášil 2019 Úvod do GIS 2 1. Geoinformační systém, místo na povrchu Země. 1.1. VYMEZENÍ POJMU GIS Začneme příkladem z poměrně nedávné historie. Příklad – Informační systém o nemovitostech Katastrální úřady evidují vlastnické vztahy fyzických a právnických osob k nemovitostem (budovám a pozemkům). Tato evidence existuje v různých formách v našich zemích již staletí. K prostorovému umístění evidovaných objektů sloužila mapa, k evidenci vlastnických vztahů kartotéka (zemské desky, pozemkové knihy, listy vlastnictví). FI MUNI, Milan Drášil 2019 Úvod do GIS 3 Uvažme informační systém o nemovitostech s datovým modelem v relačním databázovém systému, entity systému jsou navrženy v E-R diagramu: „Klasický“ IS je schopen reagovat na dotazy typu: - Kdo vlastní parcelu? - Jaké nemovitosti vlastní osoba? Jakou cenu mají parcely, které vlastní osoba? Lidsky čitelná odpověď na takové dotazy je převoditelná do textu/ů. GIS jsou navrhovány tak, aby byly navíc schopny reagovat na dotazy typu: - Kde se nalézá parcela? - Jaké typy parcel se nalézají v daném regionu? V tomto případě už lidsky čitelná (rozumná) odpověď/dotaz smysluplně převoditelná do textového tvaru není. Představme si, že si ve webovém prohlížeči zadáme do vyhledání „Restaurace Na Růžku, Brno“ a systém nám odpoví jen textovým polem 49.2104877N 16.6002348E. Pokud nemáme opravdu výjimečné schopnosti, pravděpodobně důležitou akci v tomto objektu zmeškáme. Půjdeme-li dále do historie, najdeme bezpočet příkladů používání map – geoinformací. Od navigace až po politiku a vojenství. Vždy se však jednalo o činnosti spojené s prostorovým umístěním jistých objektů na Zemi. Proto si můžeme v této chvíli intuitivně vymezit pojem GIS následujícím způsobem: FI MUNI, Milan Drášil 2019 Úvod do GIS 4 GIS je jakýkoli soubor postupů užívaných k ukládání a manipulování geograficky vztažených dat. Za geograficky vztažená data budeme považovat jakákoli data, která obsahují lokalizační složku, tedy taková data, která mají vazbu na místo na zemském povrchu. Typy GIS – tradiční dělení - Land Information System (LIS), Land Related Information System (LRIS), územně orientovaný informační systém, který zahrnuje vlastnické vztahy (nemovitosti, vlastníci). - Geoinformační systém – je systém pracující s jevy, která lze lokalizovat v území, ale ne vždy je lze považovat za geografická (umístění vodovodního šoupátka, dopravní značky). - Grafický informační systém – je systém pracující s geometrickými daty, která nemá smysl lokalizovat v nějakém (jednotném) prostoru, například schéma rozvodné sítě. Toto dělení ztratilo postupem doby smysl, po geoinformačních systémech je požadována komplexní funkcionalita. V tomto kurzu budeme termín: Geoinformační systém používat obecně– nebudeme polemizovat, zda mají objekty „geografický“ charakter (řeky, pohoří, etnicita obyvatelstva, kriminalita v obcích atd.) bude nám stačit, že „jev“ je jednoznačně vztažený k povrchu Země. FI MUNI, Milan Drášil 2019 Úvod do GIS 5 1.2. MÍSTO NA ZEMSKÉM POVRCHU Uveďme si pro začátek základní pojmy, které budeme potřebovat pro vymezení místa na Zemi. Geoid (1871 Listing) – Matematické těleso, model povrchu Země. Ekvipotenciální plocha vůči zemské gravitaci, která se nejvíce přimyká ke střední klidové hladině oceánu. Geometrické místo bodů se stejnou úrovní tíhového potenciálu způsobeného gravitací Země. Rovina místního poledníku – je rovina určená osou rotace Země a určovaným bodem. Zeměpisná délka (𝜆) - úhel, který svírá rovina místního poledníku procházejícího určovaným bodem a rovina nultého poledníku. Udává se většinou v úhlových jednotkách stupňů [−180°, 180°]. Zeměpisná šířka (𝜑) – úhel, který svírá rovina rovníku a přímka procházející středem Země a určovaným bodem. Udává se většinou v úhlových jednotkách [−90°, 90°]. Souřadnice [𝜑, 𝜆] jednoznačně určují místo na zemském povrchu. Loxodroma – křivka, která protíná poledníky pod stejným úhlem. Ortodroma – nejkratší spojnice dvou bodů, v případě kulového modelu Země kratší oblouk hlavní kružnice. 1.3. KDE JSEM? STRUČNÝ POHLED DO HISTORIE Eratosthenes (cca. 240 před n. l.) Pravděpodobně první dochovaný pokus o změření zemského obvodu s popisem metody pochází od Eratosthena z Kyrény. Použil úhlovou metodu, která spočívá ve změření maximální výšky slunce na dvou místech na stejném poledníku ve stejný čas. Vybral města Alexandrie a Syéné (dnešní Asuán). Eratosthenes totiž pobýval v Alexadrijské knihovně/univerzitě a stačily mu „dostupné“ prameny. FI MUNI, Milan Drášil 2019 Úvod do GIS 6 Alexandrie a Asuán – (maps.google.com) Poměr rozdílu úhlů dopadajícího slunečního paprsku (v poledne) a velikosti kruhového oblouku mezi dvěma místy na stejném poledníku = Poměr 2 a obvodu Země Vzdálenost Alexandrie a Syéné byla odhadnuta na 5000 stádií (záznamy obchodníků, pochodující legie kolem řeky Nil), úhlový rozdíl slunečních paprsků by změřen na cca. 2/50 Úhel byl měřen za slunovratu v Alexandrii, úhel slunce v poledne za slunovratu v Syéné byl znám jako /2. Sloupy kolmé k zemi nevrhaly stín – místo leží přibližně na obratníku Raka. Obvod Země = 50 x 5000 = 250 000 stádií 1 stadion = 177.6 metrů, obvod Země je tedy 44 400 000 m. FI MUNI, Milan Drášil 2019 Úvod do GIS 7 Geometrický princip Eratosthenovy metody  - maximální výška slunce za slunovratu v Alexandrii  - maximální výška slunce za slunovratu v Syéné (= /2)  - středový úhel mezi Syéné a Alexandrií  =  −  Jsou známy různé rozměry antického stadionu, od 157 do 185 m a není jasné se kterým údajem Eratosthenes pracoval. I přes cca chybu 10% se jedná o jeden z nejvýznamnějších výsledků antické vědy. Obvod Země byl potom v 9. století zpřesněn z příkazu chalífy Al-Mámuna (arabské civilizaci té doby do značné míry vděčíme za kontinuitu antické vědy). Pro určení míst bylo použito pozorování hvězd a výsledek byl neobyčejně přesný (chyba do 5%). Ptolemaios (90–160 n. l.) – první mapové dílo, zavedl do všeobecného používání úhlové jednotky stupně/minuty/vteřiny. Tabulkově zapsal zaměřené a odhadnuté pozice asi 8000 míst na Zemi (města, ústí řek, mysy, hory). Tato místa zobrazil na souboru map, zobrazujících téměř celý, do té doby poznaný svět. Geografie je v Ptolemaiově pojetí „rovinným zobrazením známé Země, se vším, co se na ní nachází“. FI MUNI, Milan Drášil 2019 Úvod do GIS 8 Gerardus Mercator – Geert de Kremer (1512-1594) – Vlámský matematik a kartograf. Je autorem, v té době revolučního, Mercatorova zobrazení: 1) Poledníky jsou přímky rovnoběžné s osou y, x-ová souřadnice je přímo úměrná zeměpisné délce. 2) Loxodromy zachovávají azimut, úhel přímky na mapě vůči ose 𝑌 se rovná azimutu. Z 1) okamžitě vyplývá: 3) Loxodromy jsou přímky. Zobrazení je velmi užitečné pro navigaci z místa A do místa B. Stačí je spojit na mapě přímkou a odečíst azimut úhloměrem. Mercator autorem kolekce map (první použil termín atlas). První mapy byly zpřesněné Ptolemaiovy mapy, postupně přidával i vlastní kartografická díla. V jeho díle pokračoval i jeho syn Rumold. Zobrazení Země v Mercatorově projekci (zdroj: http://en.wikipedia.org/wiki/Mercator_projection) Mercatorova projekce je dána rovnicemi: 𝑥 = 𝜆𝑅 𝑦 = 𝑦(𝜙)𝑅 Kde 𝑅 je poloměr Země (glóbu, záleží na měřítku...) Mercator druhou rovnici (pro 𝑦) v 16. století neznal, chyběl mu mocný prostředek, infinitezimální počet, který nám o století později poskytl Isaac Newton a Gottfried Wilhelm Leibniz. Podle pravidel 1), 2) a znalosti elementární trigonometrie byl však FI MUNI, Milan Drášil 2019 Úvod do GIS 9 schopen graficky vynést rozumnou aproximaci pravoúhlé sítě poledníků a rovnoběžek, do které vynášel polární souřadnice objektů. Mercator 1569 - Nova et Aucta Orbis Terrae Descriptio ad Usum Navigantium Emendate Accommodata Odvoďme tedy rovnici 𝑦 = 𝑦(𝜙) pro 𝑅 = 1 (měřítkový faktor můžeme doplnit kdykoli). Pomůže nám následující obrázek. Body A, B, C, D z polární projekce budou transformovány do bodů (roviny) A, B, C, D. Předpokládejme navíc, že zobrazované body jsou od sebe velmi málo („infinitezimálně“) vzdáleny, lokální kulovou plochu můžeme považovat za rovinu. - Délky oblouku 𝑑𝜙 na poledníku se v závislosti na poloze [𝜙, 𝜆] na kouli nemění, pro malé hodnoty můžeme jejich průmět do tečné roviny aproximovat přímo hodnotou 𝑑𝜙. FI MUNI, Milan Drášil 2019 Úvod do GIS 10 - Délky oblouku 𝑑𝜆 rovnoběžce jsou závislé na šířce 𝜙. Pro malé hodnoty je můžeme aproximovat hodnotou 𝑐𝑜𝑠(𝜙) 𝑑𝜆 (𝑐𝑜𝑠(𝜙) je poloměr rovnoběžky a šířce 𝜙). - Můžeme tedy položit: 𝑡𝑎𝑛(𝛼) = 𝑑𝜆𝑐𝑜𝑠(𝜙) 𝑑𝜙 𝑡𝑎𝑛(𝛽) = 𝑑𝑥 𝑑𝑦 - První rovnici upravíme do tvaru: 𝑑𝜆 = 𝑡𝑎𝑛(𝛼)𝑑𝜙 𝑐𝑜𝑠(𝜙) - Z definice projekce 𝑥 = 𝜆 okamžitě plyne 𝑑𝑥 = 𝑑𝜆, dosazením dostáváme: 𝑡𝑎𝑛(𝛽) = 𝑡𝑎𝑛(𝛼)𝑑𝜙 𝑐𝑜𝑠(𝜙)𝑑𝑦 - Požadavek “Loxodromy zachovávají azimut” lze vyjádřit rovností 𝛼 = 𝛽, dostaneme tedy jednoduchou diferenciální rovnici: 𝑦′ = 𝑑𝑦 𝑑𝜙 = 1 𝑐𝑜𝑠(𝜙) Tato rovnice není překvapivá. Uvědomíme-li si, že 2𝜋𝑐𝑜𝑠(𝜙) je délka rovnoběžky na šířce 𝜙, můžeme rovnici interpretovat: „Změna vertikálních vzdáleností rovnoběžek v Mercatorově projekci je nepřímo úměrná jejich délce“. Rovnice má bohatou historii. V roce 1599 Edward Wright sestavil tabulky pro hodnoty hledané funkce pomocí tehdejších numerických metod. Učitel navigace Henry Bond, ve 40-tých letech 17. století, na základě této tabulky, zformuloval domněnku: 𝑦 = 𝑙𝑛(𝑡𝑎𝑛 ( 𝜋 4 + 𝜙 2 )) Problémem se zabýval v roce 1665 i Isaac Newton. Zdá se, že v polovině 17. století byla tato domněnka, pro její roli v mapování a navigaci, skutečným „matematickým“ hitem. Její důkaz, tj. integraci ∫ 𝑠𝑒𝑐(𝑥)𝑑𝑥, podal v roce 1668 skotský matematik a astronom James Gregory. pramen: https://en.wikipedia.org/wiki/Mercator_projection a https://en.wikipedia.org/wiki/Integral_of_the_secant_function FI MUNI, Milan Drášil 2019 Úvod do GIS 11 1.4. VÝVOJ NAVIGAČNÍCH POMŮCEK Základní pojmy Substelární bod – je bod na zemském povrchu, ze kterého se objekt (hvězda) jeví v zenitu. Pro Slunce resp. Měsíc se také používají termíny subsolární bod resp. sublunární bod. Substelární body se pro nebeské objekty mění v čase. Greenwich Hour Angle, 𝐺𝐻𝐴 - je úhlová vzdálenost nebeského tělesa od nultého poledníku podél nebeského rovníku vztažená k času (GMT). Deklinace, 𝐷𝑒𝑐 - úhlová vzdálenost objektu od roviny světového rovníku. Dvě souřadnice [𝐷𝑒𝑐, 360° − 𝐺𝐻𝐴] jsou polární souřadnice [𝜙, 𝜆] substelárního bodu nebeského tělesa v daném čase. Z pohledu těchto lze Eratosthenovu metodu změření velikosti Země popsat takto: - Za slunovratu je v poledne subsolární bod v Syéné. - Syéné je (+/-) na stejném poledníku jako Alexandrie. Ze znalosti vzdálenosti Alexandrie – Syéné a úhlu Slunce nad obzorem mohu tedy dopočítat obvod Země. Poloha bodového objektu na Zemi je určena úhlovými souřadnicemi [𝜙, 𝜆]. Proto není překvapivé, že historické navigační pomůcky jsou úhloměrná zařízení. Úhloměry Kamal – dřevěná destička, v jejímž středu je upevněn provázek. Na provázku jsou uzly, které reprezentují zeměpisnou šířku. Provázek se chytne do zubů v místě uzlu s cílovou šířkou. Spodní hrana destičky se ztotožní s horizontem. Je-li Polárka pod horní hranou, jsme jižněji, naopak severněji. Splývá-li Polárka s horní hranou, jsme na správné šířce. Kamal (obrázek z http://www.lode.cz/re.php) Jakubova hůl – Posuvná destička na tyči, na které je vyznačena stupnice. Její dolní okraj ztotožníme posouváním s horizontem a horní s měřeným objektem, oko je na konci hole. Na tyči potom odečteme údaj. Byla používána od starověku až do 19. století. FI MUNI, Milan Drášil 2019 Úvod do GIS 12 Jakubova hůl, John Sellers - Practical Navigation (1672) Námořní astroláb – závěsná kruhová deska s otočným ramenem pro odměření úhlu proti horizontu. Námořní astroláb (obrázek z http://cs.wikipedia.org) Sextant – optický přístroj využívající polopropustné zrcadlo ke ztotožnění měřeného objektu s horizontem. Princip byl objeven nezávisle Isaacem Newtonem a Johnem Hadleym. Poskytoval řádově přesnější výsledky (jednotky úhlových minut), než všechny dosavadní přístroje. Poznamenejme, že sextantem lze smysluplně měřit úhel vzdálených objektů, tj. u kterých je zanedbatelná paralaxa způsobená rozdílnou polohou zrcátek. FI MUNI, Milan Drášil 2019 Úvod do GIS 13 Sextant (zdroj http://en.wikipedia.org/wiki/sextant) Určování šířky Historicky se používaly nebeské objety, u kterých se substelární bod neměnil v čase (resp. měnil málo), například Polárka. Například, změřený úhel polárky proti horizontu je (+/-) přímo šířka. Určování délky Je problém výrazně komplikovanější, metody vychází z měření místního času, resp. místního poledne. K tomu však potřebujeme „univerzální hodiny“. V historii řešení tohoto problému najdeme skutečně významné osobnosti – Gallileo Galilei, Edmont Halley, Leonhard Euler… Za zmínku také stojí, že odměna za uspokojivé vyřešení byla uzákoněna parlamentem Spojeného království v roce 1714 (Longitude Act). První zaznamenaná metoda pochází od Ptolemaia a je založena na skutečnosti, že Měsíc se pohybuje vůči hvězdám (zhruba rychlostí svého poloměru za hodinu). Stačí tedy tabulkově zaznamenat jeho pozici vůči jasným hvězdám v časech nultého poledníku a můžeme (opět úhloměrnými pomůckami), kdekoli stanovit čas nultého poledníku. Místní čas potom stanovíme okamžikem místního poledne (slunce je nejvýše nad obzorem) a zeměpisnou délku podle rozdílu místního času a času nultého poledníku (24 h = +/- 360°). Tato metoda byla značně nepřesná, proto i značné chyby v délkách v Ptolemaiových mapách. FI MUNI, Milan Drášil 2019 Úvod do GIS 14 Až do vynálezu „přesných“ hodin nebyl však problém určení uspokojivě vyřešen. Přesné hodiny (chyba v minutách za týdny provozu) byly sestrojeny Johnem Harissonem (1693–1776). Harrisonovy lodní hodiny zdroj: https://www.rmg.co.uk/discover/explore/longitude-found-john-Harrison K vyřešení problému určení zeměpisné délky výraznou měrou také přispěly výsledky Johannese Keplera a Isaaca Newtona vysvětlující pohyb Země kolem slunce. Metody určení pozice s úhloměrem a hodinami K určení pozice na Zemi tedy použijeme: 1) Hodiny ukazující přesný čas (nultého poledníku). 2) Úhloměrné zařízení (sextant). 3) Metodu pro určení substelárního bodu. Metoda místního poledne. Provedeme měření úhlu 𝛼 slunce nad horizontem a zaznamenáme čas 𝑡 (GMT) v kulminaci. Čas kulminace je roven času místního poledne. Substelární bod [𝑑𝑒𝑐, 360° − 𝐺𝐻𝐴] leží na našem poledníku. Úhlová vzdálenost k substelárnímu bodu je 𝜋 2 − 𝛼 a tedy naše šířka je: 𝜙 = 𝜋 2 − 𝛼 + 𝑑𝑒𝑐 Nabízí se použít přímo čas našeho místního poledne 𝑡 pro výpočet naší délky. Vlivem sklonu zemské osy vůči rovině ekliptiky a 2. Keplerova zákona, se však čas FI MUNI, Milan Drášil 2019 Úvod do GIS 15 místního poledne (čas kulminace) různí od středního poledne (čas na hodinkách) až o cca. 16 minut. Rozdíl času oproti nultému poledníku tedy musíme korigovat. Provedeme tedy korekci času 𝑡. Označíme-li Δ𝑡 rozdíl, mezi korigovaným časem a polednem v hodinách potom je naše délka ve stupních (za hodinu 15 stupňů): 𝜆 = 15Δ𝑡 Korekční tabulku může velmi zhruba nahradit křivka 𝑎𝑛𝑎𝑙𝑒𝑚𝑎, kde na jedné ose je zobrazena deklinace (𝑑𝑒𝑐), na druhé časová diference místního poledne proti střednímu poledni. Na křivce je potom vyznačený průběh roku. FI MUNI, Milan Drášil 2019 Úvod do GIS 16 Analema – zdroj U.S. Coast and Geodetic Survey Na první pohled jsou korekce podle analemy velmi hrubé. Nezohledňují hodiny, resp. přestupné roky a vzhledem k tomu, že korekce se mění poměrně dramaticky, v určitém období i o 30 úhlových minut za den, se pro praktickou navigaci nehodila. Používaly přesnější tabulky vydávané pro každý rok v tzv. Námořním almanachu. V tomto jsou pro jednotlivé dny rozepsány potřebné údaje po hodinách. To neznamená, že můžeme měřit jenom v celé hodiny – lineární aproximace mezi dvěma řádky v almanachu je dostačující, chyba aproximace je zanedbatelná. Námořní almanach – zdroj http://thenauticalalmanac.com Tabulky v Námořním almanachu nám zejména poskytují souřadnice substelárního (subsolárního, sublunárního) bodu planet, navigačních hvězd, Slunce a Měsíce vzhledem k času GMT a časové korekce. Zeměpisnou délku v metodě místního poledne můžeme získat také přímo ze sloupce 𝐺𝐻𝐴 – substelární bod leží na našem místním poledníku. Čas kulminace však není úplně jednoduché odhadnout. Chyba v odhadu poledne v řádu 1 minuty nám způsobí v naší zeměpisné šířce chybu polohy cca 18 km. Metoda pozičních kružnic. Opět použijeme úhloměrné zařízení, Námořní almanach a přesné hodinky. Změříme úhel navigačního objektu (hvězda, planeta, slunce) proti horizontu. Geometrické místo na zemském povrchu, ze kterého je vidět objekt pod změřeným úhlem je kružnice se středem v příslušném substelárním bodě. Její (úhlový) poloměr je 𝜋 2⁄ mínus změřený úhel. Takové kružnice získáme dvě, tři a více. Společný průsečík je naše poloha. FI MUNI, Milan Drášil 2019 Úvod do GIS 17 Rovnice pro výpočet průsečíku pozičních kružnic: Vstup: [𝜙1, 𝜆1] – pozice substelárního bodu 1. objektu (z almanachu) [𝜙2, 𝜆2] – pozice substelárního bodu 2. objektu (z almanachu) 𝛼1 – naměřený úhel 1. objektu proti horizontu 𝛼2 – naměřený úhel 2. objektu proti horizontu 𝛿1 = 𝜋 2 − 𝛼1 – úhlová vzdálenost k 1. bodu 𝛿2 = 𝜋 2 − 𝛼2 – úhlová vzdálenost k 2. bodu Převod do kartézské soustavy: Převedeme substelární body do kartézské soustavy na jednotkovou kouli se středem v nule (skutečný poloměr Země nemusíme zahrnovat, jde nám o jednotkové směrové vektory). Rovnice pro tento převod mají následující tvar: 𝐹(𝜙, 𝜆) = [𝑥, 𝑦, 𝑧] (0) kde: 𝑥 = 𝑐𝑜𝑠(𝜆) 𝑐𝑜𝑠(𝜙) 𝑦 = 𝑠𝑖𝑛(𝜆) 𝑐𝑜𝑠(𝜙) 𝑧 = 𝑠𝑖𝑛(𝜙) (pro zvídavé čtenáře: zkuste si tento převod ověřit, stačí tužka a papír, začněte od “𝑧”) Pomocí těchto rovnic převedeme oba substelární body do jednotkových vektorů [𝑎1, 𝑏1, 𝑐1] resp. [𝑎2, 𝑏2, 𝑐2]: [𝑎1, 𝑏1, 𝑐1] = 𝐹(𝜙1, 𝜆1) [𝑎2, 𝑏2, 𝑐2] = 𝐹(𝜙2, 𝜆2) Sestavení rovnic: Nyní můžeme sestavit rovnice, jejichž řešení nám přinese kartézskou reprezentaci hledaných průsečíků na jednotkové kouli. Položme nejprve: 𝑑1 = 𝑐𝑜𝑠 (𝛿1) 𝑑2 = 𝑐𝑜𝑠(𝛿2) a sestavme soustavu rovnic. Použijeme formule, která nám říká, že skalární součin dvou jednotkových vektorů se rovná kosinu úhlu jimi sevřených: 𝑥2 + 𝑦2 + 𝑧2 = 1 (1) 𝑎1 𝑥 + 𝑏1 𝑦 + 𝑐1 𝑧 = 𝑑1 (2) 𝑎2 𝑥 + 𝑏2 𝑦 + 𝑐2 𝑧 = 𝑑2 (3) (1) Bod leží na jednotkové kouli. (2) Úhel k 1. substelárnímu bodu je 𝛿1. (3) Úhel k 2. substelárnímu bodu je 𝛿2. Výpočet: Rovnice (2) a (3) jsou lineární, lze z nich tedy snadno (např. Cramerovým pravidlem) odvodit 𝑥, resp. 𝑦 v závislosti na 𝑧. FI MUNI, Milan Drášil 2019 Úvod do GIS 18 𝑥 = 𝑒1 𝑧 + 𝑓1 (4) 𝑦 = 𝑒2 𝑧 + 𝑓2 (5) kde: 𝑒1 = 𝑏1 𝑐2−𝑏2 𝑐1 |𝑆| , 𝑓1 = 𝑏2 𝑑1−𝑏1 𝑑2 |𝑆| 𝑒2 = 𝑎2 𝑐1−𝑎1 𝑐2 |𝑆| , 𝑓2 = 𝑎1 𝑑2−𝑎2 𝑑1 |𝑆| |𝑆| = 𝑎1 𝑏2 − 𝑎2 𝑏1 Substitucí (4) a (5) do (1) dostaneme kvadratickou rovnici 𝐴𝑧2 + 𝐵𝑧 + 𝐶 = 0 (6) kde: 𝐴 = 𝑒1 2 + 𝑒2 2 + 1, 𝐵 = 2(𝑒1 𝑓1 + 𝑒2 𝑓2), 𝐶 = 𝑓1 2 + 𝑓2 2 − 1 Řešení: Vyřešíme kvadratickou rovnici (6) a dostaneme dvě hodnoty 𝑧1,2 = −𝐵 ± √𝐵2 − 4𝐴𝐶 2𝐴 Z rovnic (4) a (5) dostaneme kartézské souřadnice pro průsečíky [𝑥1, 𝑦1, 𝑧1] a [𝑥2, 𝑦2, 𝑧2]. Pro zpětný převod do polárních souřadnic použijeme triviálně odvoditelnou inverzní funkci: 𝐹−1(𝑥, 𝑦, 𝑧) = [𝜙, 𝜆] 𝜙 = 𝑠𝑖𝑛−1(𝑧) 𝜆 = 𝑐𝑜𝑠−1 ( 𝑥 𝑐𝑜𝑠 (𝜙) ) (ponechávám čtenáři vyřešení případu, kdy se nacházíme na pólech a tedy 𝑐𝑜𝑠(𝜙) = 0). Přímý výpočet je pro praktickou navigaci příliš komplikovaný, představme si, že nemáme po ruce žádný „počítač“ nebo kalkulačku, který/á za nás uvedený postup provede, máme jen tabulky, tužku a papír. Dále, není možné vynášet do map kružnice „graficky“. Jednak jsou jejich poloměry hodně velké (tisíce km, tedy i v měřítku 1:1000 000 bychom museli mít několikametrové kružítko a mapu), jednak obrazem kružnice na zemském povrchu ve 2D mapě není kružnice. FI MUNI, Milan Drášil 2019 Úvod do GIS 19 Proto se používala metoda pozičních linií (interceptu), ve které jsou kružnice nahrazeny tečnami – ty už lze do mapy vynášet pravítkem. Postup je následující: - Změříme úhel nebeského objektu a zaznamenáme čas měření a z almanachu určíme substelární bod SP. - Odhadneme naši polohu bodem AP (Assumed Position). - Z almanachu odečteme souřadnice substelárního bodu SP v zaznamenaném čase. - Určíme azimut AP-SP, kolmice němu je (+/-) rovnoběžná s tečnou poziční kružnice. - Úhlová délka naší skutečné pozice k SP je doplněk do 90° námi změřeného úhlu, úhlovou délku AP-SP lze spočítat. Tečnu z předešlého bodu posuneme ve směru azimutu podle diference úhlových délek. - Postup opakujeme pro jiný objekt a tím dostaneme více tečen. Jejich „průsečík“ je naše pozice. - Azimut, resp. úhlové vzdálenosti počítáme pomocí vzorců sférické geometrie. Metoda má tu nespornou výhodu, že je použitelná v libovolném čase, ve dne můžeme měřit Slunce v různých časech. Metoda interceptu ze tří měření FI MUNI, Milan Drášil 2019 Úvod do GIS 20 Sférická geometrie – úhlová vzdálenost a azimut Úhlová vzdálenost bodů na kouli: [𝜙1, 𝜆1] – pozice 1. bodu [𝜙2, 𝜆2] – pozice 2. bodu 𝛿 – úhlová vzdálenost Metoda – spočteme kartézské souřadnice 1. a 2. bodu podle (0). Skalární součin těchto bodů (vektorů) je kosinus hledané úhlové vzdálenosti, tedy: 𝑐𝑜𝑠(𝛿) = 𝑐𝑜𝑠(𝜆1)𝑐𝑜𝑠(𝜙1)𝑐𝑜𝑠(𝜆2)𝑐𝑜𝑠(𝜙2) + 𝑠𝑖𝑛(𝜆1)𝑐𝑜𝑠(𝜙1)𝑠𝑖𝑛(𝜆2)𝑐𝑜𝑠(𝜙2) + 𝑠𝑖𝑛(𝜙1)𝑠𝑖𝑛(𝜙2) = 𝑐𝑜𝑠(𝜙1)𝑐𝑜𝑠(𝜙2)(𝑐𝑜𝑠(𝜆1)𝑐𝑜𝑠(𝜆2) + 𝑠𝑖𝑛(𝜆1)𝑠𝑖𝑛(𝜆2)) + 𝑠𝑖𝑛(𝜙1)𝑠𝑖𝑛(𝜙2) 𝛿 = 𝑐𝑜𝑠−1 (𝑐𝑜𝑠(𝜙1)𝑐𝑜𝑠(𝜙2)𝑐𝑜𝑠(𝜆2 − 𝜆1) + 𝑠𝑖𝑛(𝜙1)𝑠𝑖𝑛(𝜙2)) (7) Výchozí azimut z bodu do bodu (výchozí úhel ortodomy): Použijeme základní formuli sférické trigonometrie, větu kosinovou. Věta kosinová: Pro libovolný trojúhelník ABC na kulové ploše (tj. strany jsou ortodromy) platí: 𝑠𝑖𝑛(𝑎) 𝑠𝑖𝑛(𝛼) = 𝑠𝑖𝑛(𝑏) 𝑠𝑖𝑛(𝛽) = 𝑠𝑖𝑛(𝑐) 𝑠𝑖𝑛(𝛾) Pozn: strany 𝑎, 𝑏, 𝑐 jsou na kouli také úhly. Předpokládejme nyní, že: A – naše pozice [𝜙0, 𝜆0] B – je severní pól C – pozice substelárního bodu [𝜙1, 𝜆1] 𝛼 – výchozí azimut FI MUNI, Milan Drášil 2019 Úvod do GIS 21 Zřejmě: 𝛽 = 𝜆1 − 𝜆0 𝑏 = 𝛿 (úhlová vzdálenost AC z předešlého výpočtu) 𝑎 = 𝜋 2 − 𝜙1 𝑠𝑖𝑛(𝑎) 𝑠𝑖𝑛(𝛼) = 𝑠𝑖𝑛(𝑏) 𝑠𝑖𝑛(𝛽) tedy, 𝑠𝑖𝑛(𝛼) = 𝑠𝑖𝑛(𝑎) 𝑠𝑖𝑛(𝛽) 𝑠𝑖𝑛(𝑏) (8) 𝛼 = 𝑠𝑖𝑛−1 ( 𝑐𝑜𝑠(𝜙1) 𝑠𝑖𝑛(𝜆1 − 𝜆0) 𝑠𝑖𝑛(𝛿) ) Majáky Metody navigace byly v průběhu vývoje doplňovány pozemními navigačními body. Optické majáky pro lokální navigaci. Byly budovány od Říše římské, dosud hojně používané na pobřežích všech moří. Radiomajáky. Síť budována od druhé poloviny 20. století, hlavní navigační metoda pro leteckou dopravu. Původní nesměrový systém NDB, který vyžadoval směrovou anténu přijímače, je postupně nahrazován systémem VOR, u kterého směrová anténa nutná není. Pomocí statické a rotační antény, VOR systém vysílá směrově závislý signál, přijímač je schopen určit azimut majáku nezávisle na síle signálu. Global Positioning System GPS (Viz https://en.wikipedia.org/wiki/Global_Positioning_System) V současné době je nejmasovějším prostředkem pro určení polohy systém GPS. Je provozován Ministerstvem obrany USA, vývoj se datuje od první poloviny 70. let 20. století. Jeho využívání pro civilní účely je bohužel spjato s tragickým osudem 269 pasažérů letu KAL 007. Civilní let KAL 007, společnosti Korean Air se vlivem navigační chyby v roce 1983 odchýlil o běžné trasy a dostal se tak do sovětského vzdušného prostoru. Co způsobilo chybu není jasné, snad chybně zadaná trasa do palubního počítače. V době studené války se tato chyba stala osudnou. Civilní let byl sovětskou stíhačkou sestřelen. Po tomto hrůzném incidentu oznámil americký prezident Ronald Reagan, že po dokončení bude systém GPS bude k dispozici i pro civilní účely. GPS je založen na 32 satelitech schopných vysílat svoji přesnou polohu a přesný čas. Čas je pro všechny satelity velmi přesně synchronizován. Přijímače mají hodiny, které nejsou synchronizovány se satelity, ale jsou schopny měřit velmi přesně časový interval. FI MUNI, Milan Drášil 2019 Úvod do GIS 22 Úsměvný se z dnešního pohledu jeví počet satelitů – pro jejich identifikaci bylo rezervováno 5 bitů, přidání dalšího by znamenalo změnit komunikační protokol. Pro základní představu uvedeme analogii v 1D prostoru. Příklad – 1D GPS: Bedřich, Jan a Ctirad mají jazzovou kapelu a shodou okolností bydlí po řadě v Bedřichově Janově a Ctiradově. Tato sídla jsou na jedné rovné linii železniční trati Bedřichov–Janov–Ctiradov. Janov je třeba na půli cesty, vzdálenosti všichni dobře znají. Bedřich a Ctirad jsou dříči a pedanti, každou celou hodinu (na milisekundy!) si přehrají nutné etudy. Bedřich vždy začíná tónem B. Není divu, jeho nástroj je soprán saxofon. Ctirad hraje C-kornet, snad proto jeho cvičení začíná vždy tónem C. Zato Jan, to je bohém. Vůbec necvičí, ale má absolutní hudební sluch a geniální cit pro rytmus. A ta jeho životospráva... Jednou po zkoušce trochu přebral a ujel mu vlak. I vydal se podél trati a po cestě usnul. Po probuzení byl zoufalý, nevěděl, odkud vyšel a kterým směrem se vydat. Až najednou, Bedřich a Ctirad začali cvičit – a Jan najednou ví, kde je… Zaznamenal časový rozdíl, kdy zaslechl Bedřichův a Ctiradův nástroj. Poznámka: Uvedený příklad lze vyřešit triviální lineární rovnicí, ve více dimenzích je to obtížnější. Řešení vede na soustavu nelineárních rovnic. Zkusme si nyní tyto rovnice sestavit. Pro jednoduchost zanedbáme některé faktory, které se musí do reálného výpočtu zahrnout: - čas vyslání signálu satelitu je mírně zpožděn oproti poloze - ionosférická a troposférická refrakce elektromagnetických vln - vliv ionosféry a troposféry na rychlost světla - relativistické vlivy, satelity se pohybují nezanedbatelnou rychlostí v gravitačním poli Země Mějme tedy 4 satelity, které vyslaly svoji polohu a čas: [𝑥𝑖, 𝑦𝑖, 𝑧𝑖, 𝑡𝑖], 𝑖 = 1, . . ,4 Tyto signály zachytil přijímač v časech svých hodin: 𝜏𝑖, 𝑖 = 1, . . ,4 Hodiny přijímače nejsou přesně synchronizovány s hodinami satelitů. Předpokládejme ale, že se jejich rozdíl během měření nezměnil a označme ho symbolem 𝜏. Označme hledanou pozici: 𝑿 = [𝑥, 𝑦, 𝑧] a 𝑟𝑖 – vzdálenost hledané pozice od satelitu 𝑖. FI MUNI, Milan Drášil 2019 Úvod do GIS 23 Je 𝑟𝑖 = 𝑟𝑖(𝑿) = √(𝑥 − 𝑥𝑖)2 + (𝑦 − 𝑦𝑖)2 + (𝑧 − 𝑧𝑖)2 a 𝑟𝑖 = (𝜏𝑖 − 𝑡𝑖)𝑐 − 𝜏𝑐 kde 𝑐 je rychlost světla. Dostáváme soustavu nelineárních rovnic, 𝑖 = 1, . . ,4: √(𝑥 − 𝑥𝑖)2 + (𝑦 − 𝑦𝑖)2 + (𝑧 − 𝑧𝑖)2 + 𝜏𝑐 = (𝜏𝑖 − 𝑡𝑖)𝑐 (1) Označme levou stranu rovnic: 𝐹𝑖(𝑥, 𝑦, 𝑥, 𝜏) = √(𝑥 − 𝑥𝑖)2 + (𝑦 − 𝑦𝑖)2 + (𝑧 − 𝑧𝑖)2 + 𝜏𝑐 Soustavu rovnic: 𝐹𝑖(𝑥, 𝑦, 𝑥, 𝜏) = (𝜏𝑖 − 𝑡𝑖)𝑐 nyní linearizujeme (abychom ji vůbec mohli nějak řešit). To provedeme tak, že funkci 𝐹𝑖(𝑥, 𝑦, 𝑥, 𝜏) aproximujeme lineárním členem Taylorova rozvoje v nějakém bodě 𝑿0 = [𝑥0, 𝑦0, 𝑧0] a časovým posunem hodin přijímače od satelitů 𝜏0. Poznamenejme, že tento postup není výlučný pro výpočet polohy přijímače GPS, jedná se o aplikaci Newton-Raphsonovy metody řešení nelineárních rovnic, kdy jsou známé první parciální derivace levých stran, viz Příloha 9.1. 𝐹𝑖(𝑿, 𝜏) ≈ 𝐹𝑖(𝑿0, 𝜏0) + 𝜕𝐹𝑖 𝜕𝑥 (𝑿0, 𝜏0)(𝑥 − 𝑥0) + 𝜕𝐹𝑖 𝜕𝑦 (𝑿0, 𝜏0)(𝑦 − 𝑦0) + 𝜕𝐹𝑖 𝜕𝑧 (𝑿0, 𝜏0)(𝑧 − 𝑧0) + 𝜕𝐹𝑖 𝜕𝜏 (𝑿0, 𝜏0)(𝜏 − 𝜏0) Zřejmě (prostou parciální derivací podle 𝑥, 𝑦, 𝑧, 𝜏) je: 𝜕𝐹𝑖 𝜕𝑥 (𝑿0, 𝜏0) = 𝑥0 − 𝑥𝑖 𝑟𝑖(𝑿0) , 𝜕𝐹𝑖 𝜕𝑦 (𝑿0, 𝜏0) = 𝑦0 − 𝑦𝑖 𝑟𝑖(𝑿0) 𝜕𝐹𝑖 𝜕𝑧 (𝑿0, 𝜏0) = 𝑧0 − 𝑧𝑖 𝑟𝑖(𝑿0) , 𝜕𝐹𝑖 𝜕𝜏 (𝑿0, 𝜏0) = 𝑐 FI MUNI, Milan Drášil 2019 Úvod do GIS 24 Označíme-li Δ𝑥 = 𝑥 − 𝑥0, Δ𝑦 = 𝑦 − 𝑦0, Δz = 𝑧 − 𝑧0,Δ𝜏 = 𝜏 − 𝜏0 dostáváme soustavu lineárních rovnic: [ 𝑥0 − 𝑥1 𝑟1(𝑿0) 𝑦0 − 𝑦1 𝑟1(𝑿0) 𝑧0 − 𝑧1 𝑟1(𝑿0) 𝑐 𝑥0 − 𝑥2 𝑟2(𝑿0) 𝑦0 − 𝑦2 𝑟2(𝑿0) 𝑧0 − 𝑧2 𝑟2(𝑿0) 𝑐 𝑥0 − 𝑥3 𝑟3(𝑿0) 𝑦0 − 𝑦3 𝑟3(𝑿0) 𝑧0 − 𝑧3 𝑟3(𝑿0) 𝑐 𝑥0 − 𝑥4 𝑟4(𝑿0) 𝑦0 − 𝑦4 𝑟4(𝑿0) 𝑧0 − 𝑧4 𝑟4(𝑿0) 𝑐 ] [ Δ𝑥 Δ𝑦 Δ𝑧 Δ𝜏] = [ (𝜏1 − 𝑡1)𝑐 − 𝐹1(𝑿0, 𝜏0 ) (𝜏2 − 𝑡2)𝑐 − 𝐹2(𝑿0, 𝜏0 ) (𝜏3 − 𝑡3)𝑐 − 𝐹3(𝑿0, 𝜏0 ) (𝜏4 − 𝑡4)𝑐 − 𝐹4(𝑿0, 𝜏0 )] (2) Řešení soustavy (1) dostaneme iteračním procesem. Odhadneme počáteční hodnoty (𝑿0, 𝜏0) a řešíme systém (2). Dostaneme „opravu“ Δ(𝑿, 𝜏) = Δ𝑥, Δ𝑦, Δz, Δ𝜏 pro následující krok iterace použijeme odhad (𝑿0, 𝜏0) ∶= (𝑿0, 𝜏0) + Δ(𝑿, 𝜏). Proces opakujeme, dokud |Δ(𝑿, 𝜏)| ≈ 0. Pro počáteční odhad použijeme jakékoli místo na Zemi, nebo [0,0,0,0] metoda konverguje velmi rychle, stačí jednotky iterací a |Δ𝑿| ≤ 1 𝑚. V případě, že máme k dispozici více satelitů, je do výpočtu ještě zapojena lineární regrese (metoda nejmenších čtverců, viz Přílohy). Rovnice (2) mají více řádků, než 4 (počet viditelných satelitů), maticově: 𝑨Δ𝑥 = 𝑏 v každém kroku řešíme soustavu rovnic: 𝑨 𝑇 𝑨Δ𝑥 = 𝑨 𝑇 𝑏 (3) Funkce 𝐹𝑖 jsou nazývány „𝑝𝑠𝑒𝑢𝑑𝑜 𝑟𝑎𝑛𝑔𝑒“ – pseudovzdálenost. Náš zjednodušený příklad výpočtu polohy přijímače GPS jasně ukazuje, že bez výpočetní techniky se neobejdeme. GPS je mimo jiné i počítač. Ruční počítání rovnic (2) nebo (3) je (dnes) nepředstavitelné. Pro odvození rovnic viz např.: http://www.nbmg.unr.edu/staff/pdfs/Blewitt%20Basics%20of%20gps.pdf Měření polohy pomocí GPS dosahuje přesnost 10-20 metrů. GLONAS Ruskou variantou je systém GLONAS, původně vyvinut také pro vojenské účely, do operačního provozu byl uveden 2011. BeiDou Čínská varianta, očekávaná plná operační schopnost v roce 2020. FI MUNI, Milan Drášil 2019 Úvod do GIS 25 Galileo Projekt EU, očekávaná plná operační schopnost také v roce 2020. Pomocí diferenciálních GPS stanic a lokálních korekcí lze dosáhnout přesnost měření až 0.1 m. Princip je ten, že máme k dispozici systém stacionárních GPS stanic, u kterých je přesně zaměřená poloha. Tato GPS vysílá odchylky vypočtené polohy od své skutečné polohy a diferenciální přijímače GPS v okolí podle toho korigují svou polohu. Předpokladem je to, že v relativně malém okolí mají všechny přijímače GPS stejnou chybu. Navigátoři by však měli být schopni určit svou polohu, pro případ kolapsu technologie, čistě pomocí mechanických pomůcek a tabulek. 1.5. GEODÉZIE https://cs.wikipedia.org/wiki/Geodézie Doposud jsme se zabývali určením místa na zemském povrchu v malých měřítkách. S výjimkou diferenciálních GPS jsme nedosahovali přesnosti, která je vyžadována pro evidenci a analýzu jevů, jako jsou nemovitosti, inženýrské sítě a podobně. Určením vzájemné polohy bodů na zemském povrchu nebo v prostoru ve zvoleném souřadnicovém systému se zabývá obor geodézie. V základním dělení se geodézie dělí na vyšší a nižší. Vyšší geodézie – se zabývá územím velkého rozsahu, pro které je nutné uvažovat zakřivení zemského povrchu. Výsledkem její činnosti jsou hlavně zaměření a výpočty geodetických sítí, které tvoří základ pro podrobná měření polohopisná a výškopisná. Vyšší geodézie už nepovažuje tvar Zemského tělesa za kouli, musí počítat s přesnějším modelem – elipsoidem, viz další kapitola. Nižší geodézie – se zabývá územím, na kterém lze zakřivení zemského povrchu zanedbat. Výsledkem její činnosti je zejména podrobné měření polohopisné a výškopisné. Česká státní trigonometrická síť https://bodovapole.cuzk.cz/vyznamneTB.aspx Je síť stabilizovaných bodů, tvoří základ pro podrobná měření na území ČR. Trigonometrický bod se podle metody zaměření/výpočtu dělí na 5 kategorií Poloha bodů I. Řádu byla přesně geometricky zaměřena v rámci trigonometrické (trojúhelníkové) sítě metodou triangulace. Pomocí geodetického přístroje pro měření úhlů – teodolitu byly měřeny vodorovné úhly mezi směry na sousední body trigonometrické sítě a následně byly výpočtem určeny přesné souřadnice, které se využívají pro další geodetická měření a mapování. Po vzniku Československa bylo rozhodnuto vybudovat jednotné geodetické základy na celém území státu, proto byla vybudována Jednotná trigonometrická síť katastrální. FI MUNI, Milan Drášil 2019 Úvod do GIS 26 Na území České republiky bylo zřízeno 181 trigonometrických bodů I. řádu a přibližně 75 tisíc bodů II. – V. řádu. Česká státní trigonometrická síť I. Řádu https://bodovapole.cuzk.cz/vyznamneTB.aspx Pro podrobné měření bodů geodeti používají zejména tyto přístroje: Teodolit – Optický přístroj pro měření vodorovných a vertikálních úhlů objektů s přesností až na úhlové vteřiny. Dálkoměr – Přístroj pro měření délek. Totální stanice – Měření úhlů i délek, většinou je vybavena počítačem pro ukládání bodů a provádění výpočetních úloh. Nivelační přístroj – přístroj pro určování výšek, vytváří vodorovnou rovinu pomocí vodorovného dalekohledu. FI MUNI, Milan Drášil 2019 Úvod do GIS 27 2. Kartografická zobrazení a mapy 2.1. ZÁKLADNÍ TYPY ZOBRAZENÍ Zemský povrch, geoid, je geometricky poměrně složitý útvar, proto je modelován rotačním elipsoidem, který je určen hlavní a vedlejší poloosou. Pro různá kartografická zobrazení jsou používány různé elipsoidy. Bessel Hayford Krasovskij IAG 1967 WGS-84 rok 1841 1909 1940 1967 1984 a[m] 6377397.16 6378388 6378245 6378160 6378137 b[m] 6356078.96 6356911.95 6356863.02 6356774.52 6356752.31 Geodetické datum: model zemského tělesa (koule, elipsoid ..), jeho umístění (orientace) vůči zemskému tělesu a datum určení. Matematická kartografie – disciplína zabývající se převodem zemského povrchu do roviny. Konformní zobrazení – ponechává nezkreslené úhly, značně jsou však zde zkreslovány plochy. Ekvidistantní zobrazení – nezkresluje délky určité soustavy. Zpravidla touto soustavou bývají zeměpisné poledníky nebo rovnoběžky. Nelze definovat ekvidistantní zobrazení, které by nezkreslovalo žádné délky. Ekvivalentní zobrazení – nezkresluje plochy, zkreslení úhlů je však zde poměrně značné, což se projevuje zejména ve tvarech ploch. Geografické souřadnice – určení polohy bodu na ploše elipsoidu pomocí zeměpisné šířky 𝜑 a zeměpisné délky 𝜆. Geocentrické souřadnice [𝑋, 𝑌, 𝑍] - prostorový souřadný systém s počátkem ve středu elipsoidu, osa 𝑋 je vložena do průsečíku rovníku a roviny základního (nultého) poledníku, osa 𝑍 spojuje střed elipsoidu a severní pól a osa 𝑌 leží v rovině rovníku otočena o 900 proti směru hodinových ručiček od osy 𝑋. Rovinné souřadnice – určení polohy v rovině pomocí dvojice rovinných souřadnic [𝑋, 𝑌] v ortogonálním souřadném systému. V některých zobrazeních (zobrazení UTM) se používá symbolika E, N (Easting, Northing), aby bylo zřejmé, v jakém směru rostou rovinné souřadnice. FI MUNI, Milan Drášil 2019 Úvod do GIS 28 Základní typy převodu geografických do rovinných souřadnic Základní typy převodu jsou dány 3D geometrickými tělesy, které jsou jednoduše převoditelné do roviny. Takovými tělesy jsou: - rovina - válec - kužel Princip je „geometricky jednoduchý“. Těleso umístíme do polohy vzhledem k modelu Země (koule, elipsoid) tak, aby projekce bodů na něj splňovala co nejvíce požadavky kladené na projekci. V praxi tento „jednoduchý“ požadavek vede k velmi komplikovaným konstrukcím, problémem se zabývá celý obor, Matematická kartografie (zkuste si třeba znovu sestavit triviální diferenciální rovnici pro Mercatorovo zobrazení, aniž byste listovali zpět). Například, umístění kuželu pro projekci S-JTSK je takové, aby délkové zkreslení na území bývalé ČSR bylo minimální. Jednotlivé body potom promítáme na těleso od středu Země. Kartografický souřadný systém – je soubor těchto údajů: - Geodetické datum (elipsoid, referenční bod, datum určení). - Souřadný systém geografických souřadnic 𝜑, 𝜆 (volba základního poledníku). - Zobrazovací rovnice z geografických do rovinných souřadnic [𝜙, 𝜆] → [𝑥, 𝑦]. FI MUNI, Milan Drášil 2019 Úvod do GIS 29 Příklad: Souřadný systém „Popular Mercator“ je definován: Revision date: 2008-03-13 Ellipsoid: Popular Visualisation Sphere (𝑅 = 6378137). Prime meridian: Greenwich Remarks: Used only for Web approximate mapping and visualisation. Not recognised by geodetic authorities. 𝑦 = 𝑙𝑛 (𝑡𝑎𝑛 ( 𝜋 4 + 𝜙 2 )) 𝑅 𝑥 = 𝜆𝑅 Viz https://epsg.io/6055-datum 2.2. NEJPOUŽÍVANĚJŠÍ SOUŘADNÉ SYSTÉMY V ČR Civilní souřadnicový systém S-JTSK je určen – Besselovým elipsoidem z roku 1841 s referenčním bodem Herrmanskogel, zeměpisné délky se určují od ferrského poledníku, zobrazovací rovnice dvojitého konformního kuželového zobrazení v obecné poloze (Křovákovo zobrazení) s volbou délkového faktoru 0.9999 pro snížení vlivu délkového zkreslení. Vojenský souřadnicový systém S-42 je určen Krasovského elipsoidem z roku 1942 s referenčním bodem Pulkovo, zeměpisné délky se měří od Greenwiche, zobrazovací rovnice Gaussova–Krügerova zobrazení s opakovatelností vždy pro šestistupňové poledníkové pásy. Od r. 2005 je nahrazen WGS84 – zobrazení UTM (Universal Transversal Mercator) Souřadný systém WGS 84 - Word Geodetic System 1984 Systém byl původně definován Ministerstvem obrany USA pro obranné účely, dnes je celosvětově používanou technologií prostorové lokalizace. UTM (Universal Transversal Mercator) – Systém transverzálního válcového zobrazení používající elipsoid WGS 84. Jsou použity šestistupňové pásy. Cassini-Soldner – Válcové zobrazení na Zachově elipsoidu s referenčními body Sv. Štěpán, resp. Gusterberg. Definováno v Rakousko Uherské monarchii pro stabilní katastr v měřítkách 1:2800 a 1:1440 (dodnes používaná). ETRS-LAEA: ETRS89 (Lambert Azimuthal Equal Area) – Elipsoid ETRS89 (téměř shodný s WGS 84), azimutální zobrazení používané pro střední a malá měřítka. FI MUNI, Milan Drášil 2019 Úvod do GIS 30 2.3. TRANSFORMACE SOUŘADNÝCH SYSTÉMŮ Problém: Mapy různých kartografických zobrazení transformovat do zobrazení cílového. Základem je znalost, jak transformovat bod. Příklad: V systému, který je provozován v kartografické projekci S-JTSK (národní projekce ČR) požadujeme, aby lokalizoval objekt o souřadnicích zadaných v projekci WGS84 (např. GPS). Přímá transformace: 1. Zdrojové souřadnice [𝑥, 𝑦] převedeme na geografické souřadnice zdrojového systému [𝜑, 𝜆]. 2. [𝜑, 𝜆] transformujeme do cílových geografických souřadnic [𝜑′, 𝜆′] 3. Geografické souřadnice [𝜑′, 𝜆′] převedeme do cílového rovinného zobrazení [𝑥′, 𝑦′] Transformace geografických souřadnic: 1. Geografické souřadnice [𝜑, 𝜆] převedeme na geocentrické souřadnice[𝑥 𝑠, 𝑦 𝑠, 𝑧 𝑠] 2. Pro převod mezi dvěma systémy geocentrických souřadnic použijeme tzv. Helmertovu prostorovou transformaci: 𝑥 = (1 + 𝑚)( 𝑥 𝑠 + 𝛾𝑦 𝑠 − 𝛽𝑧 𝑠) + 𝛥𝑥 𝑦 = (1 + 𝑚)(−𝛾𝑥 𝑠 + 𝑦 𝑠 + 𝛼𝑧 𝑠) + 𝛥𝑦 𝑧 = (1 + 𝑚)( 𝛽𝑥 𝑠 − 𝛼𝑦 𝑠 + 𝑧 𝑠) + 𝛥𝑧 3. Geocentrické souřadnice z 1. převedeme na geografické.Poznámka: Parametry Helmertovy prostorové transformace jsou získávány ze sad identických bodů, které se mohou lišit vlivem nehomogenity kartografické projekce pro různá území. Lze interpretovat jako: - vektor posunu – [𝛥𝑥, 𝛥𝑦, 𝛥𝑧] - úhel (resp. 𝑐𝑜𝑠) otočení pro jednotlivé osy - [𝛼, 𝛽, 𝛾] - změna měřítka - 𝑚 Tím dostáváme velmi snadno transformaci inverzní, jednoduše změníme znaménka všech parametrů. FI MUNI, Milan Drášil 2019 Úvod do GIS 31 Poznámka: Pomocí Helmertovy transformace jednoduše odvodíme transformace geocentrických souřadnic pro „neznámé“ elipsoidy. Stačí například znát pro všechny používané elipsoidy transformační klíč pro jeden pevný elipsoid, například WGS84. Pro transformaci geocentrických souřadnic z elipsoidů: 𝐸𝐿1 → 𝐸𝐿2 použijeme postup: 𝐸𝐿1 → 𝑊𝐺𝑆84 → 𝐸𝐿2 Polynomiální transformace Ze znalosti identických bodů ve zdrojovém a cílovém rovinném zobrazení [𝑋𝑖, 𝑌𝑖] → [𝑋′𝑖, 𝑌′𝑖] určíme koeficienty polynomiální transformace a tuto potom použijeme pro jednotlivé body. Používáme polynomy do 3. stupně, vyšší stupeň vede k nestabilitě řešení (omezený počet platných cifer). Používá se v případech, kdy: - Není známá zdrojová kartografická projekce prostorových dat. - Zdrojová projekce je sice známá, avšak je zatížena takovou chybou, že obecná polynomiální transformace dává lepší výsledky. - Zdrojová projekce je známá, ale její výpočet je příliš náročný vzhledem k počtu geometrických elementů a polynomiální transformace výsledek významně nezkreslí. V tomto případě použijeme kartografickou transformaci pro generování sítě odpovídajících si bodů (např. rohové body dotazového okna) a tyto body potom použijeme pro výpočet transformačního klíče. Základní typy polynomiálních transformací: Lineárni: 𝑥 = 𝑓(𝑢, 𝑣) = 𝑎1 𝑢 + 𝑏1 𝑣 + 𝑐1 𝑦 = 𝑔(𝑢, 𝑣) = 𝑎2 𝑢 + 𝑏2 𝑣 + 𝑐2 Bilineární: 𝑓(𝑢, 𝑣) = 𝑎1 𝑢 + 𝑏1 𝑣 + 𝑐1 𝑢𝑣 + 𝑑1 𝑔(𝑢, 𝑣) = 𝑎2 𝑢 + 𝑏2 𝑣 + 𝑐2 𝑢𝑣 + 𝑑2 Obecně - Kvadratická, kubická, obecná polynomiální … Nalezení koeficientů polynomiální transformace: Vstup: Dva seznamy „odpovídajících“ si bodů: [𝑥1, 𝑦1], . . , [𝑥 𝑛, 𝑦𝑛] a [𝑢1, 𝑣1], . . , [𝑢 𝑛, 𝑣 𝑛] FI MUNI, Milan Drášil 2019 Úvod do GIS 32 Výstup: Seznam koeficientů polynomiální transformace zvoleného stupně. Transformujeme [𝑢, 𝑣] do [𝑥, 𝑦] s co nejmenší chybou, tedy: ∑ 𝑑2 ([𝑥𝑖, 𝑦𝑖], [𝑓(𝑢𝑖, 𝑣𝑖), 𝑔(𝑢𝑖, 𝑣𝑖)]) = 𝑚𝑖𝑛 𝑛 𝑖=1 kde 𝑑2 ([𝑥1, 𝑦1], [𝑥2, 𝑦2]) = (𝑥1 − 𝑥2)2 + (𝑦1 − 𝑦2)2 Lineární regrese (příklad pro lineární transformaci): Σ = ∑(𝑎1 𝑢𝑖 + 𝑏1 𝑣𝑖 + 𝑐1– 𝑥𝑖)2 + (𝑎2 𝑢𝑖 + 𝑏2 𝑣𝑖 + 𝑐2– 𝑦𝑖)2 𝑛 𝑖=1 = 𝑚𝑖𝑛 Výraz Σ derivujeme podle proměnných 𝑎1. . 𝑐2 a derivaci položíme rovnu 0 (hledání extrémů funkcí více proměnných). 𝜕Σ 𝜕𝑎1 = ∑ 2(𝑎1 𝑢𝑖 + 𝑏1 𝑣𝑖 + 𝑐1 – 𝑥𝑖). 𝑢𝑖 = 0 𝜕Σ 𝜕𝑏1 = ∑ 2(𝑎1 𝑢𝑖 + 𝑏1 𝑣𝑖 + 𝑐1 – 𝑥𝑖). 𝑣𝑖 = 0 𝜕Σ 𝜕𝑐1 = ∑ 2(𝑎1 𝑢𝑖 + 𝑏1 𝑣𝑖 + 𝑐1 – 𝑥𝑖) = 0 Dostáváme soustavu tzv. normálních rovnic (pro 𝑔(𝑢, 𝑣) analogicky, místo indexů „1“ vystupují indexy „2“ a 𝑥 nahradíme 𝑦): 𝑎1Σ𝑢𝑖 2 + 𝑏1Σ𝑢𝑖 𝑣𝑖 + 𝑐1Σ𝑢𝑖 = Σ𝑥𝑖 𝑢𝑖 𝑎1Σ𝑢𝑖 𝑣𝑖 + 𝑏1Σ𝑣𝑖 2 + 𝑐1Σ𝑣𝑖 = Σ𝑥𝑖 𝑣𝑖 𝑎1Σ𝑢𝑖 + 𝑏1Σ𝑣𝑖 + 𝑐1 𝑛 = Σ𝑥𝑖 Tedy maticově: [ Σ𝑢𝑖 2 Σ𝑢𝑖 𝑣𝑖 Σ𝑢𝑖 Σ𝑢𝑖 𝑣𝑖 Σ𝑣𝑖 2 Σ𝑣𝑖 Σ𝑢𝑖 Σ𝑣𝑖 𝑐1 𝑛 ] [ 𝑎1 𝑏1 𝑐1 ] = [ Σ𝑥𝑖 𝑢𝑖 Σ𝑥𝑖 𝑣𝑖 Σ𝑥𝑖 ] FI MUNI, Milan Drášil 2019 Úvod do GIS 33 Příklad: Ukázka zdrojového kódu pro převod geografických souřadnic do systému SJTSK. Je zřejmé, že použití např. bilineární transformace je mnohem méně výpočtově náročné, v některých případech (např. pro rastrový obrázek 1024*1024 pixelů) hodně významně. /* konstanty Besselova elipsoidu */ Double consta = 6377397.15508; Double constb = 6356078.96290; /* konstanty */ Double constalfa = 1.00059749837159;Double constr = 6380703.610617; Double constk = 0.996592486879232; /* uhel v radianech odpovidajici 45 */ Double constu45 = 0.785398163397448; /* uhel v radianech odpovidajici 78:30 */ Double constu7830 = 1.370083462815548; /* uhel v radianech odpovidajici 59:42:42.6969 */ Double constuk = 1.042168563853224; /* uhel v radianech odpovidajici 42:31:31.41725 */ Double constvk = 0.742208135432484; /* uhel v radianech odpovidajici 17:39:59.7354 */ Double constferra = 0.308340218368665; public override bool FromEllipsoid( double lat, double lon, double alt, ref double x, ref double y, ref double h) { try { Double dpom; Double de, dsfi; Double du, dv, ddeltav; Double ds, dd; Double dro, depsilon; Double xx, yy; /* prevod vzhledem k ferra */ lon += constferra; /* prevod lat,lon (Bessel) -> u,v (Gaussova koule) */ dv = constalfa * lon; de = Math.Sqrt(consta * consta - constb * constb) / consta; dsfi = Math.Sin(lat); dpom = (1 - de * dsfi) / (1 + de * dsfi); dpom = Math.Pow(dpom, de / 2); dpom = dpom * Math.Tan(lat / 2 + constu45); dpom = Math.Pow(dpom, constalfa) / constk; du = 2 * (Math.Atan(dpom) - constu45); /* prevod u,v (Gaussova koule) -> s,d (sirka,delka) */ ddeltav = constvk - dv; dpom = Math.Sin(constuk) * Math.Sin(du) + Math.Cos(constuk) * Math.Cos(du) * Math.Cos(ddeltav); ds = Math.Asin(dpom); dd = Math.Asin(Math.Sin(ddeltav) * Math.Cos(du) / (ds)); /* prevod s,d -> x,y (rovinne, JTSK) */ dpom = Math.Tan(constu7830 / 2 + constu45) / Math.Tan(ds / tu45); dpom = Math.Pow(dpom, Math.Sin(constu7830)); dro = dpom * 0.9999 * constr / Math.Tan(constu7830); depsilon = Math.Sin(constu7830) * dd; x = dro * Math.Cos(depsilon); y = dro * Math.Sin(depsilon); h = alt; return true; } catch { return false; } } FI MUNI, Milan Drášil 2019 Úvod do GIS 34 Katalog kódů projekcí a jejich hlavních vlastností lze nalézt například na: http://www.epsg-registry.org/ Poznámka: Při vytváření normálních rovnic nám hodně pomůže jejich maticová reprezentace a obecně odvozená metoda nejmenších čtverců, viz příloha 9.2 2.4. TRADIČNÍ MAPY A POJMY SOUVISEJÍCÍ S GIS - Mapování – vytváření map měřením nebo fotogrammetrickým mapováním pomocí geodetických základů – bodů geodetických sítí. Mapa je 2D obraz zemského povrchu. - Měřítko mapy – v GIS kontextu neznamená doslova poměr skutečné a zobrazované velikosti objektů. Měřítkem mapy se spíše rozumí úroveň územní podrobnosti obsahu geografického informačního systému. - Dálkový průzkum Země (DPZ) – získávání informací o zemském povrchu a jeho blízkém okolí pomocí snímacích zařízení (kamery, skenery) umístěných v letadlech nebo satelitech Země. Využitelné a využívané jsou také drony. - Topografická mapa – je grafická prezentace (zobrazení) části zemského povrchu se standardizovaným obecným obsahem (voda, lesy, komunikace, objekty viditelné na zemském povrchu...) - Tematická mapa – zobrazení geografických dat a jevů v topografickém podkladu pomocí grafické reprezentace prostorových dat: bodů, linií a areálů. Metody reprezentace: bodové značky, lokalizované kartodiagramy, kartodiagramy, symbolika čar, kartogramy. FI MUNI, Milan Drášil 2019 Úvod do GIS 35 2.5. MAPOVÉ DÍLO V ČR Mapy velkých měřítek do 1:5000 - Katastrální mapy (mapy stabilního katastru) v systému Cassini–Soldner (počátek Gusterberg v Čechách, Sv. Štěpán na Moravě) v sáhových měřítkách 1:2880, 1:1440, 1:720 (měřítko je odvozeno ze vztahu 1 jitro – 1600 čtverečních sáhů – je zobrazeno jako čtvereční palec), ale v dekadických měřítkách. - Katastrální mapy v S-JTSK (systém jednotné trigonometrické sítě katastrální – Křovák), měřítka 1:1000 ve městech (intravilán), 1:2000 v extravilánu vznikaly po roce 1928. - Státní mapa odvozená v měřítku 1:5000, systém JTSK, obsah: vlastnické hranice, polohopis (vnitřní kresba). - Digitální katastrální mapa – mapa vedená digitálně, udržovaná katastrálními úřady. Státní mapové dílo velkých měřítek v České republice vznikalo v průběhu dvou století. Mapové dílo je charakteristické svou rozmanitostí a rozdílnou kvalitou (především vzhledem k přesnosti a aktuálnosti mapy). Mapy středních měřítek 1 : 10000 až 1 : 200 000 - Základní mapa středního měřítka – v měřítkách 1:10000, 1:25000, 1:50000, 1:100000, 1:200000 s obsahem topografické mapy, v digitální formě mapové dílo ZABAGED. - Topografická mapa Generálního štábu armády ČR, měřítko 1:25000 (v některém území i 1:10000) FI MUNI, Milan Drášil 2019 Úvod do GIS 36 3. Datové sklady geoinformačních systémů 3.1. TYPY PROSTOROVÝCH DAT Vektorová data – reprezentují objekty pomocí datových struktur, jejichž základní položkou je bod 2–rozměrného spojitého (euklidovského) prostoru. Termínem “spojitý” myslíme spojitý, až na technické omezení použité počítačové aritmetiky. Rastrová data – podmnožina 2D prostoru je pravidelně rozdělena (většinou čtvercovou) sítí, každý element této sítě je nositelem tematické části (geografické) informace. Prostorová lokalizace je určena indexem elementární složky sítě, popřípadě jeho zobrazením do cílového (kartografického) souřadného systému. Grid data – Základem této reprezentace je opět pravidelná síť položená na 2D prostor. Rozdíl oproti rastrovým datům spočívá v tom, že tematická část informace je získávána na základě předem definované sítě, kterou je rozděleno zájmové území. Topologie – vymezuje vztahy mezi entitami (objekty) systému, aniž by musela obsahovat umístění objektu v prostoru. Například informační systémy o spojení míst silniční sítí nevyžadují přesné umístění uzlů v prostoru, pracují pouze s relací na množině všech uzlů. 3.2. RASTROVÉ DATOVÉ SKLADY Využívají běžné formáty obrazových dat (bmp, png, jpeg tif, gif,ecw..). Některé z nich mohou mít informace o vztahu k souřadnicím kartografické projekce přímo ve své hlavičce (ecw, tif) u ostatních se používá „doplňující“ textový soubor (*.tfw), který obsahuje parametry lineárního převodu z/do pixelových do/z kartografických souřadnic. Nechť [𝑖, 𝑗]jsou pixelové souřadnice, [𝑥, 𝑦] kartografické souřadnice, potom převod z pixelových do kartografických je reprezentován lineárními rovnicemi: 𝑥 = 𝑎 𝑥 𝑖 + 𝑏 𝑥 𝑗 + 𝑐 𝑥 𝑦 = 𝑎 𝑦 𝑖 + 𝑏 𝑦 𝑗 + 𝑐 𝑦 FI MUNI, Milan Drášil 2019 Úvod do GIS 37 TFW soubor je potom textový soubor obsahující parametry 𝑎 𝑥, 𝑏 𝑥, 𝑎 𝑦, 𝑏 𝑦, 𝑐 𝑥, 𝑐 𝑦 Příklad TFW souboru: 12.7011007 0 0 -12.7011007 -775450 -965225 Poznámka: V naprosté většině případů je 𝑏 𝑥 = 0 a 𝑎 𝑦 = 0. To znamená, že osy obrázku jsou rovnoběžné s osami kartografického zobrazení. 3.3. VEKTOROVÉ DATOVÉ SKLADY Pro fyzickou reprezentaci je možné použít vlastních datových struktur a ukládat je přímo v souborovém systému operačního systému. Přes nesporné výhody tohoto přístupu, jako je optimalizace uložení prostorové složky informace, převažují nevýhody, zejména aplikační závislost. Běžnější přístup je použití robustních databázových strojů, které vektorovou geometrii běžně používají (Oracle, Microsoft SQL-Server, Postgresql, MySQL). Není jednotný standard ukládání vektorové geometrie. Nejpoužívanější veřejné formáty: Shape file (systém ARC/INFO fy. ESRI) Geograficky vztažená informace je obsažena ve trojici souborů - *.shp prostorová informace - *.shx prostorový index - *.dbf popisná informace a topologické vazby Základní geometrické primitivy: Point bod MultiPoint body Line lomená čára MultiLine lomené čáry Polygon areál MultiPolygon areály Vše 2D a 3D varianty. Formát neobsahuje symboliku (barva, síla, styl linií, výplň polygonu). Tu obsluhuje aplikace na základě popisných informací. Norma neobsahuje „heterogenní“ kolekce geometrií a tím ani definici mapových symbolů. Jako mapové symboly jsou použity speciální znakové sady (fonty). FI MUNI, Milan Drášil 2019 Úvod do GIS 38 Shape file byl jedním z prvních obecně používaných formátů, dodnes je podporován většinou „velkých“ geoinformačních systémů. DGN file, norma IGDS (Intergraph, Bentley) Geometrická informace je obsažena v souboru *.DGN, soubor sám o sobě nenese popisnou informaci, ta je uložena v relační databázi (nebo *.dbf souboru), DGN soubor obsahuje pouze tzv. „link“ = společný klíč v souboru/databázi. Základní geometrické primitivy: CELL_HEADER_ELM Hlavička buňky LINE_ELM Úsečka LINE_STRING_ELM Lomená čára SHAPE_ELM Polygon TEXT_ELM Text TEXT_NODE Textový uzel CURVE_ELM Křivka CMPLX_STRING_ELM Komplexní linie CMPLX_SHAPE_ELM Komplexní polygon ELLIPSE_ELM Elipsa ARC_ELM Oblouk POINT_STRING_ELM Body BSPLINE_ELM B-spline DIMENSION_ELM Kóta - Komplexní linie se mohou skládat z různých segmentů – např. lomených čar a oblouků. - Geometrie obsahuje symboliku geometrických primitiv. - Typ CELL může opět obsahovat CELL. Podobné vlastnosti mají i ostatní CAD formáty (DXF, DWG...) Bohatý repertoár geometrických primitiv CAD formátů umožňuje i definici mapových symbolů. Problematické jsou interpretace: - nelineárních geometrií v různých klientských aplikacích (např. B-spline) - složitějších zobrazovacích symbolik (např, linie vzorovaná symbolem – směr toku media ve vodovodní síti). - operace na složitějšími geometriemi (např. průnik polygonů, jejichž hranice se skládají z lomených čar či eliptických oblouků …). Příklad: Definice složitého mapového symbolu (Sídla Městského úřadu v GIS pro města). Je zřejmé, že v takovém případě je výhodné mít k dispozici robustní CAD formát vektorových dat, u kterého si každý geometrický element nese i výchozí symboliku kresby. FI MUNI, Milan Drášil 2019 Úvod do GIS 39 ORACLE 7x (Spatial Data Option): Geometrie je reprezentována čtyřmi tabulkami: _SDOLAYER – obsahuje meta údaje pro prostorovou indexaci _SDODIM – obsahuje rozsah pro jednotlivé dimenze geometrie _SDOGEOM – obsahuje vlastní geometrii _SDOINDEX – obsahuje prostorové indexy objektů možné typy geometrie jsou: bod, lomená čára a areál, Jednalo se o první pokus o standardizaci geometrie – ten se vlivem denormalizace uložení souřadnic ukázal jako slepá ulička, v současné době není rozvíjen. ORACLE 8x a výše – datový typ SDO_GEOMETRY: UNKNOWN_GEOMETRY neznámá geometrie POINT bod LINESTRING lomená čára POLYGON areál (oblast) COLLECTION kolekce geometrií MULTIPOINT Body MULTILINESTRING kolekce linií MUTIPOLYGON kolekce areálů - Linie jsou tvořeny sekvencemi bodů a kruhových oblouků. - Typ COLLECTION nemůže obsahovat typ COLLECTION. - Definice neobsahuje symboliku geometrických primitiv. FI MUNI, Milan Drášil 2019 Úvod do GIS 40 Open GIS Consortium, Inc. Je sdružení soukromých, veřejných organizací (univerzit, komerčních firem…) se zájmem na vybudování „standardu“ struktur a služeb v prostorových datech. I přes byrokratickou těžkopádnost způsobenou množstvím subjektů, které se účastní vývoje standardů lze konstatovat, že standardy OGC jsou poměrně rozšířeny a podporovány (např. Oracle podporuje konverzi datových typů, prostorová informace MS-SQL a MySQL je přímou implementací OGC). Our Vision - is a world in which everyone benefits from geographic information and services made available across any network, application, or platform. Our Mission - is to deliver spatial interface specifications that are openly available for global use. OGC Well known binary: Je norma binárního ukládání vektorové geometrie včetně C-definic struktur: // Basic Type definitions // byte: 1 byte // uint32: 32 bit unsigned integer (4 bytes) // double: double precision number (8 bytes) // Building Blocks : Point, LinearRing Point { double x; double y; }; LinearRing { uint32 numPoints; Point points[numPoints]; } enum wkbByteOrder { wkbXDR = 0, // Big Endian wkbNDR = 1 // Little Endian }; WKBPoint { byte byteOrder; uint32 wkbType; // 1 Point point; }; FI MUNI, Milan Drášil 2019 Úvod do GIS 41 WKBLineString { byte byteOrder; uint32 wkbType; // 2 uint32 numPoints; Point points[numPoints]; }; WKBPolygon { byte byteOrder; uint32 wkbType; // 3 uint32 numRings; LinearRing rings[numRings]; }; WKBMultiPoint { byte byteOrder; uint32 wkbType; // 4 uint32 num_wkbPoints; WKBPoint WKBPoints[num_wkbPoints]; }; WKBMultiLineString { Byte byteOrder; uint32 wkbType; // 5 uint32 num_wkbLineStrings; WKBLineString WKBLineStrings[num_wkbLnStrgs]; }; wkbMultiPolygon { byte byteOrder; uint32 wkbType; // 6 uint32 num_wkbPolygons; WKBPolygon wkbPolygons[num_wkbPolygons]; } WKBGeometry { union { WKBPoint point; WKBLineString linestring; WKBPolygon polygon; WKBGeometryCollection collection; WKBMultiPoint mpoint; WKBMultiLineString mlinestring; WKBMultiPolygon mpolygon; } }; WKBGeometryCollection { Byte byte_order; uint32 wkbType;// 7 uint32 num_wkbGeometries; WKBGeometry wkbGeometries[num_wkbGeoms]; } FI MUNI, Milan Drášil 2019 Úvod do GIS 42 Základní datové typy WKB neumožňují vykreslit plnohodnotnou mapu, neobsahují: - Symboliku (grafická reprezentace – barva, síla, tloušťka) geometrických elementů. - Reprezentaci bodových prvků, mapové symboly. - Texty (velikost, rotace, font …). Definice WKB neobsahuje definici kruhových oblouků, což může působit obtíže při práci v měřítkách, kde platí Euklidovská geometrie a „kružítko“ běžně používáme. Rekurzivní definice WKBGeometryCollection umožňuje i definici komplexních mapových symbolů. GML – Geographic Markup Language: Původně formát WKB v XML formátu, od verze 3.2 značně rozšířená (křivky…). Geometrie mohou obsahovat i kód kartografické projekce, a to teoreticky umožňuje v rámci jedné datové sady používat projekcí více. V EU (projekt INSPIRE) byla specifikace GML převzata jako norma pro výměnu vektorových prostorových dat. Formát GML definuje XML reprezentaci geometrie, popisnou část informace nijak neomezuje. Vzhledem k masivní podpoře XML parserů mnoha vývojových prostředí se zřejmě jedná nejuniverzálnější výměnný formát. Pro ukládání prostorových informací je však nevhodný, jeho velikost je oproti binárnímu formátu až desetinásobná. Příklad GML: ¨ 541484.77 1113662.92 541484.36 1113666.36 541476.26 1113665.3 541469.89 1113664.51 541469.62 1113662.23 541453.022 1113662.019 541452.31 1113662.01 541451.77 1113664.1 541446.27 1113662.85 541446.96 1113668.85 541438.25 1113668.81 541440.74 1113661.7 541432.057 1113660.152 541427.45 1113659.33 541426.737 1113661.296 541425.54 1113664.6 541419.741 1113662.802 541401.42 1113657.12 541397.74 1113655.81 541391.63 1113653.68 541394.736 1113649.478 541395.23 1113648.81 541407.55 1113652.235 541427.96 1113657.91 541438.177 1113658.825 541438.209 1113658.828 541440.8 1113659.06 541469.46 1113660.94 541484.77 1113662.92 … FI MUNI, Milan Drášil 2019 Úvod do GIS 43 OGC specifikace pro ukládání prostorových informací v relačních databázích: (OpenGIS Simple Features Specification For SQL metamodel) GEOMETRY_COLUMNS F_TABLE_CATALOG F_TABLE_SCHEMA F_TABLE_NAME F_GEOMETRY_COLUMN G_TABLE_CATALOG G_TABLE_SCHEMA G_TABLE_NAME STORAGE_TYPE GEOMETRY_TYPE COORD_DIMENSION MAX_PPR SRID SPATIAL_REFERENCE_SYSTEMS SRID AUTH_NAME AUTH_SRID SRTEXT GEOMETRY_COLUMNS GID ESEQ ETYPE SEQ X1 Y1 … … X Y GEOMETRY_COLUMNS GID XMIN YMIN XMAX YMAX WKB_GEOMETRY Feature Table/View GID (Geometry Column) or FI MUNI, Milan Drášil 2019 Úvod do GIS 44 3.4. WEBOVÉ SLUŽBY A JEJICH STANDARDY 3.4.1. OGC – WEB MAP SERVICE (WMS): Je webová služba poskytující mapu v zadaném výřezu a zadané kartografické projekci. Obsahuje (mimo jiné) dva základní dotazy: GetCapabilities – vrací metainformace o službě ve formátu XML. Metainformace obsahují zejména výčet mapových vrstev, podporované kartografické projekce a rozsah měřítek, pro které jsou tato data poskytována. WMS KN - CUZK EPSG:102067 EPSG:32633 EPSG:32634 obrazy_parcel Obrazy parcel … FI MUNI, Milan Drášil 2019 Úvod do GIS 45 GetMap – vrátí obrázek s požadovanou mapou. Příklad: http://wms.cuzk.cz/wms.asp? request=GetMap& version=1.1.1& layers=dalsi_p_mapy_i,hranice_parcel_i, obrazy_parcel_i,OMP,parcelni_cisla_i& srs=EPSG:102067& bbox=-815350,-1096562,-815058,-1096388& width=1371& height=815& bgcolor=0x999999& Format=image/gif FI MUNI, Milan Drášil 2019 Úvod do GIS 46 3.4.2. OGC – WEB MAP TILE SERVICE (WMTS) Je služba, která poskytuje „dlaždice map“ v definované kartografické projekci. Je používána tam, kde se obsah mapy příliš nemění v čase a dlaždice je možné předem připravit. Na podobném principu pracují i robustní světové servery (Google maps, Bing maps). Služba obsahuje opět mimo jiné dva základní dotazy: GetCapabilities – vrací metadata služby ve formátu XML. Metadata obsahují informace o dlaždicových sadách (TileMatrixSet, TileMatrix), jejich měřítkách a umístění v souřadném systému. Příklad: URL?request=GetCapabilities&version=1.1.1& service=WMTS OrtoFoto ORTOMAPA EPSG:5514 Zoom_09 Zoom_09 14285.714285714286 -931496 -819102 256 256 512 512 Zoom_10 Zoom_10 7142.8571428571431 -931496 -819102 256 256 1024 1024 FI MUNI, Milan Drášil 2019 Úvod do GIS 47 GetTile – vrací mapovou dlaždici. URL&service=WMTS&request=GetTile& TileMatrixSet=ORTOMAPA&TileMatrix=Zoom_11& TileRow=1024&TileCol=1024 3.4.3. OGC WEB FEATURE SERVICE (WFS) WFS je služba, která poskytuje vektorová data v normě GML a atributy libovolné struktury. GetCapabilities – vrací XML soubor popisující služby a seznam poskytovaných datových datové sad (feature). . . CP:CadastralParcel Cadastral parcel polygons Cadastral parcel polygons urn:ogc:def:crs:EPSG::102067 urn:ogc:def:crs:EPSG::102066 urn:ogc:def:crs:EPSG::2065 text/xml; subtype=gml/3.2.1 10 43 22 55 FI MUNI, Milan Drášil 2019 Úvod do GIS 48 . . DescribeFeatureType – vrací XSD šablonu pro požadovaný feature (= typu datové sady). ListStoredQueries – vrací XML soubor se seznamem uložených dotazů a návratových typů. CP:CadastralBoundary CP:CadastralParcel CP:CadastralZoning CP:CadastralBoundary CP:CadastralParcel DescribeStoredQueries – vrací XML soubor s popisem parametrů uložených dotazů Cadastral Parcel Boundaries GetFeature – Vrací sadu prostorových (geometrických) dat. Požadavek je definován parametrem Filter, což je XML fragment definující prostorovou a/nebo logickou podmínku na požadovaná data. Je možné definovat velmi složité podmínky. FI MUNI, Milan Drášil 2019 Úvod do GIS 49 557658 106440 557489 106330 CP:zoning KU601977 CP:label 198 215 CP:beginLifespanVersion 2010-03-08T12:46:20 CP:zoning Zdrojový kód pro poskytovatele/konzumenta je většinou generován strojově. Výsledný kód má totiž řádově 104 až 105 řádků. Je téměř nemyslitelné, aby kód, která akceptuje tuto službu byl sestaven ručně a bezchybně programátorem. I tak se často musí generovaný kód ručně upravovat, zejména, když jsou třídy generovány odděleně. FI MUNI, Milan Drášil 2019 Úvod do GIS 50 Příklad – implementace tříd směrnice INSPIRE (EU) Třídy jmenného prostoru (namespace) gml byly vygenerovány dříve a jsou fyzicky v jiném sestavení (DLL), než objekty jmenného prostoru CP v konkrétní implementaci WFS. Má to logiku, nebudeme přece opakovat gml kód pro každou WFS službu zvlášť. Při vygenerování tříd pro konkrétní implementaci, se však dostaneme do situace, že třída gml:FeaturePropertyType musí akceptovat objekty ze jmenného prostoru CP. To není možné, sestavení pro gml prostě „nezná“ objekty z tohoto prostoru. Možné řešení je naznačeno v následujícím diagramu tříd: FI MUNI, Milan Drášil 2019 Úvod do GIS 51 4. Efektivní přístup k prostorovým datům 4.1. VYMEZENÍ PROBLÉMU Efektivním přístupem k prostorovým datům rozumíme vyhledávací techniky nad prostorovými daty, kdy dotazy mají opět prostorový charakter. Například, vyhledej všechny objekty, které leží uvnitř obdélníku, který je vymezen obrazovkou počítače. Problém vyhledání můžeme popsat termíny formální logiky prvního řádu (FOL). Za vybudování FOL vděčíme německému matematikovi a logikovi Gottlobu Fregemu (1848-1925), na jejím základě byl vybudován „betonový“ základ matematiky – Teorie množin. Pro naše účely postačí následující „zkratka“: - Oblast zájmu (universum) je popsána souborem relací – predikátů, některé z nich mohou být funkce. - FOL poskytuje aparát pro „správné“ sestavení výrazů a formulí s proměnnými. Proměnné mohou být kvantifikovány (univerzálně – „pro každý“, existenčně – „existuje“) Vyhledání v relaci 𝑈 podle formule 𝛷 je potom podmnožina 𝑉 ⊆ 𝑈: 𝑉 = {𝑢 ∈ 𝑈|𝛷(𝑢)} Příklad – Problém příslušnosti prvku k množině: Pro libovolnou množinu 𝑉 můžeme definovat jednoduchou formuli: 𝛷: 𝑥 = 𝑎 kde 𝑥 je proměnná a 𝑎 je libovolná konstanta (z univerza). Aplikací formule 𝑉 = {𝑢 ∈ 𝑈|𝛷(𝑢)} dostaneme odpověď na to, zda 𝑎 ∈ 𝑈, tj. 𝑈 obsahuje 𝑎. Příklad – Rozsahový dotaz na uspořádané množině: Nechť (𝑇, ≤) je úplně uspořádaná množina (její každé dva elementy lze porovnat) 𝑈 ⊆ 𝑇. Rozsahovým dotazem potom rozumíme aplikaci formule (𝑎, 𝑏 ∈ 𝑇 jsou konstanty): 𝛷: 𝑎 ≤ 𝑥 ≤ 𝑏, kde 𝑎 ≤ 𝑏 Příklad – Rozsahový dotaz na body ve 2D prostoru: Nechť (𝑇, ≤) je úplně uspořádaná množina, 𝑈 ⊆ 𝑇2 . 2D rozsahovým dotazem potom rozumíme formuli: FI MUNI, Milan Drášil 2019 Úvod do GIS 52 𝛷: 𝑥 𝑚𝑖𝑛 ≤ 𝑥 ≤ 𝑥 𝑚𝑎𝑥 ∧ 𝑦 𝑚𝑖𝑛 ≤ 𝑦 ≤ 𝑦 𝑚𝑎𝑥 kde 𝑥 𝑚𝑖𝑛, 𝑥 𝑚𝑎𝑥, 𝑦 𝑚𝑖𝑛, 𝑦 𝑚𝑎𝑥 ∈ 𝑇, 𝑥 𝑚𝑖𝑛 ≤ 𝑥 𝑚𝑎𝑥, 𝑦 𝑚𝑖𝑛 ≤ 𝑦 𝑚𝑎𝑥 jsou konstanty. Typ dotazů 𝛷 jsou obdélníky (okna) rovnoběžné s osami souřadného systému). Aplikací formule dostaneme všechny body uvnitř „okna“. Příklad – Rozsahový dotaz na obdélníky ve 2D prostoru: Nechť (𝑇, ≤) je úplně uspořádaná množina 𝑈 ⊆ 𝑇4 s vlastností: ∀[𝑥1, 𝑦1, 𝑥2, 𝑦2] ∈ 𝑈(𝑥1 ≤ 𝑥2 ∧ 𝑦1 ≤ 𝑦2) Prvky množiny 𝑈 jsou obdélníky ve 2D prostoru. Formule, která vyhledá všechny prvky z 𝑈, které mají neprázdný průnik (incidují) s 2D rozsahem (také obdélníkem) má tvar: Φ: 𝑥1 ≤ 𝑥 𝑚𝑎𝑥 ∧ 𝑦1 ≤ 𝑦 𝑚𝑎𝑥 ∧ 𝑥2 ≥ 𝑥 𝑚𝑖𝑛 ∧ 𝑦1 ≥ 𝑦 𝑚𝑖𝑛 kde 𝑥 𝑚𝑖𝑛, 𝑥 𝑚𝑎𝑥, 𝑦 𝑚𝑖𝑛, 𝑦 𝑚𝑎𝑥 jsou konstanty s vlastnostmi z předešlého příkladu. Jednotný zdroj pro prostorovou indexaci geometrických objektů: (tj. struktur pro efektivní vyhledání) je minimální omezující obdélník geometrického objektu rovnoběžný s osami souřadného systému – MBR (Minimal Bounding Rectangle): tedy minima, resp. maxima lomových (definičních) bodů [𝑥 𝑚𝑖𝑛 , 𝑦 𝑚𝑖𝑛, 𝑥 𝑚𝑎𝑥, 𝑦 𝑚𝑎𝑥] V naprosté většině případů vystačíme s obdélníkovým dotazem: FI MUNI, Milan Drášil 2019 Úvod do GIS 53 Metoda, která realizuje tento dotaz je často nazývána primárním filtrem. Metoda, která realizuje přesnou odpověď je nazývána filtrem sekundárním. Vraťme se ze světa logiky a množin do světa informatiky. Všechny uvedené příklady lze triviálně řešit jedním průchodem množiny/seznamu 𝑈, tedy v lineární časové složitosti 𝑂(|𝑈|). Uvádění jiných metod má tedy smysl pouze v případě, že tento základní odhad nějak zlepšíme. Pro rozsahové výběry se většinou studuje časová složitost „zásahu“ prvního objektu, který splňuje podmínku rozsahového výběru. Běžné indexovací metody (tj. ty které jsou implementovány v RDBMS – např. B+ stromy) poskytují efektivní aparát pro vyhledávací problémy: - příslušnosti k množině - rozsahový dotaz ale samy o sobě neposkytují aparát vhodný k prostorovým dotazům: Příklad – incidence intervalů: Představme si tabulku, kde bude za zaměstnance zaznamenán jeho nástup do firmy a konec jeho pracovního poměru. Vedení firmy bude zajímat, kdo všechno byl zaměstnán ve firmě v určitém časovém období. Formálně: Máme soubor intervalů (1D obdélníků), a dotaz bude opět interval. Odpovědí budou všechny intervaly, které s dotazem incidují. Lineární indexovací metoda (tj. uspořádání podle jednoho či více klíčů), nám nepomůže, neboť nejhorší případ dotazu vede prohledání celého souboru. Problémy vyhledávání rozdělíme na dvě hlavní třídy: - problém statický - problém dynamický FI MUNI, Milan Drášil 2019 Úvod do GIS 54 Statický: - 𝑏𝑢𝑖𝑙𝑑(𝑈) vybuduje podpůrné struktury pro množinu 𝑈 - 𝑠𝑒𝑎𝑟𝑐ℎ(𝛷, 𝑈) odpoví na vyhledávací dotaz Dynamický: - 𝑖𝑛𝑠𝑒𝑟𝑡(𝑥, 𝑈) vloží do množiny 𝑈 nový objekt 𝑥 - 𝑑𝑒𝑙𝑒𝑡𝑒(𝑥, 𝑈) vymaže z množiny 𝑈 objekt 𝑥 - 𝑠𝑒𝑎𝑟𝑐ℎ(𝛷, 𝑈) odpoví na vyhledávací dotaz Dynamické řešení problému zároveň řeší statickou variantu (opakujeme funkci insert). Funkce 𝑠𝑒𝑎𝑟𝑐ℎ bývá většinou rozdělena na dvě části, a to - 𝛷𝑖 = 𝑖𝑛𝑖𝑡(𝛷, 𝑈) inicializace dotazu - 𝑜 = 𝑓𝑒𝑡𝑐ℎ(𝛷𝑖) vrací jeden objekt z množiny 𝑈 práce potom probíhá podle jednoduchého schématu: var f=init(q,U); while((object x=fetch(f))!=null) { zpracuj_objekt(x); } close(f); 4.2. METODA „GRID“: Spočívá v pravidelném rozdělení 2D prostoru, resp. zájmového území obdélníkovou sítí. Elementy sítě lze maticově indexovat, tím prostor „linearizovat“ a využít je tím pro primární prostorový filtr. .1 .3 .4.5 .6 4 3 2 1 7654321 .7 FI MUNI, Milan Drášil 2019 Úvod do GIS 55 . 3 . 4. 5 . 6 . 1 5 4 3 2 1 7654321 Dotazový obdélník . 7 Prostorový dotaz v GRIDu“: Prohledáváme pouze čtverce incidentní s dotazem, tedy (4,2) a (5,2), pro efektivní přístup ke čtvercům použijeme libovolnou vyhledávací metodu podporující rozsahový dotaz na úplně uspořádaných množinách. Tyto metody jsou standardně implementovány ve všech databázových strojích a lze je téměř okamžitě použít. Realizace GRID metody v prostředí SQL: Každé prostorové tabulce (tj. obsahující vektorovou geometrii např. ve formátu WKB) přiřadíme indexovací prostorovou tabulku. Tabulka s prostorovými daty: create table GTABLE ( GID number, XMIN number, YMIN number, XMAX number, YMAX number, WKB_GEOMETRY blob, ... constraint GTABLE_PK primary KEY (GID)); GID = GID ID = ID GTABLE GID XMIN YMIN XMAX YMAX WKB_GEOMETRY NUMBER NUMBER NUMBER NUMBER NUMBER BLOB GTABLE_IDX GID GRID_X GRID_Y NUMBER NUMBER NUMBER SPATIAL_QUERY ID XMIN YMIN XMAX YMAX NUMBER NUMBER NUMBER NUMBER NUMBER SPATIAL_QUERY_IDX ID GRID_X GRID_Y NUMBER NUMBER NUMBER FI MUNI, Milan Drášil 2019 Úvod do GIS 56 Tabulka grid indexů: create table GTABLE_IDX ( GID number, GRID_X number, GRID_Y number ); Omezení a prostorové indexy: alter table GTABLE_IDX add constraint GTABLE_IDX_PK primary key (gid,grid_x,grid_y); alter table GTABLE_IDX add constraint GTABLE_IDX_fk1 foreign key (GID) references GTABLE(GID) ON DELETE CASCADE; create index GTABLE_IDX_I1 on GTABLE_IDX(grid_x, grid_y); Triggerem zajistíme vkládání prostorových indexů: create trigger gtable_spatial before insert or update of x,y on GTABLE for each row begin xfrom:=GET_GRID_X(:NEW.XMIN); xto :=GET_GRID_X(:NEW.XMAX); yfrom:=GET_GRID_Y(:NEW.YMIN); yto :=GET_GRID_Y(:NEW.YMAX); pro xfrom<=i<=xto a yfrom<=j<=yto begin INSERT INTO GTABLE_IDX VALUES(:NEW.GID,i,j); end; end; / (funkce GET_GRID_X/Y vrací gridové indexy) FI MUNI, Milan Drášil 2019 Úvod do GIS 57 Implementace prostorového dotazu: Vytvoříme dotazovou tabulku: create table SPATIAL_QUERY ( id int, xmin int, ymin int, xmax int, ymax int, constraint SPATIAL_QUERY_PK primary key (id) ); a ostatní objekty (indexová tabulka, integritní omezení, trigger) stejně jako u tabulek s prostorovými daty. Prostorový dotaz pro obdélník [xmin,ymin,xmax,ymax] provedeme následovně: 1. Identifikace dotazu: Z databáze získáme nový (jednoznačný) klíč dotazu id, například ze sekvence. 2. Inicializace dotazu: insert into spatial_query values (id,xmin,ymin,xmax,ymax), vlivem triggeru SPATIAL_QUERY_SPATIAL automaticky vloží identifikace gridových čtverců to tabulky spatial_query_idx 3. Prostorový dotaz: select … from gtable A, gtable_idx B, spatial_query_idx C where A.GID=B.GID AND B.grid_x=C.grid_x AND B.grid_y=C.grid_y AND C.query_id=id; Ukončení prostorového dotazu: delete from spatial_query where id=id; Jaký mechanismus odstraňuje řádky z tabulky spatial_query_idx?) FI MUNI, Milan Drášil 2019 Úvod do GIS 58 Výhody vs. nevýhody GRID metody. + - velmi snadná implementace v prostředí RDBMS - snadné rozšíření na více dimenzí (?) - relativně snadná (resp. řešitelná implementace neobdélníkových dotazů) - netriviální odhad velikosti GRIDových čtverců, špatná volba má dramatické důsledky - nepravidelné chování při řádově rozdílné velikosti geometrických objektů 4.3. MODIFIKACE BINÁRNÍCH STROMŮ PRO PROSTOROVÉ VYHLEDÁVÁNÍ, KD STROMY Definice – Binární strom: - Nechť (𝑈, ≤) je úplně uspořádaná množina, 𝑉 𝑈 a |𝑉| < ∞ - Nechť (𝑉, 𝐸) je strom ve smyslu teorie grafů, takový, že každý jeho uzel obsahuje maximálně dva syny. - Každý uzel 𝑛 vyjma kořene má vlastnost „pravý“ a “levý“, označme 𝑛 𝐿, 𝑛 𝑅 levé a pravé následníky uzlu 𝑛. - 𝑛 𝐿 ≤ 𝑛 pro libovolný uzel 𝑛. - 𝑛 ≤ 𝑛 𝑅 pro libovolný uzel 𝑛. Implementaci binárního stromu si můžeme představit jako jednoduchou strukturu (c#): class KeyType : IComparable { … … } class BinTreeNode { KeyType Key; BinTreeNode Left; BinTreeNode Right; } (Interface IComparable vynucuje „srovnatelnost“ libovolných instancí typu KeyType) FI MUNI, Milan Drášil 2019 Úvod do GIS 59 Algoritmus – Vyhledání klíče v binárním stromu: 1. Vstup: uzel stromu 𝑛, klíč 𝑘. 2. Je-li 𝑛 = ∅ (strom je prázdný), potom končíme vyhledávání ”neúspěchem”. 3. Je-li 𝑛 = 𝑘, potom končíme ”úspěchem”. 4. Je-li 𝑘 < 𝑛, pokračujeme 1. pro levý podstrom 𝑛 𝐿 5. Je-li 𝑘 > 𝑛, pokračujeme 1. pro pravý podstrom 𝑛 𝑅 Algoritmus – Vkládání klíče do binárního stromu: 1. Vstup: klíč 𝑘. 2. Procházíme strom, jako bychom hledali klíč 𝑘 dokud nenarazíme na volnou pozici, tedy končíme bodem 2 předešlého algoritmu. 3. Do volné pozice vložíme klíč 𝑘. Algoritmus – Rozsahové vyhledání v binárním stromu: 1. Vstup: interval [𝑚𝑖𝑛, 𝑚𝑎𝑥], uzel stromu 𝑛. 2. Je-li 𝑛 = ∅, potom konec. 3. Patří-li 𝑛 do intervalu [𝑚𝑖𝑛, 𝑚𝑎𝑥], pošleme jej na výstup a aplikujeme algoritmus na 𝑛 𝐿 a 𝑛 𝑅. 4. Je-li 𝑚𝑎𝑥 < 𝑛, aplikujeme algoritmus na 𝑛 𝐿. 5. Je-li 𝑚𝑖𝑛 > 𝑛, aplikujeme algoritmus na 𝑛 𝑅. Povšimněme si kroků 4. a 5. Zlepšení časové složitosti spočívá v tom, že v určitých fázích algoritmů jsme schopni rozhodnout, kterou větev stromu můžeme bez rizika vynechat. Potíže s touto jednoduchou strukturou způsobuje fakt, že v jistých případech může být strom degenerovaný (např. 𝑛 𝐿 = ∅ pro všechny uzly). Degenerace nastává tehdy, když jednotlivé prvky vstupují do stromu v nevhodném pořadí (jsou uspořádány). V případě statické verze vyhledávacích problémů lze vybudovat tzv. optimální binární strom (na vstupu procedury 𝑏𝑢𝑖𝑙𝑑 známe celou množinu 𝑉). Definice - Optimální strom: Strom nazveme optimální, liší-li se počty uzlů v podstromech |𝑛 𝐿| a |𝑛 𝑅| maximálně o 1 pro jeho každý uzel 𝑛. FI MUNI, Milan Drášil 2019 Úvod do GIS 60 Poznámka: Hloubka optimálního stromu 𝑉 obsahujícího |𝑉| klíčů je: 𝑙𝑜𝑔2(|𝑉|) Algoritmus - Vybudování optimálního binárního stromu: 1. Vstup: množina klíčů 𝑉, uzel stromu 𝑛. 2. Je-li 𝑉 = ∅, skonči. 3. Rozděl množinu 𝑉 na po dvou disjunktní množiny 𝑉1, {𝑚𝑒𝑑(𝑉)}, 𝑉2 tak, že 𝑚𝑒𝑑(𝑉) je medián množiny 𝑉, klíče z 𝑉1 jsou menší než 𝑚𝑒𝑑(𝑉) a klíče z 𝑉2 jsou větší než 𝑚𝑒𝑑(𝑉). 4. Polož 𝑛 = 𝑚𝑒𝑑(𝑉). 5. Aplikuj algoritmus na množinu 𝑉1 pro levý podstrom 𝑛 𝐿. 6. Aplikuj algoritmus na množinu 𝑉2 pro pravý podstrom 𝑛 𝑅. Definice – Vyvážené stromy: Binární strom nazveme vyvážený, liší-li se hloubky 𝑛 𝐿 a 𝑛 𝑅 maximálně o 1 pro jeho každý uzel 𝑛 (hloubkou stromu rozumíme maximální délku cesty od kořene k listu). Poznámka – Rotace ve vyvážených stromech: Podmínku vyvážení lze udržovat dynamicky pomocí tzv. rotací (AVL stromy – Adelson, Velskij, Landis). Každý uzel si může zapamatovat hloubku levého a pravého podstromu. Při vložení nového klíče se strom v těch uzlech, ve kterých došlo k porušení podmínky vyvážení přeorganizujeme, např: h h+1 hhh+1h C D EA B ED B+C A++ FI MUNI, Milan Drášil 2019 Úvod do GIS 61 Počet uzlů v podstromu s kořenem 𝑛 označíme |𝑛|. Definice – 𝐵𝐵[𝛼] stromy: Buď 0 < 𝛼 < 1 2 . Binární strom patří do třídy 𝐵𝐵[𝛼] stromů, platí-li pro jeho každý uzel 𝑛: 𝛼 < |𝑛 𝐿|/(|𝑛 𝑅| + |𝑛 𝐿|) < 1 − 𝛼 Poznámka: Podmínka platí triviálně i pro pravý podstrom: 𝛼 < |𝑛 𝑅|/(|𝑛 𝑅| + |𝑛 𝐿|) < 1 − 𝛼 Poznámka – význam parametru 𝛼: Pokud byl v nějakém okamžiku podstrom definovaný uzlem 𝑛 optimální, pak k porušení podmínky z definice 𝐵𝐵[𝛼] stromů musí dojít k minimálně 𝑐 × |𝑛| vložení/mazání uzlů do/z příslušného podstromu (𝑐 je konstanta závislá pouze na parametru 𝛼). To znamená, že čím více má podstrom uzlů, tím méně častěji dochází k porušení podmínky definice 𝐵𝐵[𝛼] a nutnosti jeho reorganizace. Definice – 𝑘𝐷 strom: Úrovní 𝑙𝑒𝑣𝑒𝑙(𝑛) uzlu 𝑛 binárního stromu rozumíme délku cesty k tomuto uzlu od kořene stromu. Buď (𝑆, ≤) úplně uspořádaná množina, 𝑘 > 0, 𝑥 = (𝑥0, . . , 𝑥𝑖, . . , 𝑥 𝑘−1), 𝑦 = (𝑦0, . . , 𝑦𝑖, . . , 𝑦 𝑘−1) ∈ 𝑆 𝑘 Řekneme, že 𝑥 ≤𝑖 𝑦, jestliže 𝑥𝑖 ≤ 𝑦𝑖 a pro 𝑎 ∈ 𝑆 𝑥 ≤𝑖 𝑎 jestliže 𝑥𝑖 ≤ 𝑎 Analogicky jsou definovány operátory <𝑖, ≥𝑖, >𝑖 𝑘𝐷 stromem nad S nazveme binární strom, jehož uzly jsou k-tice z 𝑆 𝑘 , a kde pro každý uzel 𝑛, jeho levý podstrom 𝑛 𝐿 resp. 𝑛 𝑅 platí: 𝑛𝑙 ∈ 𝑛 𝐿 ⇒ 𝑛𝑙 ≤𝑖 𝑛, 𝑛 𝑟 ∈ 𝑛 𝑅 ⇒ 𝑛 ≤𝑖 𝑛 𝑟 kde 𝑖 = 𝑙𝑒𝑣𝑒𝑙(𝑛) Algoritmus – Vyhledání bodu ve 2-d stromu: Analogicky k binárním stromům, s tím rozdílem, že na každé úrovni použijeme jinou srovnávací metodu (≤𝑖). FI MUNI, Milan Drášil 2019 Úvod do GIS 62 Algoritmus – Vložení bodu do 2-d stromu: Analogicky k binárním stromům, hledáme ve stromu „bod“ dokud nenarazíme na volnou pozici. Do ní vložíme nový klíč. Geometrická interpretace 2𝐷 stromu, každý vložený bod dělí, jistou část roviny na dvě „poloroviny“. 4 3 2 1 4 23 1 Algoritmus – Rozsahový dotaz pro body ve 2𝐷 stromu: 1. Vstup: Uzel 𝑛 a obdélník [𝑥𝑚𝑖𝑛, 𝑦𝑚𝑖𝑛, 𝑥𝑚𝑎𝑥, 𝑦𝑚𝑎𝑥]. 2. Je-li 𝑛 = ∅, konec. 3. Je-li 𝑛 (tj. „bod“ ve 2D prostoru) uvnitř dotazového obdélníku, pošli jej na výstup a aplikuj algoritmus na oba dva syny. 4. Vyber syny, pro které budeš aplikovat algoritmus, a to podle úrovně ve které se nachází uzel 𝑛, například je-li: 𝑙𝑒𝑣𝑒𝑙(𝑛) 𝑚𝑜𝑑 2 = 0 ∧ 𝑛 >0 𝑥𝑚𝑎𝑥 potom aplikuj algoritmus jen na větev 𝑛 𝐿, analogicky pro další možné případy. 4 3 2 1 Dotazový obdélník FI MUNI, Milan Drášil 2019 Úvod do GIS 63 Vyvažování multidimensionálních stromů je komplikované. Nedají se totiž provádět rotace jako v klasických binárních stromech, protože v každém patře stromu měníme srovnávací kritérium, např. (obrázek rotace) z 𝐵 >0 𝐴  𝐷 <1 𝐵 neplyne 𝐴 <0 𝐵  𝐷 >1 𝐴 Pomocí 𝐵𝐵[𝛼] techniky lze však 𝑘𝐷 stromy udržovat vyvážené pomocí tzv. částečné reorganizace, tedy „hlídat“ v každém uzlu 𝑘𝐷 stromu poměr počtu uzlů v jeho levém a pravém podstromu. V případě porušení podmínky 𝐵𝐵[𝛼] nahradit optimální podstrom. Algoritmus – Vybudování optimálního 2𝐷 stromu: 1. Vstup: množina bodů 𝑉, uzel stromu 𝑛. 2. Je-li 𝑉 = ∅, skonči. 3. 𝑙 = 𝑙𝑒𝑣𝑒𝑙(𝑛), rozděl množinu 𝑉 na po dvou disjunktní množiny 𝑉1, {𝑚𝑒𝑑𝑙(𝑉)}, 𝑉2 tak, že 𝑚𝑒𝑑𝑙(𝑉) je takový bod, že jeho 𝑙-tá souřadnice je medián množiny 𝑙–tých souřadnic z 𝑉, 𝑙-té souřadnice z 𝑉1 jsou menší než 𝑚𝑒𝑑𝑙(𝑉) a 𝑙-té souřadnice z 𝑉2 jsou větší než 𝑚𝑒𝑑𝑙(𝑉). 4. Polož 𝑛 = 𝑚𝑒𝑑𝑙(𝑉). 5. Aplikuj algoritmus na množinu 𝑉1 pro 𝑛 𝐿 a 𝑉2 pro 𝑛 𝑅. Poznámka: Rozdělení množiny z kroku 3. lze realizovat snadnou modifikací algoritmu 𝑄𝑢𝑖𝑐𝑘𝑆𝑜𝑟𝑡. Metodu 𝑘𝐷 stromů lze použít i na obdélníky, které můžeme považovat za 4𝐷 body [𝑥𝑚𝑖𝑛, 𝑦𝑚𝑖𝑛, 𝑥𝑚𝑎𝑥, 𝑦𝑚𝑎𝑥]. Použijeme tedy 4𝐷 strom. Algoritmus – Rozsahový výběr pro obdélníky ve 4𝐷 stromu: 1. Vstup: kořen stromu 𝑛, dotazový obdélník [𝑥𝑚𝑖𝑛, 𝑦𝑚𝑖𝑛, 𝑥𝑚𝑎𝑥, 𝑦𝑚𝑎𝑥]. 2. Je-li 𝑛 = ∅, skonči. 3. Je-li obdélník 𝑛 incidentní s dotazem, pošli 𝑛 na výstup a aplikuj algoritmus na 𝑛 𝐿 a 𝑛 𝑅. 4. Jinak, podle úrovně 𝑙𝑒𝑣𝑒𝑙(𝑛), ve které se nacházíš ve stromu, se rozhodni, zda můžeš vynechat nějakou větev. Je-li např.: 𝑙𝑒𝑣𝑒𝑙(𝑛𝑜𝑑𝑒) 𝑚𝑜𝑑 4 = 0 ∧ 𝑛 >0 𝑥𝑚𝑎𝑥 vynechej 𝑛 𝑅. Analogicky pro další úrovně, v každé se dá za jistých podmínek jedna větev vynechat. FI MUNI, Milan Drášil 2019 Úvod do GIS 64 5. Aplikuj 1. pro vybrané následníky. 4.4. QUAD TREE – KVADRANTOVÉ STROMY Kvadrantové stromy jsou prostorové vyhledávací struktury založené na pravidelném dělení části 2D prostoru (většinou čtverce) na čtyři stejné části. Každý uzel reprezentuje (čtvercovou) část oblasti a je opatřen příznakem, zda se areálový jev v této oblasti vyskytuje. Trik spočívá v tom, že je-li uzel označen jako prázdný/plný potom nemusíme budovat prázdný/plný podstrom – všechny jeho uzly „dědí“ tuto vlastnost. V základní verzi tedy QuadTree mohou přímo definovat polygonální geometrii. 𝑧– 𝑖– 𝑗 na obrázku označuje: 𝑧 – hloubka uzlu 𝑖, 𝑗 maticové indexy vzhledem všem uzlům v dané úrovni. FI MUNI, Milan Drášil 2019 Úvod do GIS 65 Definice – kvadrantový strom: Kvadrantový strom hloubky ℎ je kvartérní strom (tj. každý uzel má maximálně 4 následníky) s těmito vlastnostmi: - Uzel 𝑛 je tvořen obdélníkem 𝑅 = [𝑥𝑚𝑖𝑛, 𝑦𝑚𝑖𝑛, 𝑥𝑚𝑎𝑥, 𝑦𝑚𝑎𝑥] a příznakem 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛) ∈ {𝑒𝑚𝑝𝑡𝑦, 𝑓𝑢𝑙𝑙, ℎ𝑎𝑙𝑓}. - List 𝑛 má příznak 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛) ∈ {𝑒𝑚𝑝𝑡𝑦, 𝑓𝑢𝑙𝑙}. - Nelistový uzel má čtyři následníky a příznak 𝑐𝑜𝑛𝑡𝑒𝑛𝑡 = ℎ𝑎𝑙𝑓. - Následníci uzlu dělí obdélník rodičovského uzlu na čtyři „stejné“ části (podle středu). - Maximální délka od kořene k listu je ℎ. Podle toho, jakou část rodičovského obdélníku reprezentuje uzel, označíme jej 𝑇𝐿, 𝑇𝑅, 𝐵𝐿, 𝐵𝑅 (Top/Bottom Left/Right). Poznámka: Struktura je velmi dobře použitelná v těch případech, kdy je prostor pevně rozdělen na „dlaždice“, například pro služby typu WMTS (viz datové sklady GIS). Každá úroveň podrobnosti 𝑛 je reprezentována maticí dlaždic 2 𝑛 × 2 𝑛 .Struktura potom může odpovídat na dotaz, zda je dlaždice obsazená (je na ní mapa/jev). Strukturu lze implementovat jednoduchým objektem: public class QtreeNode { public QtreeNode TopLeft; public QtreeNode TopRight; public QtreeNode BottomLeft; public QtreeNode BottomRight; public QtreeContent Content;} Struktura nemusí obsahovat souřadnice čtverce/obdélníku, který reprezentuje, stačí si pamatovat nultý čtverec/obdélník kořene stromu. Souřadnice ostatních uzlů jsou potom definovány cestou od kořene. Poznámka: Vlastnost ℎ𝑎𝑙𝑓 je redundantní. Je ekvivalentní s „uzel není list“. Je zavedena pouze pro přehlednost. FI MUNI, Milan Drášil 2019 Úvod do GIS 66 Algoritmus – vyhledání dlaždice v QuadTree: 1. Vstup sloupec 𝑐𝑜𝑙, řádek 𝑟𝑜𝑤, podrobnost 𝑧𝑜𝑜𝑚 a uzel 𝑛. Nechť 𝑐𝑖, 𝑖 = 1, . . , 𝑧𝑜𝑜𝑚 je binární reprezentace sloupcového indexu 𝑐𝑜𝑙 a 𝑟𝑖, 𝑖 = 1, . . , 𝑧𝑜𝑜𝑚 je binární reprezentace řádkového indexu 𝑟𝑜𝑤. 2. 𝑖 = 0 3. Je-li 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛) ≠ ℎ𝑎𝑙𝑓 vrať 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛) a konec. 4. Je-li 𝑖 = 𝑧𝑜𝑜𝑚 vrať 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛) a konec. 5. Vyber pokračovací uzel 𝑛𝑒𝑥𝑡 takto: 𝑟𝑖 = 0 - nahoru 𝑟𝑖 = 1 - dolů, 𝑐𝑖 = 0 - doleva, 𝑐𝑖 = 1 – doprava 6. 𝑛 = 𝑛𝑒𝑥𝑡, 𝑖 = 𝑖 + 1 a pokračuj 3. Na struktuře lze velmi snadno a elegantně provádět množinové operace. Stačí implementovat operaci inverze a např. průniku, ostatní množinové operace lze provést pomocí těchto dvou. Algoritmus – inverze QuadTree 1. Vstup uzel 𝑛. 2. Je-li 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛) = 𝑓𝑢𝑙𝑙 potom 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛): = 𝑒𝑚𝑝𝑡𝑦 a návrat. 3. Je-li 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛) = 𝑒𝑚𝑝𝑡𝑦 potom 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛): = 𝑓𝑢𝑙𝑙 a návrat. 4. Je-li 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛) = ℎ𝑎𝑙𝑓 potom proveď kroky 2. a 3. pro všechny následníky. Algoritmus – průnik QuadTree: 1. Vstup uzly dvou stromů 𝑛1, 𝑛2 (reprezentující stejnou plochu). 2. Je-li 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛1) = 𝑒𝑚𝑝𝑡𝑦 ∨ 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛2) = 𝑒𝑚𝑝𝑡𝑦, potom je výsledek uzel 𝑛 s obsahem 𝑒𝑚𝑝𝑡𝑦. 3. Je-li 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛1) = 𝑓𝑢𝑙𝑙, potom je výsledek uzel 𝑛2. 4. Je-li 𝑐𝑜𝑛𝑡𝑒𝑛𝑡(𝑛2) = 𝑓𝑢𝑙𝑙, potom je výsledek uzel 𝑛1. 5. Jinak je výsledek nový uzel 𝑛 s obsahem ℎ𝑎𝑙𝑓 a následníky z kroků 1. – 4. FI MUNI, Milan Drášil 2019 Úvod do GIS 67 Pevný kvartérní strom (non-pointer Quad Tree) Zájmové území je postupně děleno na obdélníkové části a podle nich je jim přidělován „klíč“ 4300 4400 42003000 2000 4100 1000 Obr. - Číslování obdélníků-dlaždic v non-pointer Quad Tree. Obdélník bude mít index takové dlaždice „pevné struktury“ která je jeho nadmnožinou a je nejmenší s touto vlastností. - Pro libovolný obdélník 𝑅 označme 𝑄(𝑅) jeho klíč v non-pointer QuadTree. - Pro libovolný klíč 𝐾 označme jeho „nenulovou“ část, tedy levý podřetězec symbolem 𝑁𝑍(𝐾). - Délku znakového řetězce 𝐾 označme 𝑙𝑒𝑛(𝐾). - Podřetězec řetězce 𝐾 z levé strany délky 𝑙 označme 𝑠𝑢𝑏(𝐾, 𝑙). Tvrzení – incidence obdélníků n non-pointer QuadTree: Buďte 𝐴, 𝐵 libovolné obdélníky, jejichž strany jsou rovnoběžné s osami souřadného systému. Nechť dále 𝐴 ∩ 𝐵 ≠ ∅. Označíme-li, 𝑙 = 𝑚𝑖𝑛{𝑙𝑒𝑛(𝑁𝑍(𝑄(𝐴))), 𝑙𝑒𝑛(𝑁𝑍(𝑄(𝐵)))} potom: 𝑠𝑢𝑏(𝑄(𝐴), 𝑙) = 𝑠𝑢𝑏(𝑄(𝐵), 𝑙) tj. klíče 𝐴 a 𝐵 mají stejnou společnou nenulovou část. FI MUNI, Milan Drášil 2019 Úvod do GIS 68 Algoritmus - Vyhledání obdélníků v non-pointer QuadTree: 1. Vstup – obdélník 𝑆 = [𝑥𝑚𝑖𝑛, 𝑦𝑚𝑖𝑛, 𝑥𝑚𝑎𝑥, 𝑦𝑚𝑎𝑥]. 2. Pošli na výstup všechny obdélníky 𝐴, pro které: 𝑠𝑢𝑏 (𝑄(𝐴), 𝑙𝑒𝑛 (𝑁𝑍(𝑄(𝑆)))) = 𝑁𝑍(𝑄(𝑆)) ∧ 𝐴 ∩ 𝑆 ≠ ∅ 3. Pošli na výstup všechny obdélníky 𝐴, pro které: 𝑄(𝐴) = 𝑃 ∧ 𝐴 ∩ 𝑆 ≠ ∅ kde 𝑃 jsou všechny klíče, které jsou na cestě od 𝑄(𝑆) ke kořenu, tj. v 𝑄(𝑆) zprava postupně nahrazujeme nenulové číslice nulami. Tento postup má jednu nevýhodu. V případě, že například dotaz inciduje se středem území, potom procházíme v bodě 2 všechno. Této nevýhodě se vyhneme dekomponování dotazu. Algoritmus - Dekompozice dotazu v non-pointer QuadTree: 1. Vstup – dotazový obdélník 𝑆. 2. Rozděl obdélník 𝑆 na obdélníky 𝑆1 a 𝑆2 (𝑆1 ∪ 𝑆2 = 𝑆) podle takové souřadnice 𝑥 resp. 𝑦, která způsobila klíčování v non-pointer quadTree, tj. takovou, která ohraničuje nějaký čtverec v non-pointer quadTree a prochází dotazovým obdélníkem 𝑆. V případě, že taková souřadnice neexistuje potom obdélník 𝑆 neděl a konec. 3. Aplikuj krok 2. na čtverce 𝑆1 a 𝑆2 podle druhé souřadnice. Tímto postupem získáme maximálně 4 obdélníky na které aplikujeme algoritmus vyhledání obdélníků Výhody této metody: - Velmi snadná implementace v prostředí SQL – tedy relačních databází. Například norma OGC WKB nepředepisuje metodu efektivního výběru, tímto způsobem ji můžeme doplnit. - „jeden objekt“ reprezentuje „jeden klíč“ znamená, že prostorová indexace je zabezpečena přímo v geometrické tabulce. Prostorový výběr nevyžaduje součin, či spojení s dalšími tabulkami. FI MUNI, Milan Drášil 2019 Úvod do GIS 69 Dekompozice dotazu – 4 obdélníky s klíči: 2330000000,2343000000,4110000000,4120000000 SELECT ID FROM KM_ALL WHERE ( (SPAT_KEY BETWEEN '2330000000' AND '2335000000') OR (SPAT_KEY BETWEEN '2343000000' AND '2343500000') OR (SPAT_KEY BETWEEN '4110000000' AND '4115000000') OR (SPAT_KEY BETWEEN '4120000000' AND '4125000000') ) OR SPAT_KEY IN ('0000000000', '2000000000', '2300000000', '2340000000', '4000000000', '4100000000' ) AND (xmax>=-642646042) AND (ymax>=-1114990337) AND (xmin<=-569087654) AND (ymin<=-1070777051) 4.5. R-STROMY Definice – B+-stromy: B-strom řádu 𝑚 je strom s těmito vlastnostmi: - Každý uzel má maximálně 𝑚 synů. - Každý uzel, s výjimkou kořene a listů, má minimálně 𝑚/2 synů. - Kořen má minimálně 2 syny, pokud není list. - Všechny listy jsou na stejné úrovni. - Nelistový uzel s 𝑘 syny obsahuje 𝑘 − 1 klíčů - Klíče v uzlu 𝑘𝑒𝑦1, . . , 𝑘𝑒𝑦 𝑘 jsou vzestupně uspořádány - Uzly B+ stromu mají tvar 𝑝0 𝑘𝑒𝑦1 𝑝1 … 𝑝 𝑘−1 𝑘𝑒𝑦 𝑘 𝑝 𝑘. - ukazatel 𝑝𝑖 ukazuje na uzel, jehož všechny klíče jsou v intervalu [𝑘𝑒𝑦𝑖, 𝑘𝑒𝑦𝑖+1], formálně předpokládáme, že 𝑘𝑒𝑦0 = −∞ a 𝑘𝑒𝑦 𝑘+1 = ∞ FI MUNI, Milan Drášil 2019 Úvod do GIS 70 Algoritmus vložení klíče do B+ stromu: 1. Vstup klíč 𝑘𝑒𝑦 a kořen B+ stromu. 2. Není-li uzel list vyber dvojici klíčů 𝑘𝑒𝑦𝑖, 𝑘𝑒𝑦𝑖+1 s vlastností 𝑘𝑒𝑦𝑖−1 < 𝑘𝑒𝑦 < 𝑘𝑒𝑦𝑖 a pokračuj uzlem na který ukazuje 𝑝𝑖. Opakuj tento krok. Dokud uzel není list. 3. Vlož uzel do listu na správnou pozici. Je-li počet uzlů větší než 𝑚, rozděl uzel. Algoritmus dělení uzlů v B+ stromu řádu 𝑚 : 1. Vstup uzel s 𝑚 + 1 klíči. 2. Nechť 𝑘 = 𝑚/2 + 1, vytvoř dva nové uzly: 𝑝0 𝑘𝑒𝑦1 𝑝1 … 𝑝 𝑘−2 𝑘𝑒𝑦 𝑘−1 𝑝 𝑘−1 a 𝑝 𝑘 𝑘𝑒𝑦 𝑘+1 𝑝 𝑘+1 … 𝑝 𝑚 𝑘𝑒𝑦 𝑚+1 𝑝 𝑚+1 a ukazatele na ně 𝑞 resp. 𝑟 3. Je-li uzel kořen, vytvoř nový prázdný kořen: 𝑞 𝑘𝑒𝑦 𝑘 𝑟 Jinak, nechť 𝑝𝑖 je ukazatel z otcovského uzlu. Vlož klíč 𝑘𝑒𝑦 𝑘 do otcovského uzlu na pozici: . . 𝑘𝑒𝑦𝑖 𝑝𝑖 𝑘𝑒𝑦𝑖+1 → 𝑘𝑒𝑦𝑖 𝑞 𝑘𝑒𝑦 𝑘 𝑟 𝑘𝑒𝑦𝑖+1 4. Má-li otcovský uzel 𝑚 + 1 klíčů, opakuj pro něj kroky 1.-3. Definice R-strom: Analogie k B-stromům, klíče jsou obdélníky. 𝑀 – maximální počet klíčů v uzlu, 𝑚 ≤ 𝑀/2 – minimální počet klíčů v uzlu - Každý uzel obsahuje minimálně 𝑚 klíčů a maximálně 𝑀 klíčů, pokud není kořen. - Klíče v R-stromech jsou obdélníky s ukazateli na synovské uzly, v listech obdélníky s ukazateli na geometrické prvky. Následníka obdélníku 𝐾 označíme 𝑛𝑒𝑥𝑡(𝐾), rodičovský uzel uzlu 𝑛 označíme 𝑝𝑟𝑒𝑣(𝑛). - Pro synovské uzly platí, že jejich klíče (tj. obdélníky) jsou uvnitř “otcovského” obdélníku. FI MUNI, Milan Drášil 2019 Úvod do GIS 71 - Listy stromu jsou na téže úrovni. - Kořen obsahuje minimálně dva klíče, pokud není list. Algoritmus – Vyhledání klíčů obdélníků v R- stromech: 1. Vstup uzel R-stromu 𝑛. Dotazový obdélník 𝑄. 2. Je-li 𝑛 list, potom všechny klíče z 𝑛 incidentní s 𝑄 na výstup. 3. Jinak aplikuj algoritmus na syny takových klíčů z 𝑛, pro které je klíč incidentní s 𝑄. Algoritmus – Vkládání klíčů do R- stromů: 1. Vstup, klíč 𝐾 (obdélník) 2. Vyhledej list 𝑛: a. Polož 𝑛: = 𝑟𝑜𝑜𝑡 (kořen stromu). b. Je-li 𝑛 list pokračuj 3., jinak c). c. Nechť 𝑘 je klíč v 𝑛, jehož obdélník vyžaduje nejmenší rozšíření takové, aby obsahoval 𝐾. Rozšiř 𝑘 o klíč 𝐾, polož 𝑛: = 𝑛𝑒𝑥𝑡(𝑛) a pokračuj b. 3. Přidej 𝑘 do vybraného listu 𝑛. Q1 Q2 R1 R2 S1 S2 S3 R3 2 2 4 6 8 10 12 14 16 4 6 8 R1 R2 S1 S2 S3 R3 Q1 Q2 FI MUNI, Milan Drášil 2019 Úvod do GIS 72 4. Je-li počet klíčů v 𝑛 menší, nebo roven 𝑀 konec. Jinak rozděl uzel 𝑛 na dva nové uzly. Je-li 𝑛 kořen, vytvoř nový kořen se dvěma novými klíči, jinak odstraň z rodičovského uzlu původní klíč a nahraď jej dvěma novými klíči a polož 𝑛 = 𝑝𝑟𝑒𝑣(𝑛). 5. Opakuj 4. Problém rozdělení uzlu v R-Tree: Najdi dva obdélníky (možná incidentní) s následujícími vlastnostmi (NP – úplný problém): - Sjednocení obou obdélníků je původní obdélník - Oba obdélníky obsahují zhruba stejný počet klíčů, splňují podmínku z definice R-tree - Oba obdélníky se překrývají co nejméně Algoritmus dělení uzlu R-stromu (kvadratická složitost): 1. Vyber první dva obdélníky a) Pro každou dvojici klíčů 𝑘, 𝑙 vytvoř minimální obdélník 𝑗 obsahující oba klíče a polož: 𝑝(𝑘, 𝑙) = 𝑃𝑙𝑜𝑐ℎ𝑎(𝑗) − 𝑃𝑙𝑜𝑐ℎ𝑎(𝑘) − 𝑃𝑙𝑜𝑐ℎ𝑎(𝑙) b) Vyber dvojici obdélníků 𝑘, 𝑙 s maximem 𝑝(𝑘, 𝑙), zařaď je do první a druhé skupiny. 2. V případě, že jedna skupina obsahuje tak málo obdélníků, že pro zachování podmínky minima 𝑚 musí obsahovat všechny nezařazené obdélníky, zařaď do ní zbývající obdélníky a konec. 3. Pro všechny nezařazené obdélníky spočítej rozdíl ploch, o které se zvětší obdélníky první a druhé skupiny začleněním nezařazeného obdélníku. 4. Vyber obdélník z 3. který má maximální rozdíl ploch a zařaď ho do té skupiny, jejíž celkový obdélník se rozšíří méně. Pokračuj krokem 2. FI MUNI, Milan Drášil 2019 Úvod do GIS 73 Algoritmus dělení uzlu R-stromu (lineární složitost): 1. Vyber první dva obdélníky: a. Pro každou dimenzi najdi klíče s maximem minima a minimem maxima, stanov „separační vzdálenost“ mezi těmito klíči (minimum minus maximum). b. Normalizuj separační vzdálenost tak, že vzdálenost intervalů podělíš rozsahem všech klíčů v dané dimenzi. c. Vyber dvojici 𝑘, 𝑙 s největší normalizovanou separační vzdáleností, zařaď je do první a druhé skupiny. 2. V případě, že jedna skupina obsahuje tak málo obdélníků, že pro zachování podmínky minima 𝑚 musí obsahovat všechny nezařazené obdélníky, zařaď do ní zbývající obdélníky a konec. 3. Vezmi další nezařazený klíč a zařaď jej do takové skupiny, jejíž MBR vyžaduje menší rozšíření. R-Tree je nejpoužívanější metoda prostorové indexace, je nezávislá na velikosti objektů, nemusíme znát rozsah území, poměrně snadno je modifikovatelná na polygonální dotazy (viz algoritmus dotazu, dotazový obdélník můžeme nahradit dotazovým polygonem). FI MUNI, Milan Drášil 2019 Úvod do GIS 74 5. Funkce a operace nad geometrickými objekty 5.1. KONVERZNÍ FUNKCE (OGC): Funkce umožňující konverze formátů, například vkládání do RDBMS obecným klientem SQL. AsText(g Geometry) : String AsBinary(g Geometry) : Binary Používají pro vkládání a výběr geometrických objektů. 5.2. MĚŘÍCÍ FUNKCE Délka: Length(l LineString) : double Length(l Polygon) : double Plocha areálu: Předpokládáme, že polygony (hranice) jsou “správně” orientovány, Area(l Polygon) : Double Precision 𝐴 = 1/2 ∑(𝑥𝑖 − 𝑥𝑖+1)(𝑦𝑖 + 𝑦𝑖+1) 𝑖 5.3. POLOHOVÉ FUNKCE Určení orientace polygonu: Algoritmus - Určení orientace polygonu: 1. Vyber z hranic oblastí takovou hranici a tři po sobě jdoucí její body tak, aby střední bod měl minimální souřadnici y (ze všech souřadnic 𝑦 v polygonu) a první bod měl souřadnici 𝑦 větší než bod prostřední. 2. Rotuj souřadnou soustavu tak, aby orientovaná úsečka definovaná prvními dvěma body splynula s osou 𝑥 v kladném směru. 3. Znaménko souřadnice y posledního bodu určuje orientaci polygonu. FI MUNI, Milan Drášil 2019 Úvod do GIS 75 y>0 a a b b (0,0) Poloha bodu vůči polygonu: Algoritmus - Bod v polygonu: 1. Najdi v bodech polygonu bod, jehož 𝑦 souřadnice je různá od souřadnice bodu, který testujeme. Nechť 𝑝𝑝 je polopřímka vycházející z testovaného bodu rovnoběžná s osou 𝑥 v kladném směru. 𝑛𝑃𝑟𝑢𝑠: = 0. Od vybraného bodu postupně procházej všechny úsečky a proveď body 2 – 5. 2. Leží-li testovaný bod na úsečce, ukonči proceduru s výsledkem 𝑁𝐴_𝐻𝑅𝐴𝑁𝐼𝐶𝐼. 3. Má-li úsečka vlastní průsečík s polopřímkou 𝑝𝑝, potom 𝑛𝑃𝑟𝑢𝑠: = 𝑛𝑃𝑟𝑢𝑠 + 1. 4. Končí-li úsečka na polopřímce a začíná-li mimo polopřímku, stanov podle počátku úsečky 𝑜𝑑𝑘𝑢𝑑: = 𝑃𝑂𝐷 nebo 𝑜𝑑𝑘𝑢𝑑: = 𝑁𝐴𝐷. 5. Začíná-li úsečka na polopřímce a končí-li mimo polopřímku a pokračuje-li do jiné poloroviny, než je stav proměnné 𝑜𝑑𝑘𝑢𝑑, potom 𝑛𝑃𝑟𝑢𝑠: = 𝑛𝑃𝑟𝑢𝑠 + 1. 6. Je-li 𝑛𝑃𝑟𝑢𝑠 sudý, ukonči proceduru s výsledkem 𝑉𝑁Ě. 7. Je-li 𝑛𝑃𝑟𝑢𝑠 lichý, ukonči proceduru s výsledkem 𝑈𝑉𝑁𝐼𝑇Ř. FI MUNI, Milan Drášil 2019 Úvod do GIS 76 Poloha bodu vůči úsečce: (0,0) Rotujeme souřadnou soustavu tak, aby její počátek splynul s počátečním bodem úsečky a koncový bod ležel na x-ové souřadnici v kladném směrů. Určení polohy je potom triviální operace (nalevo, napravo, minimální vzdálenost…) 𝑥′ = 𝑥. 𝑐𝑜𝑠(𝛼) − 𝑦. 𝑠𝑖𝑛(𝛼) 𝑦′ = 𝑥. 𝑠𝑖𝑛(𝛼) + 𝑦. 𝑐𝑜𝑠(𝛼) Vzdálenost geometrických objektů: Distance(g1 Geometry, g2 Geometry) : Double Kombinace metod podle typu geometrií: Vzdálenost bod-bod, poloha bodu vůči úsečce, průsečík úseček, bod v polygonu … 5.4. GEOMETRICKÉ OPERÁTORY Konvexní obal množiny bodů: Je nejmenší konvexní polygon s takovou vlastností, že všechny vstupní body leží uvnitř něj nebo na jeho hranici. FI MUNI, Milan Drášil 2019 Úvod do GIS 77 Konvexní polygon je polygon s touto vlastností: Každý vnitřní bod úsečky, jejíž krajní body leží na hranici polygonu, je i vnitřním bodem polygonu. Nechť polygon 𝐶 = [𝑠0, . . , 𝑠 𝑛] je konvexní obal množiny bodů 𝑆, 𝑠0 má minimální 𝑦 souřadnici a maximální 𝑥 souřadnici (v případě více minim 𝑦). Nechť dále [𝛼1, . . , 𝛼 𝑛] jsou úhly mezi 𝑥 osou a úsečkou [𝑠0, 𝑠 𝑛]. Potom: - Bod 𝑠𝑖 leží nalevo od přímky definované body [𝑠𝑖−2, 𝑠𝑖−1], (úhel [𝑠𝑖−2, 𝑠𝑖−1, 𝑠𝑖] je konvexní). - Polygon 𝐶 je tvořen nejvíce body z 𝑆 s předešlou vlastností - 𝛼𝑖  𝛼𝑖+1. Algoritmus – konvexní obal (Graham scan - 𝑂(𝑁𝑙𝑜𝑔(𝑁)): 1. Vstup – soubor 𝑁 bodů 𝑆. 2. Vyber z 𝑆 bod 𝑃0 s minimální 𝑦-souřadnicí, ten nejvíce vpravo (max. 𝑥- souřadnice). 3. Setřiď body podle úhlu, které svírá přímka procházející daným bodem a bodem 𝑃0 do pole 𝑃[], bod 𝑃0 bude prvním bodem pole. 4. Vlož do zásobníku 𝑅 bod 𝑃[0] = 𝑃0 a bod 𝑃[1]. 5. Pro 1 < 𝑖 < 𝑁, kde 𝑁 je počet bodů v 𝑃: 6. Nechť 𝑝1 je první bod v zásobníku 𝑅, 𝑝2 druhý bod 7. Je-li 𝑃[𝑖] napravo od přímky 𝑝1 → 𝑝2 vlož 𝑃[𝑖] do zásobníku a 𝑖: = 𝑖 + 1, jinak odstraň 𝑝1 ze zásobníku a znovu 6. 8. Zásobník 𝑅 obsahuje konvexní obal bodů. FI MUNI, Milan Drášil 2019 Úvod do GIS 78 Triangulace polygonu: Je postup, kterým získáme sadu trojúhelníků, které pokrývají vstupní polygon. Algoritmus triangulace „Ear clipping“: (David Eberly - http://www.geometrictools.com/) 1. Vstup polygon bez děr. Vytvoř prázdný seznam trojúhelníků 𝑇 a pokračuj 2. 2. Vstup polygon bez děr, seznam trojúhelníků 𝑇. 3. Zorientuj polygon tak, aby měl kladnou orientaci. 4. Je-li polygon trojúhelník, přidej ho do seznamu a konec – výstup 𝑇. 5. Procházej postupně trojice po sobě jdoucích lomových bodů 𝑝𝑖, 𝑝𝑖+1, 𝑝𝑖+2. Zjisti, zda je úhel bodů 𝑝𝑖, 𝑝𝑖+1, 𝑝𝑖+2 konvexní, použijeme algoritmus polohy bodu 𝑝𝑖+2 vůči úsečce, [𝑝𝑖, 𝑝𝑖+1]), po posunu a rotaci musí mít poslední bod kladnou 𝑦 souřadnici. Není-li tomu tak, pokračuj další trojicí bodů. 6. Otestuj, zda pro trojúhelník 𝑡 tvořený zkoumanou trojicí platí – bod 𝑝𝑖 je „viditelný“ z bodu 𝑝𝑖+2): a) Žádný lomový bod vstupního polygonu neleží uvnitř 𝑡. b) Žádná úsečka z hranice vstupního polygonu neprotíná žádnou úsečku hranice 𝑡. 7. V případě 5. je splněno, přidej do seznamu 𝑇 trojúhelník 𝑡. Odstraň bod 𝑝𝑖+1 z hranice polygonu a pokračuj 2. (snížili jsme počet bodů v hranici a rekurzivně aplikujeme týž postup). FI MUNI, Milan Drášil 2019 Úvod do GIS 79 Algoritmus „Ear clipping“ – pro polygony s dírami Postupujeme tak, že polygon zorientujeme, vytvoříme z hranic jedinou hranici vhodným zřetězením a aplikujeme předešlý algoritmus. Například: Vhodným zřetězením hranice rozumíme takový postup, že můžeme vkládat ‚spojnice‘ mezi jednotlivé hranice pouze mezi body, které jsou vzájemně přímo viditelné, tj. novou spojnici neprotíná žádná úsečka z řetězených hranic polygonu. Poznámka: Pro velké polygony vytvoříme prostorový index na všechny úsečky hranic a lomové body. Indexových struktur využijeme v kroku 6. Poznámka: Zřetězením hranic polygony s dírami obdržíme ‚nekorektní‘ hranici, je tečná sama se sebou, to ale pro navazující algoritmus není podstatné, vlivem zorientování zůstává zachována základní vlastnost konvexity úhlů. Množinové operace: Geometry Intersection (g1 Geometry, g2 Geometry) Geometry Difference (g1 Geometry,g2 Geometry) Geometry Union (g1 Geometry, g2 Geometry) Geometry SymDifference(g1 Geometry, g2 Geometry) Geometry Buffer (g1 Geometry, d Double) Bod – bod: Triviální operace. Pozor, pro příslušnost bodu k množině je nutné použít vhodnou přístupovou metodu k prostorovým datům. FI MUNI, Milan Drášil 2019 Úvod do GIS 80 Bod – lomená čára: Vzájemná poloha úsečka X bod. Bod – oblast: Poloha bodu vůči polygonu Lomená čára – lomená čára: Poloha dvou úseček. Lomená čára – oblast: Algoritmus – Průnik lomené čáry s oblastí. 1. Vstup: oblast a lomená čára. 2. Ze vstupní lomené čáry vytvoř seznam 𝑃 segmentů lomené čáry takových, které buď neprotínají hranice oblasti, nebo jsou celé tečné. 3. Ze seznamu 𝑃 vytvoř seznam 𝑆 ⊆ 𝑃 takový, že libovolný vnitřní bod každého segmentu z 𝑆 leží uvnitř oblasti. 4. Zřetěz segmenty z 𝑆 do “co nejdelších” lomených čar, a výsledek pošli na výstup Oblast – oblast: Množinové operace. Pro ty stačí algoritmizovat např. operaci doplňku a průniku, ostatní lze jimi vyjádřit. Doplněk: Prostá změna orientace hranic oblasti. Průnik dvou oblastí: FI MUNI, Milan Drášil 2019 Úvod do GIS 81 Předpokladem je, že hranice vstupních polygonů jsou správně orientovány. Tyto doplníme o vnitřní body tak, že mimo lomové body se hranice vstupních polygonů neprotínají. Leží-li část hranice oblasti uvnitř oblasti druhé, potom je i hranicí průniku. V jejím levém okolí se totiž nachází vnitřní body první oblasti a v pravém nikoli, zatímco v jejím levém i pravém okolí se nachází vnitřní body oblasti druhé. Část hranice je tedy součástí hranice průniku. Obdobnou úvahou odvodíme, že pokud část oblasti splývá s částí oblasti druhé, potom je tato část hranicí operace průniku právě když splývající části mají totožnou orientaci. Algoritmus – Průnik dvou oblastí: 1. Vstup: dvě orientované oblasti. 2. Všechny hrany hranic oblastí modifikuj tak, aby se vzájemně neprotínaly, mohou však splývat s hranami hranic z druhé oblasti. Potom mají tyto vlastnosti − hrana splývá s jinou z druhé oblasti − hrana leží celá uvnitř druhé oblasti − hrana leží celá vně oblasti 3. Do seznamu zařaď ty hrany, které buď, leží celé v druhé oblasti, nebo splývají s nějakou hranou z druhé oblasti, se kterou mají stejnou orientaci, (totožné hrany jen jednou). 4. Z vybudovaného seznamu zřetěz hranice výsledné oblasti a výsledek pošli na výstup. (Nástin důkazu, že 4. Je uskutečnitelný…) Obalová zóna poloměru 𝑚 (buffer zone): Je polygon, jehož vnitřní body mají vzdálenost od vstupního objektu maximálně 𝑚. FI MUNI, Milan Drášil 2019 Úvod do GIS 82 Algoritmus – Obalová zóna linie: 1. Obalíme jednotlivé segmenty (úsečky) lomené čáry. 2. Obaly segmentů sjednotíme. Algoritmus – Obalová zóna oblasti: 1. Obalíme všechny hranice předešlým algoritmem. 2. Výsledek sjednotíme se vstupem (polygonem, oblastí) FI MUNI, Milan Drášil 2019 Úvod do GIS 83 Příklad: V GIS systému máme (mimo jiné...) vrstvu lesů reprezentovanou jako areály a vrstvu venkovních úseků vysokého napětí. Zajímá nás, kde lesy zasahují do bezpečnostního pásma 10 metrů kolem úseků vysokého napětí. Jednotlivé úseky obalíme buffer zónou o daném poloměru. Obalové zóny úseků nemusí být disjunktní, sjednotíme je. Výsledek pronikneme s vrstvou lesů. FI MUNI, Milan Drášil 2019 Úvod do GIS 84 6. Rastrová data v GIS Typy rastrových dat používaných pro GIS technologie jsou stejná jako v počítačové grafice: − Binární − Polotónová − Víceúrovňová FI MUNI, Milan Drášil 2019 Úvod do GIS 85 Pro další práci očíslujeme sousedy pixelu P následujícím způsobem: 3 2 1 4 P 0 5 6 7 Sousedé 0,2,4,6 se nazývají přímí (d-) sousedé pixelu p. Sousedé 1,3,5,7 se nazývají nepřímí (i-) sousedé pixelu p. Pro polotónová obraz 𝑓 označme 𝑓(𝑖, 𝑗) barvu jeho pixelu na pozici 𝑖, 𝑗. 6.1. KVANTITATIVNÍ CHARAKTERISTIKY OBRAZU Definice - Histogram obrazu: Nechť 𝑓 je polotónový obraz barev 1. . 𝑀. Jeho histogramem rozumíme konečnou posloupnost ℎ(𝑓) = (ℎ1. . ℎ 𝑀), kde, ℎ𝑖 je počet pixelů s barvou 𝑖. 2550 počet hodnota Obr. 18 - Histogram obrazu Definice - Matice sousednosti: Nechť 𝑓 je polotónový obraz barev 1. . 𝑀. Jeho maticí sousednosti rozumíme čtvercovou 𝑀𝑥𝑀 matici 𝐶𝑀(𝑓) = {𝑐𝑚𝑖,𝑗}, kde, 𝑐𝑚𝑖,𝑗 je počet (přímo) sousedících pixelů o barvě 𝑖 a 𝑗. FI MUNI, Milan Drášil 2019 Úvod do GIS 86 6.2. OPERACE NAD OBRAZY Lineární filtrace: Buď 𝑓 polotónový obraz, 𝑀 > 0. Položme 𝑔(𝑥, 𝑦) = 𝐻(𝑀, 𝑥, 𝑦) kde 𝐻 je libovolná funkce, která v konstantním čase počítá novou hodnotu pixelu 𝑔(𝑥, 𝑦) z okolí pixelu 𝑥, 𝑦 o rozměru 𝑀. Funkce 𝐻 bývá většinou vyjádřena váženým průměrem pixelů a lze ji vyjádřit: 𝐻(𝑀, 𝑥, 𝑦) = Σ𝑖=−𝑀,𝑀Σ𝑗=−𝑀,𝑀 ℎ(𝑖, 𝑗) ∗ 𝑓(𝑥 + 𝑖, 𝑦 + 𝑗) Zhlazující filtry: a) b) 1 1 1 1 2 1 1 1 1 2 4 2 1 1 1 1 2 1 FI MUNI, Milan Drášil 2019 Úvod do GIS 87 Zostřující filtry: -1 -1 -1 -1 n -1 -1 -1 -1 FI MUNI, Milan Drášil 2019 Úvod do GIS 88 Detektory hran: Základním filtr této třídy přiřadí novému pixelu největší absolutní hodnotu ze dvou výsledků: 1 2 1 1 0 -1 0 0 0 2 0 -2 -1 -2 -1 1 0 -1 Původní obraz FI MUNI, Milan Drášil 2019 Úvod do GIS 89 Zhlazující filtr Zostřující filtr Detektor hran: Kontura (hranice) oblasti v binárním obrazu: Nechť O je libovolná oblast (množina složená z jedniček) v binárním obrazu, Konturou (hranicí) oblasti O rozumíme všechny pixely patřící této oblasti, které mají nulového d-souseda. 1 1 1 11 1 1 1 000 0 0 0 0 0 rastrověvektorově FI MUNI, Milan Drášil 2019 Úvod do GIS 90 6.3. TRANSFORMACE RASTROVÝCH DAT Velmi častá úloha, zejména pro umisťování obrazu do dané kartografické projekce. 1. Určíme bodové objekty ve zdrojovém obrazu, jejichž (kartografické) souřadnice jsou známé. Například přesně odečtené z mapy, geodeticky zaměřené apod. 2. Určíme transformační funkce z nového do starého obrazu komerční produkty většinou poskytují polynomiální transformaci založené na metodě nejmenších čtverců, je však možné použít i přesnou kartografickou transformaci. 𝑖 = 𝐹(𝑥, 𝑦) 𝑗 = 𝐺(𝑥, 𝑦) kde (𝑖, 𝑗) značí souřadný systém originálního obrazu, (𝑥, 𝑦) souřadný systém obrazu nového. 3. Procházíme cílový obraz, pomocí transformačních funkcí 𝐹 a 𝐺 se „díváme“ do originálu a počítáme hodnotu pixelu. Podle toho, z jakého okolí zdrojového pixelu určujeme výslednou hodnotu pixelu nového, se hovoří o metodách: − nejbližší soused − bilineární transformace − konvoluce okolí 𝑀𝑥𝑀 První případ prostě přenese hodnotu pixelu do nového obrazu, v dalších případech se výsledná hodnota se počítá z jistého okolí pixelu v originálním obrazu. Velmi často se setkáváme s následující situací (například v klientech WMS služby). - Známe omezující obdélník cílové kartografické projekce a rozměr obrazu v pixelech - Známe sadu projekcí, které poskytuje (WMS) server - Potřebujeme dostat mapovou kompozici do “okna” na klientskou stanici. [i,j] [x,y] FI MUNI, Milan Drášil 2019 Úvod do GIS 91 V tomto případě postupujeme takto: 1. Vybereme vhodnou serverovou projekci (např. 𝑊𝐺𝑆84, tu poskytuje každý WMS server). 2. Transformujeme dotazový obdélník do serverové projekce, upravíme pixelový rozměr obrazu. K tomuto účelu použijeme ”přesný” převod souřadnic: 𝑟𝑜𝑣𝑖𝑛𝑎1 → 𝑒𝑙𝑖𝑝𝑠𝑜𝑖𝑑1→ 𝑐𝑒𝑛𝑡𝑟1 → 𝑐𝑒𝑛𝑡𝑟2 → 𝑒𝑙𝑖𝑝𝑠𝑜𝑖𝑑2 → 𝑟𝑜𝑣𝑖𝑛𝑎2 (viz transformace souřadných systémů). 3. Provedeme dotaz, získáme obraz v souřadnicích (𝑖, 𝑗). 4. Vybereme sadu zdrojových pixelových souřadnic (většinou pravidelná mřížka 𝑁𝑥𝑁, nebo rohy obrazu). Tyto převedeme do zdrojových souřadnic, poté “přesnou” kartografickou transformací do souřadnic cílové projekce a cílových souřadnic pixelových (𝑥, 𝑦). Takto získáme sadu dvojic “odpovídajících si” pixelů, kterou použijeme pro získání parametrů polynomiální transformace obrazu. Tento postup “šetří” výpočtovou náročnost tím, že nahradí náročné přepočty kartografických souřadnic polynomiální transformací. Skelet binárního obrazu: Při vektorizaci liniových prvků je prvním krokem zúžení rastrových linií tak, aby byly maximálně dva pixely široké. ..... ..... .*.*.. .*.*. ..*.*.. ..*.. ..*..*.. ..*.. ..*...*.. ..*.. ..*....*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*.. ..*....*.. ..*.. ..*...*.. ..*.. ..*..*.. ..*.. ..*.*.. .*.*. ..*.*. ..... ..... Definice – Skelet: FI MUNI, Milan Drášil 2019 Úvod do GIS 92 Nechť 𝑅 je množina pixelů, 𝐵 její hranice (kontura), 𝑃 bod v 𝑅. Nejbližší soused bodu 𝑃 na hranici 𝐵 je bod 𝑀 z 𝐵 takový, že pro každý bod 𝑀′ z 𝐵 je vzdálenost 𝑃𝑀′ větší nebo rovna vzdálenosti 𝑃𝑀. Má-li bod 𝑃 více než jednoho nejbližšího souseda, nazývá se bodem skeletu množiny 𝑅. Sjednocení všech bodů skeletu tvoří skelet množiny 𝑅. Algoritmus - Určení skeletu: 1. Urči hranici (konturu) 𝐵(𝑅) množiny 𝑅. 2. Urči množinu násobných pixelů 𝑀(𝑅) v hranici 𝐵(𝑅) 3. Je-li 𝐵(𝑅) = 𝑀(𝑅), skonči. 4. Polož 𝑅 = 𝑅 − (𝐵(𝑅) − 𝑀(𝑅)) a pokračuj krokem 1. (a) (b) (c) A A A A A A A A C 0 P 0 A P 0 0 P 2+ B B B A 0 2 B B C Pixel označený 2 je hraniční pixel, pixel označený 2+ je hraniční nebo násobný pixel. (a), (b): Alespoň jeden pixel ze skupiny pixelů A, B je nenulový (c): Alespoň jeden pixel ze skupiny C musí být nenulový. Pokud jsou oba nenulové, pak může být hodnota pixelů ve skupinách A i B libovolná. Pokud je jeden pixel skupiny C nulový, musí být alespoň jeden pixel skupiny A i skupiny B nenulový. FI MUNI, Milan Drášil 2019 Úvod do GIS 93 7. Topologické úlohy v GIS Byly efektivně řešeny před tím, než byly vůbec GIS technologie vymezeny, přesto je můžeme považovat za součást analytického jádra topologicky orientovaných GIS. Základní datová struktura síťového grafu: Základní datová struktura areálového grafu: Nejčastější úlohy: Trasování grafu: Vyber všechny uzly/hrany, které jsou “napájeny” z daného uzlu, hodně používaná úloha pro dispečery sítí, modelování situací „co se stane, když?“. Nejkratší cesta z uzlu do uzlu: Používá se klasický Dijkstrův algoritmus. Lze výhodně využít vlastnosti, že uzly grafu jsou prostorově lokalizovány. Uzel Vedení Přípojka HranaHrana Uzávěr Odběrné místo začíná končí Areál Hr. obce Hr. katas HranaHrana Obec Katastr 1.areál 2.areál FI MUNI, Milan Drášil 2019 Úvod do GIS 94 Algoritmus „Minimální cesta (Best search): Nechť (𝑈, 𝐻) je graf s nezáporně ohodnocenými hranami, váhu libovolné hrany ℎ𝜖𝐻 označme 𝑤(ℎ). Nechť dále každý uzel 𝑢𝑖 𝜖𝑈 má 2D souřadnice (𝑥𝑖, 𝑦𝑖). Pro libovolné uzly 𝑢, 𝑣 označme 𝑑(𝑢, 𝑣) jejich vzdálenost v 𝐸2. Pro libovolnou hranu ℎ𝑖 = (𝑢𝑖, 𝑣𝑖) nechť dále platí trojúhelníková nerovnost metrického prostoru 𝐸2. 𝑑(𝑢𝑖, 𝑣𝑖) ≤ 𝑤(ℎ) Potom pro libovolné uzly 𝑢0, 𝑣𝜖𝑈, následující postup vede k nalezení minimální cesty z 𝑢0 do 𝑣. 1. Inicializace: Pro každý uzel 𝑢𝑖 𝜖𝑈, 𝑢𝑖 ≠ 𝑢0 položme 𝑑_𝑝𝑎𝑡ℎ𝑖 = ∞, 𝑑_𝑝𝑎𝑡ℎ0 = 0, dále položme každý uzel 𝑢𝑖 𝜖𝑈 𝑐_𝑝𝑎𝑡ℎ𝑖 = 𝑛𝑢𝑙𝑙. 2. Výběr pivota: Položme 𝑐_𝑝𝑎𝑡ℎ𝑖 = 𝑑_𝑝𝑎𝑡ℎ𝑖 pro takový 𝑢𝑖, pro který je 𝑐_𝑝𝑎𝑡ℎ𝑖 = 𝑛𝑢𝑙𝑙, 𝑑_𝑝𝑎𝑡ℎ𝑖 < ∞ a pro který je 𝑑_𝑝𝑎𝑡ℎ𝑖 + 𝑑(𝑢𝑖, 𝑣) minimální. Když neexistuje – končíme, cesta neexistuje. Je-li 𝑢𝑖 = 𝑣, potom konec, 𝑐_𝑝𝑎𝑡ℎ𝑖 je délka minimální cesty a provedeme „zpětný“ chod. 3. Expanze: Pro všechny uzly 𝑢 𝑘 𝜖𝑈 takové, že existuje hrana ℎ 𝑘 𝜖𝐻, ℎ 𝑘 = (𝑢𝑖, 𝑢 𝑘) (𝑢𝑖 je pivot z předchozího kroku) položme 𝑑_𝑝𝑎𝑡ℎ 𝑘 = 𝑚𝑖𝑛{𝑑_𝑝𝑎𝑡ℎ 𝑘, 𝑐_𝑝𝑎𝑡ℎ𝑖 + 𝑤(ℎ 𝑘)} a pokračujeme 2. Problém obchodního cestujícího: Hamiltonovská kružnice v grafu, extrémně obtížná úloha, dosud nebyla uspokojivě vyřešena (jedná se o NP úplný problém). Topologicko-geometrické úlohy: - Vytváření topologických vazeb na základě geometrických vlastností objektů; jedná se o vytvoření příslušného typu grafu (uzel-hrana, hrana-hrana, areálový graf). Hojně se využívá přístupových metod pro geometrické objekty. - Generování oblastí z hran areálového grafu. - Identifikace hran areálového grafu. - Generování vyšších územních celků areálového grafu. - Kontrola konzistence geometrických a topologických vlastností dat (kontrola shody umístění uzlu a konce hrany s ním incidentní, kontrola křížení hran, kontrola stupňů uzlových bodů a další kontroly). FI MUNI, Milan Drášil 2019 Úvod do GIS 95 8. 3D geometrie v GIS 8.1. 3D GEOMETRICKÉ PRIMITIVY Jsou zcela analogické 2D s tím rozdílem, že vycházejí z 3 dimensionálních bodů [𝑥, 𝑦, 𝑧]. 8.2. ODHAD NORMÁLY MNOŽINY 3D BODŮ, KOVARIANČNÍ MATICE BODŮ 1) 𝑎𝑥 + 𝑏𝑦 + 𝑐𝑧 + 𝑑 = 0 je rovnice roviny v prostoru. 2) [𝑥𝑖 𝑝 , 𝑦𝑖 𝑝 , 𝑧𝑖 𝑝 ], 𝑖 = 1, . . , 𝑛 jsou body vstupní body 3) Úloha: Pro body [𝑥𝑖 𝑝 , 𝑦𝑖 𝑝 , 𝑧𝑖 𝑝 ], 𝑖 = 1, . . , 𝑛 hledáme 𝑎, 𝑏, 𝑐, 𝑑 , co nejlepší aproximaci roviny z 1), resp. její normálový vektor [𝑎, 𝑏, 𝑐]. 4) Nechť bod [𝑥̅, 𝑦̅, 𝑧̅] = [ ∑ 𝑥 𝑖 𝑝 𝑛 , ∑ 𝑦𝑖 𝑝 𝑛 , ∑ 𝑧𝑖 𝑝 𝑛 ] je centroid bodů z 2). 5) Položíme [𝑥𝑖, 𝑦𝑖, 𝑧𝑖] = [𝑥𝑖 𝑝 − 𝑥̅, 𝑦𝑖 𝑝 − 𝑦̅, 𝑧𝑖 𝑝 − 𝑧̅], 𝑖 = 1, . . , 𝑛. Triviálně, posunem bodů „do centroidu“ se normálový vektor nezmění. Dále, po této úpravě máme: ∑ 𝑥𝑖 = ∑ 𝑦𝑖 = ∑ 𝑧𝑖 = 0 6) Použijeme metodu nejmenších čtverců: 𝑟 = ∑(𝑎𝑥𝑖 + 𝑏𝑦𝑖 + 𝑐𝑧𝑖 + 𝑑)2 = 𝑚𝑖𝑛 7) Derivujeme 𝑟 z 6) podle 𝑎, 𝑏, 𝑐 a 𝑑 a dostaneme soustavu rovnic (4x4): 𝜕𝑟 𝜕𝑎 = 𝜕𝑟 𝜕𝑏 = 𝜕𝑟 𝜕𝑐 = 𝜕𝑟 𝜕𝑑 = 0 Maticově: [ ∑ 𝑥𝑖 𝑥𝑖 ∑ 𝑦𝑖 𝑥𝑖 ∑ 𝑧𝑖 𝑥𝑖 ∑ 𝑥𝑖 ∑ 𝑥𝑖 𝑦𝑖 ∑ 𝑦𝑖 𝑦𝑖 ∑ 𝑧𝑖 𝑦𝑖 ∑ 𝑦𝑖 ∑ 𝑥𝑖 𝑧𝑖 ∑ 𝑦𝑖 𝑧𝑖 ∑ 𝑧𝑖 𝑧𝑖 ∑ 𝑧𝑖 ∑ 𝑥𝑖 ∑ 𝑦𝑖 ∑ 𝑧𝑖 𝑛 ] × [ 𝑎 𝑏 𝑐 𝑑 ] = [ 0 0 0 0 ] FI MUNI, Milan Drášil 2019 Úvod do GIS 96 8) Z 5) a posledního řádku matice soustavy máme 𝑑 = 0. Tedy hledáme netriviální řešení (tj. [𝑎, 𝑏, 𝑐] ≠ [0,0,0]): [ ∑ 𝑥𝑖 𝑥𝑖 ∑ 𝑦𝑖 𝑥𝑖 ∑ 𝑧𝑖 𝑥𝑖 ∑ 𝑥𝑖 𝑦𝑖 ∑ 𝑦𝑖 𝑦𝑖 ∑ 𝑧𝑖 𝑦𝑖 ∑ 𝑥𝑖 𝑧𝑖 ∑ 𝑦𝑖 𝑧𝑖 ∑ 𝑧𝑖 𝑧𝑖] × [ 𝑎 𝑏 𝑐 ] = [ 0 0 0 ] Povšimněme si, že tato soustava má vždy triviální řešení, netriviální řešení dostaneme pouze tehdy, je-li matice této soustavy singulární. O triviální řešení ale zájem nemáme, zajímá nás řešení netriviální. 9) Je-li [𝑎, 𝑏, 𝑐, 0] řešením 1), potom je triviálně i 𝑘[𝑎, 𝑏, 𝑐, 0], 𝑘𝜖𝑹 řešením 1) 10) Dále, v netriviálním řešení musí platit, že buď 𝑎 ≠ 0, nebo 𝑏 ≠ 0, nebo 𝑐 ≠ 0. Předpokládejme tedy, že 𝑐 ≠ 0 a vzhledem k 9) můžeme předpokládat 𝑐 = 1. Z 8) tedy dostáváme: [ ∑ 𝑥𝑖 𝑥𝑖 ∑ 𝑦𝑖 𝑥𝑖 ∑ 𝑧𝑖 𝑥𝑖 ∑ 𝑥𝑖 𝑦𝑖 ∑ 𝑦𝑖 𝑦𝑖 ∑ 𝑧𝑖 𝑦𝑖 ∑ 𝑥𝑖 𝑧𝑖 ∑ 𝑦𝑖 𝑧𝑖 ∑ 𝑧𝑖 𝑧𝑖 ] × [ 𝑎 𝑏 1 ] = [ 0 0 0 ]¨ a tedy: [ ∑ 𝑥𝑖 𝑥𝑖 ∑ 𝑦𝑖 𝑥𝑖 ∑ 𝑥𝑖 𝑦𝑖 ∑ 𝑦𝑖 𝑦𝑖 ] × [ 𝑎 𝑏 ] = [ − ∑ 𝑧𝑖 𝑥𝑖 − ∑ 𝑧𝑖 𝑦𝑖 ] V případě nulového determinantu této matice zaměníme předpoklad z 10) 𝑐 = 1 za 𝑏 = 1, resp. 𝑎 = 1. Máme řešení úlohy z 3). 8.3. 3D POLYGON Ve 2D případě požadujeme po hranicích polygonů, aby se vzájemně neprotínaly. Ve 3D variantě požadujeme, aby tato vlastnost byla splněna pro ortogonální projekci do ‚nejlepší‘ proložené roviny z předešlého odstavce. Tato vlastnost nám především umožňuje jejich kresbu pomocí triangulace. Postup je následující: 1) Posuneme polygon do jeho centroidu. 2) Provedeme odhad „nejlepší“ normály. 3) Rotujeme polygon v rovině XY a poté v rovině XZ tak, aby jeho normála z 2) byla rovnoběžná s osou Z. 4) Provedeme triangulaci, zohledňujeme jen souřadnice XY. FI MUNI, Milan Drášil 2019 Úvod do GIS 97 5) Pro výsledné trojúhelníky provedeme zpětnou transformaci z kroků 3) a 1). 8.4. MRAČNA BODŮ Mračna bodů jsou moderním zdrojem pro získávání geometrické (geograficky vztažené) informace. Vznikají snímáním (leteckým/pozemním) zemského povrchu laserovými scannery, výsledkem jsou 3D body s doplňující informací (intensita odrazu, barva). Jsou definovány standardy formátů, ve kterých jsou mračna poskytována (LAS). Zobrazení mračna bodů v centrální projekci 8.5. OSVĚTLENÍ MRAČNA PODLE NORMÁL V některých případech je například vhodné pro vizualizaci mračna použít jeho přebarvení podle prostorových vlastností. Základní metodou pro zvýraznění je ‚osvětlení‘ z určitého zdrojového bodu (vektoru), tj. jednotlivé body obarvit podle úhlu normály lokální roviny vůči vektoru osvětlení. Lokální rovinou rozumíme opět nejlepší proložení rovinou jistým (relativně malým) okolím jednotlivých bodů. Zdůrazněme, že tyto metody vyžadují velmi výkonnou prostorovou indexaci mračna bodů. Není výjimkou, že zpracováváme řádově 107 bodů a bez prostorové indexace, a tedy s kvadratickou časovou složitostí se dostáváme mimo rozumný výkon (1014 je hodně, i když počítače jsou výkonné). Pro každý bod mračna totiž potřebujeme hodně rychle vybrat všechny body z jeho okolí. FI MUNI, Milan Drášil 2019 Úvod do GIS 98 Zobrazení mračna bodů obarveného podle úhlu normálových vektorů 8.6. VLASTNÍ HODNOTY A VEKTORY KOVARIANČNÍ MATICE Připomeňme si z lineární algebry, co jsou to vlastní vektory a vlastní hodnoty matice. Buď 𝐴 čtvercová matice, vlastním vektorem 𝑢 nazveme takový vektor, pro který platí: 𝐴𝑢 = 𝜆𝑢 kde 𝜆 je číslo, které nazýváme vlastní hodnotou matice 𝐴. Vlastní vektor je tedy takový vektor, který aplikací transformace reprezentované maticí 𝐴 nemění směr, mění se pouze jeho velikost faktorem 𝜆. Některé užitečné vlastnosti vlastních čísel a vektorů matic: - Nula je vlastním číslem matice právě tehdy, když je matice singulární. - Je-li matice symetrická a reálná (tj. obsahuje pouze reálná čísla), pak všechna její vlastní čísla jsou reálná. (*) - Vlastní vektory reálné symetrické matice jsou po dvou ortogonální (tj. jejich skalární součin je nulový, jsou na sebe kolmé). (*) (*) jsou základních vlastností Hermitovských (samosdružených) operátorů v Hilbertově vektorovém prostoru. FI MUNI, Milan Drášil 2019 Úvod do GIS 99 Matici nazveme positivně definitní, resp. positivně semidefinitní, jsou-li všechna její vlastní čísla kladná, resp. nezáporná. Položme, viz 8. 2 8): 𝐶 = [ ∑ 𝑥𝑖 𝑥𝑖 ∑ 𝑦𝑖 𝑥𝑖 ∑ 𝑧𝑖 𝑥𝑖 ∑ 𝑥𝑖 𝑦𝑖 ∑ 𝑦𝑖 𝑦𝑖 ∑ 𝑧𝑖 𝑦𝑖 ∑ 𝑥𝑖 𝑧𝑖 ∑ 𝑦𝑖 𝑧𝑖 ∑ 𝑧𝑖 𝑧𝑖] × [ 1 𝑛 − 1 0 0 0 1 𝑛 − 1 0 0 0 1 𝑛 − 1] 𝐶 se nazývá kovarianční maticí. Platí např. následující tvrzení. 0 je vlastní hodnota matice 𝐶 právě když [𝑥𝑖, 𝑦𝑖, 𝑧𝑖], 𝑖 = 1, . . , 𝑛 leží v rovině (tj. existuje jejich dvouprvková báze). Dále, všechny vlastní hodnoty matice 𝐶 jsou reálné, neboť 𝐶 je symetrická a reálná. Kovarianční matice 𝐶 je navíc positivně semidefinitní, takže každá její vlastní hodnota 𝜆 ≥ 0. Poznámka: Vlastní hodnoty matice lze spočítat např. Jacobiho iterační metodou, viz např.: https://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm Z pohledu zpracování mračna bodů má kovarianční matice 3D bodů a zejména její vlastní hodnoty a vektory jednu velmi cennou vlastnost. Její nenulové vlastní vektory jsou totiž „rozumnou“ ortogonální bází vstupních bodů, jinými slovy tyto body lze vyjádřit pomoci lineární kombinace této báze. Jsou-li navíc odpovídající vlastní hodnoty příslušné těmto vektorům relativně malé vůči ostatním, můžeme je zanedbat. Vlastní vektory kovarianční matice okolí bodu, velikosti podle poměru odpovídajících vlastních hodnot FI MUNI, Milan Drášil 2019 Úvod do GIS 100 8.7. VARIANCE POVRCHU Nechť 𝜆0, 𝜆1, 𝜆2 jsou vlastní hodnoty kovarianční matice 𝐶 z bodů z okolí bodu 𝑝, 𝜆0 ≤ 𝜆1 ≤ 𝜆2. Položme: 𝜎(𝑝) = 𝜆0 𝜆0 + 𝜆1 + 𝜆2 Operátor 𝜎(𝑝) nazveme povrchovou variancí bodu 𝑝. 𝜎(𝑝) = 0 právě když body v okolí 𝑝 leží v (ideální) rovině. Čím je hodnota 𝜎(𝑝) větší, tím je větší „nerovnost“ jeho okolí. Triviálně: 𝜎(𝑝)𝜖 [0, 1 3 ] Vizualizace variance povrchu – čím světlejší, tím větší míra nerovnosti 8.8. AUTOMATICKÉ TRASOVÁNÍ LINIÍ V MRAČNU BODŮ Tyto metody se snaží z mračna bodů, ve kterém jsou klasifikovány jisté body jako kandidáti potenciální trasy linie, převést ‚mračnovou‘ linii do 3D lomené linie. Kandidáty lze získat metodami variance povrchu (viz předešlý odstavec), výškovým omezením v případě nadzemního vedení aj. Princip spočívá v tom, že metody FI MUNI, Milan Drášil 2019 Úvod do GIS 101 klasifikace kandidátů musí od sebe jednotlivé linie oddělit (nejjednodušším případem jsou nadzemní, např. trolejová vedení – ty jsou od sebe oddělena přirozeně). Opět využijeme vlastností kovarianční matice okolí zkoumaného bodu. Na rozdíl od metod variance povrchu, kde jsme hledali dva dominantní vektory tvořící bázi plochy, nyní hledáme jeden dostatečně dominantní vektor, který určuje lokální směr linie. Zkoumaný bod posunujeme podle takto zjištěných vektorů. Zdrojové mračno pro trasování linií Trasované linie FI MUNI, Milan Drášil 2019 Úvod do GIS 102 Reference ke 3D: [1] Pauly, M., Gross, M., Kobbelt, L. Efficient Simplification of Point-Sampled Surfaces. Visualization, VIS. IEEE, pp. 163-170, 2002. [2] Pauly, M., Keiser, R., Gross, M. Multi-scale feature extraction on pointsampled surfaces. Computer Graphics Forum, 22(3), pp. 281-289, 2003. [3] Dena Bazazian, Josep R. Casas, Javier Ruiz-Hidalgo. Fast and Robust Edge Extraction in Unorganized Point Clouds, https://imatge.upc.edu/web/sites/default/files/pub/cBazazian15.pdf FI MUNI, Milan Drášil 2019 Úvod do GIS 103 9. Přílohy 9.1. METODA ŘEŠENÍ NELINEÁRNÍCH ROVNIC, NEWTON–RAPHSON Snažíme se vyřešit rovnici 𝑓(𝑥) = 0 kde 𝑓(𝑥) má první derivaci a my ji navíc dokážeme spočítat. Uvažme Taylorův rozvoj funkce: 𝑓(𝑥) = 𝑓(𝑎) + 𝑓′(𝑎) 1! (𝑥 − 𝑎) + 𝑓″(𝑎) 2! (𝑥 − 𝑎)2 … V bodě 𝑎 můžeme nahradit funkci 𝑓 její tečnou 𝑓(𝑥) ≈ 𝑓(𝑎) + 𝑓′(𝑎) 1! (𝑥 − 𝑎) Rovnici: 𝑓(𝑎) + 𝑓′(𝑎) 1! (𝑥 − 𝑎) = 0 𝑥 = 𝑎 − 𝑓(𝑎) 𝑓′(𝑎) už řešit umíme. Postup: Iterace lineární aproximace rovnice. Odhadneme první hodnotu 𝑎 = 𝑥0 a řešíme lineární rovnici, výsledek bude 𝑥1. Položíme 𝑎 = 𝑥1 a postup opakujeme tak dlouho dokud rozdíl |𝑥𝑖 − 𝑥𝑖−1| nebude dostatečně malý. Proces NEMUSÍ konvergovat! U některých funkcí hodně záleží na počátečním odhadu! Iterace Newton–Raphson metody funkce jedné proměnné FI MUNI, Milan Drášil 2019 Úvod do GIS 104 Příklad 1: 𝑓: 𝑦 = 3𝑥 + 8 𝑓′ : 𝑦 = 3 Položme 𝑥0 = 𝑎 = 10 a dosazením do 𝑥1 = 𝑎 − 𝑓(𝑎) 𝑓′(𝑎) = 10 − 38 3 = − 3 8 V případě lineární funkce je tečna totožná z funkcí, stačí jedna iterace. Příklad 2: 𝑓: 𝑦 = 𝑥2 4 + √ 𝑥 − 𝑥 − 2 𝑓′ : 𝑦 = 1 2 𝑥 + 1 2√ 𝑥 − 1 S počáteční hodnotou 𝑥0 = 12345 potřebujeme cca 18 iterací: a a-f(a)/f´(a) f(a) f´(a) 12345 6173.486983 38087520.36 6171.5045 6173.486983 3087.725365 9521888.466 3085.749855 3087.725365 1544.837609 2380477.824 1542.87168 1544.837609 773.3844645 595123.2765 771.4315258 773.3844645 387.6458883 148783.3078 385.7102115 387.6458883 194.761798 37197.37652 192.8483394 194.761798 98.3030846 9300.2334 96.41672663 98.3030846 50.05846472 2325.485817 48.20197201 50.05846472 25.93059835 581.4792084 24.09990173 25.93059835 13.88925597 145.2605939 12.06348838 13.88925597 7.956261183 36.06543106 6.078790288 7.956261183 5.202262665 8.689946344 3.155392528 5.202262665 4.18901228 1.844468463 1.820348149 4.18901228 4.006273118 0.244651359 1.33880092 4.006273118 4.000007362 0.00785062 1.252940754 4.000007362 4 9.20245E-06 1.250003451 4 4 1.27027E-11 1.25 4 4 0 1.25 FI MUNI, Milan Drášil 2019 Úvod do GIS 105 V případě funkcí více proměnných a systému nelineárních rovnic 𝑓1(𝑥, 𝑦, 𝑧, . . ) = 0 𝑓2(𝑥, 𝑦, 𝑧, . . ) = 0 . je postup stejný. Funkce linearizujeme 𝑓𝑖(𝑥, 𝑦, 𝑧, . . ) ≈ 𝑓𝑖(𝑥0, 𝑦0, 𝑧0, . . ) + 𝜕𝑓𝑖 𝜕𝑥 (𝑥0, 𝑦0, 𝑧0, . . )(𝑥 − 𝑥0) + 𝜕𝑓𝑖 𝜕𝑦 (𝑥0, 𝑦0, 𝑧0, . . )(𝑦 − 𝑦0) + 𝜕𝑓𝑖 𝜕𝑧 (𝑥0, 𝑦0, 𝑧0, . . )(𝑧 − 𝑧0) … V každém kroku iterace potom řešíme systém lineárních rovnic. Označíme-li 𝑿𝒊 = [𝑥𝑖, 𝑦𝑖, 𝑧𝑖, … ] a Δ𝑿 = [𝑥 − 𝑥𝑖, 𝑦 − 𝑦𝑖, 𝑧 − 𝑧𝑖, … ], potom v každém kroku řešíme soustavu: [ 𝜕𝑓1 𝜕𝑥 (𝑿𝒊) 𝜕𝑓1 𝜕𝑦 (𝑿𝒊) 𝜕𝑓1 𝜕𝑧 (𝑿𝒊) … 𝜕𝑓2 𝜕𝑥 (𝑿𝒊) 𝜕𝑓2 𝜕𝑦 (𝑿𝒊) 𝜕𝑓2 𝜕𝑧 (𝑿𝒊) … 𝜕𝑓3 𝜕𝑥 (𝑿𝒊) 𝜕𝑓3 𝜕𝑦 (𝑿𝒊) 𝜕𝑓3 𝜕𝑧 (𝑿𝒊) … ⋮ ⋮ ⋮ ⋱] × Δ𝑿 = [ −𝑓1(𝑿𝑖) −𝑓2(𝑿𝑖) −𝑓3(𝑿𝑖) ⋮ ] Je-li řešením této soustavy Δ𝒙𝑖, položíme 𝑿𝑖+1 = 𝑿𝑖 + Δ𝒙𝑖 a pokračujeme, dokud |Δ𝒙𝑖| ≈ 0. FI MUNI, Milan Drášil 2019 Úvod do GIS 106 9.2. METODA NEJMENŠÍCH ČTVERCŮ – LINEÁRNÍ REGRESE Hledáme „co nejlepší“ řešení soustavy lineárních rovnic (𝑀 > 𝑁): [ 𝑎11 𝑎12 ⋯ 𝑎1𝑁 𝑎21 𝑎22 ⋯ 𝑎2𝑁 ⋮ ⋮ ⋮ ⋮ 𝑎 𝑀1 𝑎 𝑀2 ⋯ 𝑎 𝑀𝑁 ] [ 𝑥1 𝑥2 ⋮ 𝑥 𝑁 ] ≈ [ 𝑏1 𝑏2 ⋮ 𝑏 𝑀 ] 𝑨𝑥 ≈ 𝑏 (1) To provedeme tak, že budeme hledat takové řešení, pro které je součet čtverců odchylek jednotlivých rovnic co nejmenší (Carl Friedrich Gauss 1795), tedy: Δ𝑥 = (𝑨𝑥 − 𝑏) 𝑇(𝑨𝑥 − 𝑏) = 𝑚𝑖𝑛 (2) kde operátor 𝑨 𝑇 označuje transpozici matice 𝑨, připomeňme jeho základní vlastnosti: (𝑨 𝑇) 𝑇 = 𝑨, (𝑨 + 𝑩) 𝑇 = 𝑨 𝑇 + 𝑩 𝑇 , (𝑨𝑩) 𝑇 = 𝑩 𝑇 𝑨 𝑇 roznásobením (2) dostaneme: Δ𝑥 = (𝑨𝑥) 𝑇 𝑨𝑥 − (𝑨𝑥) 𝑇 𝑏 − 𝑏 𝑇 𝑨𝑥 + 𝑏 𝑇 𝑏 = 𝑥 𝑇 𝑨 𝑇 𝑨𝑥 − 𝑥 𝑇 𝑨 𝑇 𝑏 − 𝑏 𝑇 𝑨𝑥 + 𝑏 𝑇 𝑏 matice 𝑏 𝑇 𝑨𝑥 má rozměr 1 × 1, tedy 𝑏 𝑇 𝑨𝑥 = (𝑏 𝑇 𝑨𝑥) 𝑇 = (𝑨𝑥) 𝑇 𝑏 = 𝑥 𝑇 𝑨 𝑇 𝑏 a: Δ𝑥 = 𝑥 𝑇 𝑨 𝑇 𝑨𝑥 − 2𝑥 𝑇 𝑨 𝑇 𝑏 + 𝑏 𝑇 𝑏 Hledáme minimum Δ𝑥, tedy derivujeme Δ𝑥 podle 𝑥 = [𝑥1, … , 𝑥 𝑁] a výraz položíme roven 0. Položme tedy: 𝜕Δ𝑥 𝜕𝑥 = [ 𝜕Δ𝑥 𝜕𝑥1 ⋮ 𝜕Δ𝑥 𝑥 𝑁 ] = 2𝑨 𝑇 𝑨𝑥 − 2𝑨 𝑇 𝑏 = [ 0 ⋮ 0 ] (derivace 𝜕Δ𝑥 𝜕𝑥 není na první pohled úplně zřejmá, dá se však po jistém úsilí odvodit). Dostáváme soustavu tzv. normálních rovnic, jejímž řešením je hledaná hodnota 𝑥. 𝑨 𝑇 𝑨𝑥 = 𝑨 𝑇 𝑏 Vraťme se zpátky k zadání (1): 𝑨𝑥 ≈ 𝑏 „Nejlepší“ řešení soustavy rovnic 𝑨𝑥 ≈ 𝑏 dostaneme, když její každou stranu zleva násobíme transponovanou maticí soustavy 𝑨 𝑇 . FI MUNI, Milan Drášil 2019 Úvod do GIS 107 Příklad: Body [1,0], [2,3], [3,2], [4,5] prolož přímku 𝑦 = 𝑎𝑥 + 𝑏. Maticový zápis této úlohy je: [ 1 1 2 1 3 1 4 1 ] [ 𝑎 𝑏 ] = [ 0 3 2 5 ] 𝑨 𝑇 𝑨 = [ 30 10 10 4 ] , 𝑨 𝑇 𝑏 = [ 32 10 ] Řešením: [ 30 10 10 4 ] [ 𝑎 𝑏 ] = [ 32 10 ] dostaneme výsledek 𝑎 = 7 5 ; 𝑏 = −1 a tedy 𝑦 = 1.4𝑥 − 1 FI MUNI, Milan Drášil 2019 Úvod do GIS 108 Příklad: Máme dvě „odpovídající si“ sady bodů v rovině [𝑥𝑖, 𝑦𝑖] a [𝑢𝑖, 𝑣𝑖]. [𝑢𝑖, 𝑣𝑖] → [𝑥𝑖, 𝑦𝑖] [1,1] → [0,5] [3,1] → [4,11] [1,3] → [0,11] [5,4] → [20,50] [3,5] → [12,39] Vypočtěte koeficienty 𝑎1, 𝑏1, 𝑐1, 𝑑1, 𝑎2, 𝑏2, 𝑐2, 𝑑2 pro bilineární transformaci: 𝑥 = 𝑎1 𝑢 + 𝑏1 𝑣 + 𝑐1 𝑢𝑣 + 𝑑1 𝑦 = 𝑎2 𝑢 + 𝑏2 𝑣 + 𝑐2 𝑢𝑣 + 𝑑2 Maticový zápis úlohy je: [ 1 1 1 1 3 1 3 1 1 3 3 1 5 4 20 1 3 5 15 1] [ 𝑎1 𝑏1 𝑐1 𝑑1 ] = [ 0 4 0 20 12] a [ 1 1 1 1 3 1 3 1 1 3 3 1 5 4 20 1 3 5 15 1] [ 𝑎2 𝑏2 𝑐2 𝑑2 ] = [ 15 11 11 50 39] Aplikací metody nejmenších čtverců 𝑨 𝑇 𝑨𝑥 = 𝑨 𝑇 𝑏 dostáváme: [ 45 42 158 13 42 52 168 14 158 168 644 42 13 14 42 5 ] [ 𝑎1 𝑏1 𝑐1 𝑑1 ] = [ 148 144 592 36 ] a [ 45 42 158 13 42 52 168 14 158 168 644 42 13 14 42 5 ] [ 𝑎2 𝑏2 𝑐2 𝑑2 ] = [ 416 444 1656 116 ] Řešením těchto dvou systémů lineárních rovnic dostaneme: 𝑎1 = 1, 𝑏1 = −1, 𝑐1 = 1, 𝑑1 = −1, 𝑎1 = 1, 𝑏2 = 1, 𝑐2 = 2, 𝑑1 = 1 a 𝑥 = 𝑢 − 𝑣 + 𝑢𝑣 − 1 𝑦 = 𝑢 + 𝑣 + 2𝑢𝑣 + 1