Přesnost výpočtů v Geoinformatice

prováděných procesorem Intel

Doc. Dr. Vladimír Homola, Ph.D.

Abstrakt

Článek se zabývá otázkou přesnosti uložení dat a výpočtů prováděných na počítačích třídy PC s procesory s množinou instrukcí procesoru Intel. Uvádí formáty dat zpracovávaných matematickým koprocesorem a podporovaných programovacími jazyky. To vše směřuje k diskusi o uložení a zpracování dat v geoinformatice a dovoluje odpovědět na otázku, zda je zpracování geoinformatických dat ještě únosně přesné. Nejpodstatnější části článku však mají obecnou platnost i pro jiné obory. Pro demonstraci přesnosti výpočtů je použit tabulkový procesor Excel.

Abstract

The aim of this article is the discussion of the precision of data storage and calculation results performed by PC processors with Intel – based instruction set. Introduces the data formats processed by mathematical co-processor and supported by programming languages. All of this is pointed to the discussion of data storing and processing in Geoinformatics. The article is then able to answer the question: Is the result of geoinfmomatics data processing still sufficiently precise? The fundamentals of the article may be similarly used in other branches. The enumeration accuracy is presented by spreadsheet Excel.

Popis problému

Problém přesnosti výpočtů je dnes poněkud zastřen množstvím specializovaného programového vybavení, které své výsledky předestírá s neochvějnou suverenitou jako (díky nezbytné reklamě dokonce jako jediné) stoprocentně správné. Uživatel takových programů si evidentně otázkou přesnosti hlavu neláme: perfektní programy nelhou.

Ovšem zkušený uživatel - geoinformatik je na tom trochu jinak. Především si vzpomene, že jest cosi shnilého už v programu Excel:

Když pak stojí před volbou, jaký typ zvolit pro uložení svých čísel v databázi či jinde, je informován o nabízených možnostech (např. databázový program Access) a zjistí, že jedna z možností je třeba jednoduchá přesnost na příjemných 4 bytech - s poznámkou: 7 dekadických cifer. Ihned si tuto informaci promítne do své problematiky:

• Nejobecněji se poloha na Zemi udává v zeměpisných souřadnicích: zeměpisná délka, zeměpisná šířka.

• Maximální absolutní hodnota je zde 180 stupňů.

• Na desetinnou (spíše necelou) část zbývají 4 platné cifry. Lze se tedy spolehnout na 0.0001 stupně.

• Pro zjednodušení: poloměr Země ve sférickém modelu je zhruba 6378 km. Obvod je tedy kolem 40 000 km a to je 360 stupňů. Jeden stupeň je tedy 111 km a tedy 0.0001 stupně je asi 11 metrů.

• Při uložení zeměpisné polohy ve stupních se tedy lze spolehnout na stovky metrů, na desítky už ne vždy.

• Závěr: při modelování na úrovni Alp by problémy být nemusely; při modelování na úrovni ČOV (čistírny odpadních vod) rozhodně ano.

Protože geoinformatika obecně a GISy (geografické informační systémy) zvlášť stojí a padají s lokalizací místa na zemském povrchu, na tomto místě se zkušený uživatel - geoinformatik zarazí a začne se o problematiku zajímat hlouběji. Přijde přibližně na to, co popisují následující odstavce.

Obraz necelých čísel v binárním hardware

Počítačový hardware je k srpnu 2003 (co bude za dva roky?!) založen na různých fyzikálních principech, svou podstatou vždy binárních (dvojkových, dvoustavových): doména zmagnetovaná - nezmagnetovaná, zrcadlo odráží - neodráží, tranzistor vede - nevede apod. Cokoliv v tomto prostředí tedy musí být zobrazeno pomocí známých "nul a jedniček" jako generalizací oněch dvou stavů. Zde se soustřeďme na necelá čísla.

Poznámka: Raději zde používáme termín "necelá" místo běžně používaného termínu "desetinná", tzn. mající desetiny, setiny, tisíciny atd. Bude totiž řeč o zobrazení ve dvojkové soustavě, a jaký pak je binární ekvivalent slova desetinná - slovo označující "mající poloviny, čtvrtiny, osminy" snad neexistuje.

Modelů. kterak necelé číslo zobrazit, může být celá řada. V podstatě každý si může vymyslet svůj vlastní. Důležité však je, co v této oblasti vymyslel výrobce té součástky v počítači, která jako jediná necelá čísla zpracovává - a tou je (v počítačích označovaných jako PC) tzv. matematický koprocesor. Koprocesor je počínaje vyšší řadou procesorů 486 integrován v čipu procesoru, není tedy samostatnou jednotkou jako dříve ve dvojici procesor 386 - koprocesor 387.

Poznámka: Zpracovávat hodnoty - tj. provádět definované operace s hodnotami - dovedou samozřejmě programové sekvence. Programy dovedou zpracovávat hodnoty jakéhokoliv logického formátu, ovšem vykonáváním posloupností strojových instrukcí nad elementárními typy dat, do kterých podprogramy zpracovávaný typ převádí. Zde však jde o zpracování přímo strojovými instrukcemi daného procesoru. Ty evidentně musí být konstruovány pro jistý pevně definovaný formát hodnot; a o něm je řeč.

Binární minimum

Nejmenší paměťovou jednotkou je byte (kdysi překládáno jako slabika, dnes se "bajt" většinou vůbec nepřekládá): posloupnost osmi dvojkových cifer = bitů. Jsou-li bity vyjádřeny jako 0 a 1, pak je jeden byte uspořádanou osmicí nul a jedniček. Těchto (různých) osmic je 256, počínaje např. 00000000, pak 00000001, ..., až po 11111110 a 11111111. Uvedené osmice lze chápat jako zápisy čísel ve dvojkové soustavě. Dvojková čísla uvedená shora (ve dvojkové soustavě) odpovídají následujícím číslům (v desítkové soustavě): 0, 1, ..., 254, 255. Graficky je možno jeden byte znázornit např. takto:
 

Byte:
Obsah	O	O	O	I	O	O	I	O
bitu č.	7	6	5	4	3	2	1	0

 

Byty nějakého paměťového prostoru jsou dnes přístupny linearizovanou adresou. Na paměť jako uspořádanou množinu bytů je možno pohlížet např. takto:
 

Paměť:
Obsah	OIIOIOIOO	OIIOIOIOI	OOIIIIIO	OIIOIOII	OIIIIOOIO	OIIIOOIO	IIIOOOOOI	...
bytu s adresou:	0	1	2	3	4	5	6	...

 

Hodnota je pak zobrazena na jednom, dvou nebo více bezprostředně po sobě jdoucích bytů. Adresa hodnoty je adresou jejího nejnižšího bytu.

Protože většina lidí je zvyklá chápat hodnoty vyjádřené zápisem čísla jen v desítkové soustavě, pro upřesnění způsob převodu zápisu nějaké (např. 1001) číselné hodnoty ze dvojkové do desítkové soustavy:
 

Konečně (pro diskutovanou problematiku nesmírně důležitý) zápis necelého čísla - např. 1001,011₂:
 

Celá část

,

"Necelá" část

celkem

23 = 810

22 = 410

21 = 210

20 = 110

2-1 = 1/210 = 0,510

2-2 = 1/410 = 0,2510

2-3 = 1/810 = 0,12510

x 1

+

x 0

+

x 0

+

x 1

+

x 0

+

x 1

+

x 1

8

+

0

+

0

+

1

+

0

+

0,25

+

0,125

=

9,37510

Důležitost spočívá v tomto faktu: nelze konečným zápisem necelého čísla ve dvojkové soustavě zapsat přesně jednu desetinu (v desítkové soustavě, tj. 0,1₁₀)!

Celá čísla

Přestože se článek týká necelých čísel, je nutno popsat především zobrazení celých čísel. Ta - na rozdíl od necelých - zpracovávají přímo procesory, není zapotřebí koprocesorů. Dále: celá čísla jsou zobrazena vždy přesně. Uveďme princip zobrazení celých čísel se znaménkem na dvou bytech - analogicky se v procesorech Intel používá zobrazení na jednom, čtyřech a osmi bytech.

Datový typ často označovaný jako integer (také celé číslo se znaménkem) je takové chápání obsahu dvou po sobě jdoucích bytů, při kterém jsou nazírány jako patnáctice dvojkových cifer - bitů předcházených jedním bitem nazíraného jako znaménko (nula jako znaménko +, jednička jako znaménko -). Patnácticiferný zápis dvojkové hodnoty tvoří nezáporné celé číslo; takových různých hodnot je 32 768, proto hodnotami jsou všechna čísla od 0 do 32 767 . Hodnotami datového typu integer jsou všechna čísla (desítkově) od -32 768^* do +32 767. Graficky lze jednu hodnotu typu integer znázornit např. takto:

 

Celé číslo se znaménkem:
Obsah	O	O	O	O	O	I	O	O	O	O	O	I	O	O	O	I
bitu č.	15 = znam.	14	13	12	11	10	9	8	7	6	5	4	3	2	1	0

 

Druhým způsobem zobrazení (avšak pouze nezáporných) celých čísel je ten, kdy se všechny bity chápou jako významové cifry čísla - tedy včetně "prvního zleva".

^*) Nejmenší hodnota by správně měla být hodnota -32 767. Při zobrazení doslova popsaném shora tomu tak skutečně je. Pak se ovšem nula v takové množině hodnot vyskytne jednou jako "kladná" nula, jednou jako "záporná" nula. Nula je ovšem nulou, ať má jakékoliv znaménko. Proto hardware počítačů (procesory Intel) umí šetřit, záporná čísla jakoby "o 1 posunou" a získají tak jednu další číselnou hodnotu navíc: právě onu hodnotu -32 768. Navíc - z důvodů čistě hardwarových - mají záporná čísla "prohozené" všechny nuly a jedničky; toto zobrazení se nazývá dvojkový doplňkový kód.

Necelá čísla

Protože počítače dvacátého století vznikly původně proto, aby počítaly (zvláště tak hanebnou věc jako je atomová bomba), je problematika necelých čísel v počítačovém prostředí hodně stará. Kupodivu způsob jejich uložení (tj. jejich logický formát) je až na nepodstatné drobnosti stále stejný. Procesory typu Intel pracují s logickým formátem definovaným ve standardu IEEE 754 (Institute of Electrical and Electronics Engineers).

Vychází se z toho, že každé necelé číslo C lze zapsat ve tvaru

C = mantisa x 10^exponent

Bývá zvykem označovat takový zápis jako semilogaritmický zápis číselné hodnoty.

Dále se bere v úvahu, že každé nenulové číslo lze (vynásobením vhodnou mocninou základu) zapsat tak, že celá část mantisy je jednociferné číslo různé od nuly. Např.

234,567 x 10¹² = 2,34567 x 10¹⁴

Takový zápis se nazývá normovaný semilogaritmický zápis.

Výrobci hardware berou uvedené skutečnosti v potaz, ale pracují zásadně ve dvojkové soustavě s normovaným semilogaritmickým zápisem. Necelá čísla jsou tedy zapsána ve tvaru

C = mantisa x 2^exponent

kde jak mantisa, tak exponent jsou opět dvojková čísla. Pro nenulové číslo tedy platí, že jeho normovaný semilogaritmický tvar má ve dvojkové soustavě před "dvojkovou čárkou" vždy cifru jedna - protože to je jediná nenulová cifra dvojkové soustavy. Protože tato jednička tam je vždy, nemusí být fyzicky zapsána (je jen "myšlena") a tím se ušetří jeden bit. Příklad zobrazení čísla 4,125₁₀ = 100.001₂ = (1.00001 x 2¹⁰)_{vše 2} na dvou bytech by mohl tedy být tento:

Přesnost zobrazení je zcela zřejmá, ale jen ve dvojkové soustavě: na tolik dvojkových cifer plus jedna, kolik je bitů části pro mantisu.

Poznámka: Exponent je celé číslo se znaménkem. Je uložen ve formátu tzv. posunuté nuly na rozdíl od dvojkového doplňkového kódu používaného pro celá čísla popsaná shora.

Jednoduchá přesnost

Zobrazení datového typu označovaného jako single (také krátké desetinné číslo nebo číslo v jednoduché přesnosti) je takové chápání obsahu čtyř po sobě jdoucích bytů (= 32 bitů), při kterém je 31. bit považovaný za znaménko čísla, bity 30-23 za exponent a bity 22-0 za mantisu. Přesnost zobrazení je 24 platných dvojkových cifer (přesně), což odpovídá uváděným 7 platným cifrám dekadickým. 

U procesorů Intel je exponentová část e chápána jako celé číslo bez znaménka. Skutečná hodnota H čísla se pak určí takto (znaménko čísla udává 31. bit, dále jen pro absolutní hodnotu):

• pro 0<e<255 je H = 1.m x 2^e-127

• pro e=0 a m # 0 je H = 0.m x 2^-126

• pro e=0 a m=0 je H = 0

Z toho plynou orientační hodnoty: číslo s největší absolutní hodnotou je přibližně 3,4 x 10³⁸, číslo s nejmenší absolutní hodnotou ještě různé od nuly je přibližně 1,4 x 10^-45.

Dvojnásobná přesnost

Zobrazení datového typu označovaného jako double (také double precision, long real nebo dlouhé desetinné číslo, číslo ve dvojnásobné přesnosti) je takové chápání obsahu osmi po sobě jdoucích bytů (= 64 bitů), při kterém je 63. bit považovaný za znaménko čísla, bity 62-52 za exponent a bity 51-0 za mantisu. Přesnost zobrazení je 53 platných dvojkových cifer (přesně), což odpovídá asi 16 platným cifrám dekadickým.

U procesorů Intel je exponentová část e chápána jako celé číslo bez znaménka. Skutečná hodnota H čísla se pak určí takto (znaménko čísla udává 63. bit, dále jen pro absolutní hodnotu):

• pro 0<e<2047 je H = 1.m x 2^e-1023

• pro e=0 a m # 0 je H = 0.m x 2^-1022

• pro e=0 a m=0 je H = 0

Z toho plynou orientační hodnoty: číslo s největší absolutní hodnotou je přibližně 1,8 x 10³⁰⁸, číslo s nejmenší absolutní hodnotou ještě různé od nuly je přibližně 4,9 x 10^-324.

Rozšířená přesnost

Zobrazení datového typu označovaného jako extended (také extended precision, rozšířené desetinné číslo) je takové chápání obsahu deseti po sobě jdoucích bytů (= 80 bitů), při kterém je 79. bit považovaný za znaménko čísla, bity 78-64 za exponent a bity 63-0 za mantisu. Přesnost zobrazení je 64 platných dvojkových cifer (přesně), což odpovídá asi 20 platným cifrám dekadickým.

U procesorů Intel je exponentová část e chápána jako celé číslo bez znaménka a nepracuje v tomto případě s myšlenou 1 jako s celou částí mantisy; jediná cifra celé části je udána bitem 63 a bit 62 je první cifrou necelé části mantisy. Skutečná hodnota H čísla se pak určí takto (znaménko čísla udává 79. bit, dále jen pro absolutní hodnotu):

• pro 0<e<2047 je H = bit₆₃.bity_62...0 x 2^e-16383

• pro e=0 a m=0 je H = 0

Z toho plynou orientační hodnoty: číslo s největší absolutní hodnotou je přibližně 1,1 x 10⁴⁹³², číslo s nejmenší absolutní hodnotou ještě různé od nuly je přibližně 3,4 x 10^-4932.

Přesnost v geoinformatických výpočtech

Na rozdíl od obecně platných skutečností popsaných shora se nyní článek začíná soustřeďovat na specifika oboru geoinformatiky ve vazbě na hardwarové (Intel) i softwarové (Microsoft) prostředí počítačů třídy PC. Není určen profesionálním programátorům, ale aplikačním uživatelům, u nichž se v oboru předpokládají alespoň minimální znalosti programování. Článek používá jako programovacího nástroje dnes nejrozšířenějšího jazyka, kompilačně i interpretačně Microsoftem podporovaného Basicu. Díky poměrně zdařilé syntaxi jde o moderní, objektově orientovaný jazyk lehce zvládnutelný širokou škálou uživatelů. Lze se s ním setkat ve všech aplikacích, které mají mít naděli na komerční úspěch - počínaje produkty Microsoftu (Excel, Access, Word aj.) až po tak renomované softwarové producenty, jako je Autodesk ve svém Autocadu, Corel ve svém Draw i Photo, Golden Software ve svém geostatistickém Surferu atd.

Skutečná přesnost

Rozeberme nyní skutečnou přesnost např. uložení hodnoty 180^o v jednoduché přesnosti.

Jest

180 = 128 + 32 + 16 + 4 = 2⁷ + 2⁵ + 2⁴ + 2² = IOIIOIOO₂

přesně. V normovaném semilogaritmickém tvaru to je

I.OIIOIOO₂ x 2⁷

Je proto

e = 134₁₀ = IOOOOIIO₂,  m = OIIOI₂

Pro názornost by měla být mantisa psána včetně bezvýznamných nul zprava na celkový počet 23 bitů

m = OIIO IOOO OOOO OOOO OOOO OOO₂

aby vyniklo, že poslední nula je na místě mocniny základu 2^-16, tj. zaokrouhleno 1,526 x 10^-5 = 0,00001526. Polovina této hodnoty, 0,00000763, je tedy maximální ztrátou přesnosti při převodu zaokrouhlením reálného čísla z okolí hodnoty 180 na zobrazení v jednoduché přesnosti. Dvojkový hardware však desítkově zaokrouhlovat pochopitelně neumí, část přesahující zobrazitelné cifry "usekne", proto je maximální ztráta přesnosti skutečně 0,00001526. Nejbližší menší číslo než 180, přesně zobrazitelné, je 179,99998474. Ve sférickém modelu zeměkoule je to v přepočtu na délku poledníku přibližně o 1,694 [m] méně než 180 a to je maximální chyba zobrazení intervalu <-180;+180> reálných čísel do konečné množiny racionálních čísel přesně zobrazitelných popisovaným způsobem.

Kumulovaná chyba

Opakovanými výpočty se zaokrouhlovanými číselnými hodnotami chyba postupně narůstá. Klasický příklad je cyklické přičítání hodnoty ke kumulativní proměnné. Mějme v modulu Excelu např. takto definovanou funkci Souc (N-krát přičte hodnotu X):

Function Souc(N As Long, X As Single) As Single
   Dim i As Long, S As Single
   S = 0
   For i = 1 To N
      S = S + X
   Next i
   Souc = S
End Function

Použijme ji někde v listu Excelu; výsledek pro N=100 je zřejmý z následujícího obrázku:
 

Výsledek je zde větší než očekávaná hodnota 10, chyba je asi +2x10^-5 [%]. Podstatný je však následující příklad:

Výsledek je zde menší než očekávaná hodnota 100, chyba je asi -0,001 [%]. Pro zajímavost, při N=100 000 je chyba -0,014 [%], při N=1 000 000 je chyba +0,96 [%] a při N=10 000 000 je chyba už +8,8 [%]!

Poznámka: Různá znaménka chyb jsou dána optimalizačními algoritmy Excelu a interpretu Basicu. Absolutní hodnota chyb narůstá bez ohledu na její znaménko.

Z příkladu je zřejmé, že ani v tak renomovaném a poměrně povedeném produktu, jakým je Excel, nelze s chybou pracovat jako se systematickou a dále algoritmicky zpracovatelnou.

Uvedený příklad demonstroval přesnost výpočtů poměrně uměle - je zřejmé, že mnoho obdobných cyklů lze zapsat kombinací celočíselných (tj. přesných) a neceločíselných proměnných s větší eliminací chyb. Opakované výpočty jsou však použity v tak zásadních momentech, jakým je např. výpočet goniometrických a podobných funkcí pomocí jejich rozvoje v řady. Zde nezbývá, než spoléhat na erudici programátorů takových funkcí a nanejvýš otestovat, zda přesnost je únosná pro problém, který je aktuálně řešen. Právě o takový test se pokouší tento článek.

Přepočty souřadnic

Poslední předchozí odstavec naznačuje, v čem je právě v geoinformatice jádro problému. Veškerá geometrie v tomto oboru není ani sférická, ale bohužel "rotačně elipsoidální". Goniometrické funkce, eliptické integrály a další numericky náročné konstrukce jsou zde běžnými nástroji.

Nejobvyklejší operací tohoto typu jsou převody souřadnic mezi různými soustavami souřadnými. Soustavy souřadné jsou zavedeny na plochách, do kterých se promítá, často přímo projekcemi samotnými. Jedna z "přirozených" soustav souřadných je dána na zvoleném elipsoidu rovníkem a nultým poledníkem a je vyjádřena známou dvojicí [zeměpisná délka, zeměpisná šířka]. Z této zeměpisné soustavy souřadné do souřadného systému projekcí a naopak jsou převody nejčastější a tam také je kritické místo z hlediska přesnosti.

Zkoumání přesnosti lze založit na v praxi zcela logickém požadavku:

T^-1 ( T ( [x,y] ) = [x,y]

kde [x,y] je transformovaný bod, T je transformace a T^-1 je transformace k ní inverzní. Požadavek lze hovorově vyjádřit takto: transformace tam a hned zase zpět souřadnice nezmění. Jestliže však transformace T provedená v počítačovém prostředí s (obecně nepřesně uloženými) hodnotami [x,y] pomocí (obecně nepřesných) matematických funkcí odevzdá nepřesnou hodnotu, pak evidentně inverzní transformace (z hlediska aplikace matematických funkcí nikoliv inverzním algoritmem!) nemůže obecně tuto chybu nezvětšit.

Lze tedy klást dvě otázky:

• Nakolik se liší vypočítané transformované souřadnice T ([x,y]) od teoretických, tj. nezatížených chybou výpočtu?

• Nakolik se liší zpětně transformované transformované souřadnice T^-1 (T([x,y])) od původních [x,y]?

Na obě otázky dá vyčerpávající odpověď teorie chyb. Problém je však v tom, že její použití vylučuje absence informací o vnitřní konstrukci algoritmů knihovních matematických podprogramů používaných konkrétním software. Je ale možno se pokusit o ryze praktické zhodnocení výsledků transformací vlastním zápisem algoritmu a jeho uplatnění na známou dvojici souřadnic téhož bodu.

Praktické zhodnocení přesnosti transformací pro Geoinformatiku

Tento článek nehodnotí nějaké konkrétní programové systémy třídy GIS nebo Geostatistiky. Ty je totiž nutno přijímat tak, jak jsou. Článek diskutuje otázku přesnosti v geoinformatice obecně s přihlédnutím k hardware zvláště proto, že řada uživatelem požadovaných výpočtů (ať přípravných nebo finálních) ani žádným konkrétním software pokryta není. Nezbývá tedy než vytvořit software vlastní.

Odborníci v geovědách však nemusí být stoprocentní programátoři, alespoň minimální povědomí o vlastní programové podpoře by však mít měli. Nejméně bolestným nástrojem se ukazuje Basic pro aplikace tak, jak ho dnes nutí firma Microsoft všude a všem. Velkou výhodou pak je to, že i neprogramátor zvládne zařazení hotového modulu do své aplikace - a využívá pak funkcí modulu stejně jako ostatních funkcí své aplikace.

Nejčastější obecnou aplikací využívanou pro technické i jiné výpočty je v našich končinách bezesporu Excel. Právě do sešitů jím zpracovávaných lze vlastní programové moduly zařazovat velmi jednoduše, stejně jednoduché je i využívání připravených doplňků.

Zmíněné praktické zhodnocení přesnosti transformací je provedeno funkcemi v modulu, zařazeného do sešitu Excelu. Pro relativní složitost (a tedy možnost největších teoretických chyb) stejně jako pro nejčetnější výskyt v mapových dílech ČR byly zkoumány transformace mezi zeměpisnou soustavou souřadnou na Besselově elipsoidu a projekcí JTSK, známou jako "Křovák".

Typy neceločíselných dat zpracovávaných Basicem

V Basicu lze pro výpočty s necelými čísly použít principielně tři typy:

• Single

• Double

• Decimal

Zatímco Single a Double jsou typy popsané výše a jsou přímo deklarovatelné a zvláště přímo zpracovatelné procesorem, typ Decimal je možno použít jen jako subtyp typu Variant; přiřazovat hodnoty lze jen použitím konverzní funkce CDec. Typ Decimal odpovídá modelově hardwarovému formátu Packed Decimal (označované také jako BCD - Binary Coded Decimal) uloženým na 14 bytech. Basic pracuje s těmito hodnotami s přesností 28 platných cifer tak, že staví nad BCD aritmetikou procesoru Intel své vlastní podprogramy - proto tento typ není přirozeným typem, ale subtypem. Pro svou značnou přesnost je důležité se zabývat i tímto subtypem.

Projekce JTSK

Geoinformatické veřejnosti je dostatečně známa. Jde o kuželovou projekci v obecné poloze, matematická podstata je podána např. zde.

Kromě značného množství použití goniometrických funkcí je záludnost dána aproximacemi funkcí řadami. Ze shora uvedeného textu je patrné, proč tomu tak je. Sčítají se totiž hodnoty řádově velmi odlišné. V binární aritmetice necelých čísel popsané shora se operace sčítání provádí nad základem tvořeným větší ze sčítaných hodnot, a tedy přesnost menší ze sčítaných hodnot se dále zmenšuje o počet bitů daný zhruba rozdílem dvojkových exponentů normovaných semilogaritmických tvarů sčítanců. Existuje tedy člen řady, který ještě má vliv na změnu součtu, přičemž (strojovým!) přičtením následujícího členu se už součet nezmění. Kolik takových počátečních členů řady má na součet vliv, záleží na hardwarové velikosti zobrazení necelého čísla (viz výše).

Pro zkoumání přesnosti byla připravena jednak trojice funkcí s hlavičkou

Public Function be2kr(LambdaStup As ttt, FiStup As ttt) As ttt()

pro převod ze zeměpisných souřadnic na Křovákovy, a trojice funkcí

Public Function kr2be(Xkr As ttt, Ykr As ttt) As ttt()

pro převod opačný.

Symbolem ttt je označen jeden z typů Single, Double a Variant pro zpracování hodnot různých typů: Single, Double a Decimal.

Funkce používající typů Single a Double jsou na zápis jednoduché a v podstatě přepisují postup výpočtu. Např. fragmenty funkce pro typ Double podává následující výpis:

Public Function be2kr(LambdaStup As Double, FiStup As Double) As Double()

   ' Všechny úhly vstupují ve stupních. Všechny dálky vystupují v metrech.

   Dim X As Double
   Dim Y As Double

   . . .

   Const R = 6380703.6105# ' Poloměr koule, na kterou se konformě zobrazí Besselův elipsoid
   Const S0 = (78# + 30# / 60#) * gToRad ' Dotyková konstrukční rovnoběžka
   Const Fi0 = (49# + 30# / 60#) * gToRad ' Jediná délkově zachovaná rovnoběžka elipsoidu

   . . .

   Dim Fi As Double: Fi = FiStup * gToRad
   Dim Lambda As Double: Lambda = LambdaStup * gToRad
   Dim V As Double: V = (Lambda + Ferro) * Alfa
   Dim dF As Double: dF = (Fi - Fi0) * gToDeg

   . . .

   X = Ro * cE
   Y = Ro * sE

   DelZkres = (n * Ro) / (R * cS)
   MeridKonv = E - ArcCos((sUk - sS * sU) / (cS * cU))

   ReDim Vysl(1 To 4)

   Vysl(1) = X
   Vysl(2) = Y
   Vysl(3) = DelZkres
   Vysl(4) = MeridKonv

   be2kr = Vysl

End Function

Pro tytéž výpočty nad subtypem Decimal je však zapotřebí volit jiné obraty. Zachování přesnosti totiž teoreticky zaručují pouze přímé binární operace nad subtypem Decimal. Proto zápis algoritmu je zdlouhavější:

Public Function kr2be(Xkr As Variant, Ykr As Variant) As Variant()

   ' Všechny úhly vstupují ve stupních. Všechny dálky vystupují v metrech.

   Dim V1 As Variant
   Dim V2 As Variant
   ...
   Dim V8 As Variant
   Dim V9 As Variant

   . . .

   Dim n As Variant: n = CDec(Sin(S0))
   Dim D As Variant: D = CDec(E / n)

. . .

   V1 = CDec(Ro / R0)
   V2 = CDec(1 / n)
   V3 = CDec(V1 ^ V2)

   . . .

   V9 = CDec(k2 * V8)

   Dim S As Variant: S = CDec(V9)

   Dim sS As Variant: sS = CDec(Sin(S))
   Dim cS As Variant: cS = CDec(Cos(S))
   Dim sUk As Variant: sUk = CDec(Sin(Uk))
   Dim cUk As Variant: cUk = CDec(Cos(Uk))

   . . .

End Sub

Pro úplnost však byla připravena i funkce, která sice pracuje nad typem Decimal, ale používá obecný zápis výrazů analogicky funkcím pro Single a Double a ponechává na lexikálním analyzátoru jazyka, jak s nimi naloží. V následujících tabulkách je označen jako "typ dat" Dex, i když jde spíše o odlišný způsob programování.

Test projekce JTSK

Mějme bod na Besselově elipsoidu o souřadnicích [ 18^o 15', 48^o 50'] - jde o místo v Ostravě - Mariánských horách. V zobrazení JTSK má bod souřadnice [ X ; Y ] v jednotkách [m].

Vytvořme v modulu sešitu Excelu shora popsané funkce pro převod souřadnic. Jako zdroj parametrů připravme do nějakých buněk listu hodnoty souřadnic (X,Y jsou zde spočteny jako Double a zobrazeny jen na 4 desetinná místa):
 

Zapišme nyní 4 vzorce pro převod Besselových souřadnic na JTSK s parametry [Lam,Fi] pro každý z typů parametrů popsaných výše. Výsledkem jsou spočtené souřadnice označené [X',Y'] v JTSK. Další 4 vzorce pro převod souřadnic JTSK na Besselovy budou jako parametry používat [X',Y'] - jejich výsledek označme [Lam',Fi']. Konečně nakonec zapišme cyklus N-krát provádějící přepočet typu tam - zpět. Následující tabulky podávají výsledek takové dvojí transformace:

1. Spočtené souřadnice [X',Y']
 

2. Odchylka spočtených souřadnic od skutečných:
 

3. Zpětný přepočet na zeměpisné souřadnice:
 

4. Odchylka od výchozích souřadnic:
 

5. Převod tam - zpět 1000x:
 

Výsledek testu projekce JTSK

Při testech je především nutno brát v úvahu prvotní chybu při ukládání výchozích hodnot testu. Zatímco 18^o15' se uloží přesně při všech typech dat (=18,25^o = 18 + 1/2²), hodnota 49^o50' je vždy zaokrouhlena. Stejně tak je tomu se "skutečnými" souřadnicemi [X,Y], které jsou porovnány s výsledkem výpočtu, ale které už samy o sobě jsou uloženy nepřesně. Na druhé straně tyto zaokrouhlovací chyby jsou chybami na posledních max. 2 bitech mantisy.

Z tabulky uvedené jako bod 2. plyne jednoznačný závěr. Jedině při použití typu Single jsou odchylky nejen teoreticky významné, ale významné mohou být i při mnoha praktických aplikacích vyžadujících např. decimetrovou přesnost. Tento typ dat je ovšem často používaný v databázových aplikacích. Zabírá polovinu místa než Double a při dvojici souřadnic je to úspora 8 bytů na bod terénu. Aplikace typu GIS ovšem s určením bodu terénu stojí a padají a celkem běžných sto tisíc bodů pak reprezentuje téměř megabyte paměťového prostoru.

Použití ostatních typů dat dává výsledky se zcela zanedbatelnou odchylkou.

Při zpětném převodu (bod 4) jsou výsledky obdobné. Nezanedbatelnou odchylku vykazuje jen typ Single - asi 0,1" je na poledníku asi 3 metry. Použití ostatních typů dat dává výsledky se zcela zanedbatelnou odchylkou.

Zajímavé výsledky podává bod 5. Programově se vykonal N-krát převod tam - zpět. Tabulka v předchozím odstavci podává výsledky pro N=1000. Je dobré porovnat odchylku při 1000 převodech u typu označeného jako Dex s odchylkou při jediném převodu: jsou totožné. Větší odchylku vykazuje typ double, největší typ Single. Tento test byl proveden s různými N: 10, 100, 100 000. Odchylky pro typ označený jako Dex jsou odchylky stejné pro každé N. Odchylky pro typ označený jako Dec jsou odchylky stejné pro každé N>25. Odchylky pro typ označený jako Dbl jsou odchylky stejné pro každé N>280.

Vysvětlení této zajímavosti spočívá ve způsobu zobrazení desetinných čísel - viz shora. Souřadnice bodu tvoří dvojici desetinných čísel, všechny dvojice tvoří kartézský součin

R x R

kde R je konečná množina racionálních čísel zobrazitelných přesně v počítačovém prostředí. Graficky zobrazeno, jde o dvourozměrnou "mřížku" a každý reálný bod daný svými souřadnicemi se zaokrouhlí na nějaký uzel této mřížky. Do výpočtu pak vstupuje ne skutečný bod, ale uzel mřížky. Výsledkem výpočtu je opět nikoliv reálný bod, ale uzel mřížky. Jeho vzdálenost od skutečného očekávaného bodu je pak dána zaokrouhlovacími chybami operací použitých ve výpočtu. Po jistém počtu kroků typu tam - zpět se pak dosáhne toho, že všechny hodnoty, se kterými se ve výpočtu operuje, jsou už jen uzly mřížky a tedy dávají vždy stejné výsledky.

Další projekce

Analogicky, avšak již ne tak důkladně, byly zkoumány projekce UTM jak na WGS84, tak na Krasovkého elipsoidu ("Old Soviet"), prováděny přepočty z geocentrických na eliptické souřadnice ap. Výsledky byly zcela analogické, dokonce i řádově, přestože jde o relativně jednodušší výpočty než v JTSK. Autor tohoto článku se domnívá, že hlavní důvod je

• v tom, že k největším zaokrouhlovacím nepřesnostem dochází v počátky výpočtu (při "uchytávání na mřížku")

• v tom, že pro výpočet standardních funkcí použili autoři nejpřesnější typ, na nějž se argumenty rozšíří a jejichž výsledek se na typ argumentu zaokrouhlí.

Závěr

Ukazuje se, že v úvodu zmíněný zkušený geoinformatik si může oddechnout. Zdá se, že s přesností výpočtů to není nijak zlé. Ovšem velký pozor je zapotřebí dát na čtyřbytový typ Single a používat ho pouze v přesně daném kontextu budované datové základny. Článek ukázal na nebezpečí při jeho použití pro přímé uchování zeměpisných souřadnic. Tento typ však lze bez většího nebezpečí použít aplikacích místního rozsahu ne jako přímá souřadnice, ale jako offset - doplněk k jistému počátku; ten však už musí být minimálně s dvojnásobnou přesností.

Příklad: Zkoumáme Ostravu - Mariánské Hory, tj. okolí místa použitého shora v testu: [ 18^o 15', 48^o 50']. Toto místo nechť je počátek lokální soustavy souřadné a jeho souřadnice nechť jsou uchovány s dvojnásobnou přesností. Místo o 300 [m] na sever má pak skutečné souřadnice  [ 18^o 15', 48^o 50' 10"] nebo relativní  [ 0^o 0', 0^o 0'10"] - a tento doplněk k počátku už může být skutečně v jednoduché přesnosti. Ve vzorci obsahující součet Počátek + Doplněk dojde k rozšíření Single na Double a výsledek je pak ve dvojnásobné přesnosti. Tímto způsobem samozřejmě nelze zkoumat terén o rozloze stovek kilometrů.

Rev. 07 / 2003