ČESKÉ VYSOKÉ UČENÍ TECHNICKÉ V PRAZE FAKULTA STAVEBNÍ KATEDRA MAPOVÁNÍ A KARTOGRAFIE Studijní obor: GEODÉZIE A KARTOGRAFIE
VÝPOČTY V SYSTÉMU MGRS BAKALÁŘSKÁ PRÁCE
Vypracovala: Adéla Kratochvílová Vedoucí práce: Prof. Ing. Bohuslav Veverka, DrSc. 2007
CZECH TECHNICAL UNIVERSITY IN PRAGUE FACULTY OF CIVIL ENGINEERING DEPARTMENT OF MAPPING AND CARTOGRAPHY Branch of study: GEODESY AND CARTOGRAPHY
COMPUTATIONS IN THE MILITARY GRID REFERENCE SYSTEM BACHELOR WORK
Author: Adéla Kratochvílová Supervisor: Prof. Ing. Bohuslav Veverka, DrSc. 2007 2
Prohlašuji: Tuto práci jsem vypracovala samostatně. Veškeré literární prameny a informace, které jsem k práci použila, jsou uvedeny v seznamu použité literatury.
V Praze dne 10.8. 2007
Adéla Kratochvílová
3
Abstrakt: Tato bakalářská práce je zaměřena na výpočty v rámci hlásného systému MGRS (Military Grid Reference System). MGRS je alfanumerickým kódem rovinných souřadnic, který využívají armády NATO (North Atlantic Treaty Organisation). Práce představuje základní informace, potřebné pro sestavení algoritmu převodu, o systému WGS 84 (World Geodetic System 84) projekci UTM (Universal Transverse Mercator) a hlásném systému MGRS. Stěžejní částí celé práce je program na převod souřadnic mezi systémy UTM, WGS 84 a MGRS, který byl vpracován pomocí vývojového prostředí Microsoft Visual C++ 6.0.
Abstract: This bachelor thesis is focused on computations in Military Grid Reference
System.
MGRS
(Military
Grid
Reference
System)
is
an
alphanumeric character set for plane coordinate system, which is using by NATO (North Atlantic Treaty Organisation) army. The work present the basic information about system WGS 84 (World Geodetic System), projection UTM (Universal Transverse Mercator) and system MGRS which are required to compile algorithm of conversion. The basic part of this work is a programme for conversion among systems WGS 84, UTM and MGRS creating in development environment Microsoft Visual C++ 6.0.
Klíčová slova C++, elipsoid, MGRS, NATO, rovinné souřadnice, UTM, WGS 84 Key words: C++, ellipsoid, MGRS, NATO, plane coordinate, UTM, WGS 84 4
Poděkování: Chtěla bych poděkovat panu Prof. Ing. Bohuslavu Veverkovi, DrSc. a mému otci Ing. Miroslavu Kratochvílovi, CSc., za jejich trpělivost, odbornou pomoc a cenné rady k náležitostem mé práce a poskytnutí informací pro tvorbu výsledného programu.
5
Obsah: Použité zkratky
7
Použité symboly a značení
8
Úvod
9
1. NATO a kartografie
10
1.1 Kartografické standardy NATO 2. WGS 84 – Světový geodetický referenční systém 2.1 Charakteristika
10 11 11
2.1.1 Primární parametry
12
2.1.2 Sekundární parametry
12
2.1.3 Odvozené geometrické parametry
13
2.2 Systém souřadnic na elipsoidu WGS 84
14
2.3 WGS 84 v ČR
16
2.4 Přesnost WGS 84
17
2.5 Perspektiva využití WGS 84
17
2.6 Přechod UTM na WGS 84
17
3. UTM – Transverzální Mercatorovo zobrazení
19
3.1 Charakteristika
19
3.2 Popis zobrazení
19
3.2.1 Značení
20
3.3 Průběh zkreslení
20
3.4 Mercator
21
4. MGRS 4.1 Značení 5. Matematický úvod
25 25 28
5.1 Použitá symbolika
28
5.2 Hodnoty vybraných veličin
29
5.3 Matematický vztah pro UTM a elipsoid
29
6. Popis programu 6.1 Microsoft Visual C++
32 35
Závěr
37
Seznam použité literatury
38
Seznam příloh
38
6
Použité zkratky: •
AČR … Armáda České Republiky
•
BIH … Bureau International de I´Heure
•
CTP … konvenční terestrický pól (Conventional Terestrial Pole)
•
CTRS … konvenční terestrický systém (Convention Terestrial Reference System)
•
DMA (dnes NIMA) … Obranná mapovací agentura armády USA (Defense Mapping Agency Æ dnes National Imagery and Mapping Agency)
•
EGM96 … Gravitační model Země (Earth Gravity Model 1996)
•
ETRF89 … Europian Terestrial Reference Frame 1989
•
ETRS … Europian Terestrial Reference System
•
EUREF-CS/H-91, CS-NULRAD-92, CS-BRD-93, DOPNUL, VGSN-92, VGSN-99 … GPS kampaně na území ČR
•
GIS ... Geografický informační systém
•
GPS … Global Positioning System
•
OCS ... Operational Control System , kontrolní (pozemní) segment GPS
•
IFOR ... Implementation Force, mírová operace NATO
•
ITRF (ITRS) … International Terestrial Reference Frame (System)
•
MGRS … Military Grid Reference System
•
MNČ… metoda nejmenších čtverců
•
NASA GSFC … NASA Goddar Space Flight Center
•
NATO ... Severoatlantická aliance (North Atlantic Treaty Organization)
•
NNSS … Námořní navigační družicový systém (Navy Navigation Satellite System)
•
S-42 ... Souřadnicový systém 3°a 6°pásů Gaussova zobrazení poledníkových pásů Krasovského elipsoidu
•
SFOR ... Stabilization Force, mírová operace NATO
•
SIS ... Státní informační systém
•
TRANSIT ... navigační družice, poprvé vyslána na oběžnou dráhu r. 1960 pro zájmy amerického námořnictva
•
TS ... Topografická služba
7
•
UPS … Universal polar Stereographic (Polární stereografické zobrazení)
•
UTM … Universal Transverse Mercator (Mercatorovo transverzální zobrazení)
•
VGIS ... Vojenský geografický informační systém
•
VTIS ... Vojenský topografický informační systém
•
VÚGTK … Výzkumný ústav geodetický, topografický a kartografický
•
WGS 84 … World Geodetic System 1984 (Světový geodetický systém z r. 84)
Použité symboly, značení •
ϕ (= B), λ (= L)… zeměpisné souřadnice na elipsoidu (něm. B – Breite, L – Lange, angl. Ltatitude, Longitude)
•
N,E … (Northing, Easting)rovinné souřadnice v zobrazení UTM
•
mB, mL ... střední kvadratická chyba souřadnic (šířky, délky)
•
mh ... střední kvadratická chyba geodetické výšky (nad elipsoidem)
8
Úvod Práce se zabývá převodem zeměpisných a rovinných souřadnic do hlásného systému MGRS. Tento hlásný systém patří mezi kartografické standardy NATO a byl přijat ČR. Hlavním cílem bylo vytvořit funkční a uživatelsky přívětivý program, který by byl schopen obousměrného převodu souřadnic mezi elipsoidem WGS 84 (φ, λ) a zobrazení UTM (N, E) a jejich vyjádření v alfanumerickém kódu MGRS. Pro sestavení algoritmů vedoucích k těmto výpočtům je nutné znát definiční parametry jednotlivých systémů a zobrazení. Veškeré charakteristiky jsou uvedeny v teoretické části práce, včetně krátkého textu o Severoatlantické alianci (NATO) a historického úvodu týkajícího se kartografa Mercatora, po němž je zobrazení UTM pojmenováno. V matematické části je uveden použitý algoritmus a matematický postup pro výpočet převodu souřadnic mezi uvedenými systémy. V závěru jsou uvedeny přílohy v podobě grafického znázornění značení projekce UTM. Práce má za úkol osvětlit základní pojmy matematické kartografie a dát je do souvislostí, které byly pro tento typ výpočtů použity. Hlavním cílem práce byl funkční výpočetní program, na základě dosažených znalostí v rámci výuky na fakultě a osobních možností a schopností se dále rozvíjet. Seznam použitého software a důvod jeho výběru je součástí textové části stejně jako popis programu.
Adéla Kratochvílová
9
1. NATO a kartografie v ČR V roce 1999 byla Česká republika přijata do vojenského svazku NATO (North Atlantic Treaty Organization), což se v dlouhodobějším časovém horizontu dotkne i vojenských kartografů, resp. Topografické služby Armády ČR (AČR). Ta bude přizpůsobovat svá mapová díla novým standardům NATO. Podle vládního nařízení 116/1995 Sb. jsou i vojenské mapy státním mapovým dílem, přičemž ve své základní verzi jsou dílem veřejným, tj. neutajovaným. Lze předpokládat, že se budou s kartografickými standardy NATO setkávat i civilní uživatelé těchto map. Navíc zřejmě sehrají tyto standardy roli sjednocujícího faktoru v oblasti geodetických a kartografických základů a ovlivní i vývoj civilní kartografie v oblasti státních mapových děl. (Civilní zeměměřiči využívají pro svoji práci méně přesný souřadný systém SJTSK, zatímco vojenští topografové využívají systém S-42/83, jehož rámcem je Astronomicko-geodetická síť (AGS), která se svojí kvalitou zaujímá přední místo ve světě). [1]
1.1 Kartografické standardy NATO Za standardy NATO lze považovat následující: •
Geodetické datum ED 50 (European Datum 1950, referenční plocha Hayfordův elipsoid), s použitím zejména v Západní Evropě či WGS 84 (referenční plocha elipsoid WGS84), s použitím zejména v technologiích sběru dat na bázi GPS (Global Positioning System),
•
Mercatorovo příčné válcové konformní zobrazení (UTM) v šestistupňových poledníkových pásech (viz obrázek), [1]
10
2. WGS 84 – Světový geodetický referenční systém WGS84 je světový geodetický geocentrický systém armády USA, ve kterém pracuje globální systém určování polohy GPS a který je zároveň standardizovaným geodetickým systémem armád NATO (STANAG 2211IGEO, ed.5 - Geodetic Datums, Ellipsoids, Drids and Grid References, 1991). 2.1 Charakteristika WGS84 je konvenční terestrický sytém (CTRS), realizovaný na základě modifikace Námořního navigačního družicového systému NNSS. Modifikace spočívá v posunu počátku souřadnicové soustavy, rotaci a změně měřítka dopplerovského systému NSWC 92-2 tak, aby byl systém geocentrický a referenční nultý poledník byl identický se základním poledníkem definovaným BIH (tento nultý poledník je posunut asi 100m východně oproti tradičnímu poledníku Greenwich). Definice systému WGS se vyvíjela od počátečního WGS 60, přes následující WGS 66, WGS 72 a WGS 84. V devadesátých letech byl systém WGS 84 zpřesněn dvakrát, poprvé v roce 1994 a podruhé v roce 1996. Tyto zpřesněné systémy byly označeny WGS 84(G730) a WGS 84 (G873), kde písmeno G značí zkratku GPS a číslo 873 znamená pořadové číslo týdne od zahájení provozu GPS, ve kterém byla do užívání zavedena zpřesněná varianta geodetického systému WGS 84. Přesná data, kdy byly tyto zpřesněné systémy zavedeny do GPS OCS, jsou 29.červen 1994 a 29. leden 1997. Společně s těmito zpřesněními pracovali na společném projektu NIMA, NASA GSFC a Státní univerzita v Ohiu. Výsledkem tohoto projektu je nový gravitační model Země: Earth Gravitational Model 1996 (EGM96).
EGM96
je
definován
Stokesovými
koeficienty
sférického
harmonického rozvoje tíhového potenciálu. Stokesovy koeficienty jsou vypočteny do stupně n = 360 a řádu k = 360. Využitím parametrů modelu gravitačního pole WGS 84 nazývaného zkráceně EGM96 je možné například vypočítat k bodu, u něhož známe souřadnice, výšku geoidu (tím i nadmořskou výšku). Také ale můžeme dostat tíhové zrychlení nebo tížnicové odchylky. [2] Systém WGS84 je definován fyzikálně. Jako normální zemské těleso byl zvolen hladinový rotační elipsoid, jehož základní parametry, velikost, tvar, hmotnost a rychlost rotace, byly odvozeny z pozorování, zejména pak 11
z pozorování družic. Známe-li 4 základní parametry referenčního elipsoidu, můžeme z nich odvodit vše potřebné: tíhový potenciál a složky tíhového pole, jakož i odvozené parametry. 2.1.1 Primární parametry definují rozměry referenčního elipsoidu přiřazeného systému. Jsou jimi: •
Velká poloosa referenčního elipsoidu a - jde o poslední geometrický parametr. Již delší dobu se předpokládá, že bude nahrazen jiným prvkem, prvkem fyzikálním, totiž hodnotou tíhového potenciálu W0 na referenčním elipsoidu. Důvodů je více, především však to, že hodnotu W0 lze odvozovat přímo z družicových měření a lze ji v přírodě realizovat, což u hodnoty a nelze.
•
Úhlová rychlost rotace vůči nebeskému referenčnímu systému ω .
•
Parametr mající vztah ke tvaru referenční plochy. Dnes je to převrácená hodnota zploštění referenčního elipsoidu 1/f . Tento parametr byl zvolen nedávno. Předtím se používal jiný prvek, např. čtverec prvé excentricity e2, dynamický tvarový faktor J2 nebo gravitační koeficient druhého stupně C2, 0 .
•
Parametr definující hmotnost referenčního tělesa M a gravitační konstantu G, což je geocentrická gravitační konstanta GM.
Definiční parametry jsou uvedeny v tabulce č.1. Tabulka č.1 Primární parametry systému WGS84 Název Symbol Hodnota Velká poloosa 6 378 137 a Převrácená hodnota zploštění 298,257 223 563 1/f Úhlová rychlost rotace Země 7,292 115 x 10-5 ω Zemská gravitační konstanta GM 3,986 004 418 x 1014 (vliv hmoty atmosféry započítán)
Rozměr m rad/s m3 . s-2 )*
)* Původní hodnota GM = 398 600,5 km3 s-2 byla nahrazena tabelovanou, aby se dosáhlo shody v definici se standardy IERS.
2.1.2 Sekundární parametry definují model struktury zemského gravitačního pole (Earth Gravity Model , EGM), definovaný pomocí rozvoje geopotenciálu do sférických harmonických funkcí do stupně a řádu 360. Model gravitačního pole WGS84 EGM96 je možné využít pro výpočet průběhu
12
plochy geoidu WGS84, tížnicových odchylek, středních hodnot tíhových anomálií v síti 10´x 15´s využitím tohoto modelu.
2.1.3 Odvozené geometrické parametry WGS 84 Základní parametry systému nestačí na všechny numerické výpočty geodetických úloh. Proto z nich odvodíme další hodnoty. Přitom se budeme řídit zásadou, že základní parametry byly stanoveny zcela přesně (to znamená, že za hodnotami uvedenými v tab.č.1 následují samé nuly) a odvozené parametry z nich vypočteme na tolik desetinných míst, kolik budeme pro naše účely potřebovat. Ve světové literatuře bylo odvozeno několik postupů dalšího zpracování. Čtverec prvé excentricity určíme ze vztahu e2 = 2 f − f 2 (1 ) kde
f =
1 1 f
.
Čtverec druhé excentricity (e´)2 je dán rovnicí
(e′)2 =
e2 1 − e2
(2 ) Pro čtverec lineární excentricity E2 platí
E 2 = a 2 − b2 , (3 ) kde hodnotou malé poloosy b dostaneme z výrazu b = a (1 − f ) .
(4 ) Starší definice referenčního elipsoidu používala ještě gravitační koeficient druhého stupně C2, 0 , respektive dynamický tvarový faktor
J 2 = − 5 ⋅ C´´2,0 . Pro jejich vzájemné vztahy platí rovnice
13
C
´ 2, 0
J =− 2 , 5
1 ⎡ 2 2ω 2 a 2e3 ⎤ J 2 = ⎢e − ⎥ 3⎣ 15GMq ⎦
(5 )
q=
kde
3 1 ⎡⎛ ⎢⎜⎜1 + 2 ⎢⎣⎝ (e′)2
⎞ 3⎤ ⎟arctg (e′) − ⎥ ⎟ e′ ⎥⎦ ⎠ [3]
Odvozené parametry jsou uvedeny v tabulce č.2 . Tabulka č.2 Odvozené parametry systému WGS 84 Člen f e2 e (e´)2 e´ E2 E b C2 , 0
J2 q
Hodnota Rozměr Rovnice 0,003 352 810 664 747 481 0,006 694 379 990 141 317 (1) 0,081 819 190 842 621 5 0,006 739 496 742 276 435 (2) 0,082 094 437 949 695 7 272 331 606 107,554 726 970 m2 (3) 521 854,008 423 385 330 012 m 6 356 752,314 245 179 497 564 m (4) -4 (5) -4,841 667 749 599 43 x 10 -3 1,082 629 821 257 28 x 10 (5) -5 7,334 625 786 725 72 x 10
Obr.1: Schéma geocentrického souřadného systému WGS 84, [3]
2.2 Systém souřadnic na elipsoidu WGS 84 WGS 84 používá souřadnice zeměpisné, jednotlivé body jsou definované zeměpisnou délkou, šířkou a výškou. Zápis zeměpisné šířky a délky může 14
být ve stupních, ve stupních a minutách nebo ve stupních, minutách a vteřinách (s jednotkami v desetinném tvaru). Je třeba vždy správně uvést, zda se jedná o severní nebo jižní šířku a o východní nebo západní délku. Zápis výšky je uveden v metrech - není to skutečná nadmořská výška, ale elipsoidická výška (neboli vzdálenost od elipsoidu). Počátek souřadnicové soustavy a směry souřadnicových os jsou definovány následovně:
•
počátek je umístěn do těžiště Země,
•
osa Z má směr ke konvenčnímu terestrickému pólu (CTP) definovanému BIH na základě souřadnic stanic definujících systém BIH,
•
osa X je definována průsečnicí referenčního poledníku WGS 84 a roviny rovníku vztažené k CTP, referenční poledník je nultý poledník definovaný BIH,
•
osa Y doplňuje systém na pravotočivý pravoúhlý souřadnicový systém, směr kladné části osy leží v rovině rovníku 90° východně vzhledem k ose X. Počátek souřadnicového systému WGS 84 je totožný se středem
elipsoidu WGS 84 a osa Z je rotační osou elipsoidu. Takto definovaný systém je s reálnou zemí spojen prostřednictvím pozemních stanic kontrolního segmentu GPS. Původní referenční rámec (z roku 1987) byl realizován prostřednictvím NNSS nebo TRANSIT (Doppler). [2]
Obr.2 : Souřadné systémy, [2]
15
2.3 WGS 84 v ČR Na území bývalého Československa bylo započato s realizací systému WGS 84 na základě kampaně VGSN'92, organizované DMA (Defense Mapping Agency = Obranná mapovací agentura armády USA, dnes NIMA National Imagery and Mapping Agency) a Topografickou službou (TS) Armády ČR (AČR). Od 1.1.1998 je WGS 84 zaveden ve vojenském a civilním letectvu a v AČR je běžně používán v rámci kooperace a armádami NATO a standardizace v geodézii a kartografii. Světový geodetický referenční systém 1984 je na území České republiky určen:
•
technologiemi kosmické geodézie, které jsou součástí programů monitorovacího a zpracovatelského centra správce systému,
•
souborem rovinných souřadnic bodů vztažených ke světovému geodetickému referenčnímu systému 1984 (World Geodetic System 1984), epoše G873,
•
elipsoidem světového geodetického systému 1984 s konstantami a = 6378137m, f = 1:298,257223563, kde "a" je délka hlavní poloosy a "f" je zploštění.“
Jiná definice uvádí, že geodetický systém WGS 84 je definován:
•
polohou počátku a orientace pravoúhlé prostorové souřadnicové soustavy,
•
parametry referenčního elipsoidu,
•
gravitačním modelem Země a geoidem Počátek a orientace souřadnicových os jsou prakticky realizovány
souřadnicemi X, Y a Z dvanácti stanic, které monitorují dráhy GPS družic. Souřadnice těchto monitorovacích stanic byly (do roku 1994) určeny na základě dopplerovských měření systému TRANSIT na základě zpracování časově dlouhodobých "absolutních" observačních kampaní. Od 1.1.1994 jsou WGS 84 souřadnice 10 sledovacích stanic GPS, zpřesněny na WGS 84 (G730) (Malys, Slater, 1994) a připojeny přesným relativním měřením pomocí technologie GPS k systému ITRF-91, později byl systém rozšířen na 12 stanic v dále zpřesněném systému WGS 84 (G873). [3]
16
2.4 Přesnost WGS 84 Přesnost geocentrických souřadnic bodů přímo určených v systému WGS 84 na základě technologie GPS, s použitím odpovídajících efemerid a relativního měření ve statickém módu, je charakterizována středními kvadratickými chybami v zeměpisné šířce B, zeměpisné délce L, a geodetické výšce h: mB = mL < 0,4 m mh < 0,5 m. Tyto chyby neobsahují pouze měřické chyby, ale především chybu v realizaci počátku souřadnicového systému (v současnosti cca 10 cm v každé souřadnici) a v určení rozměru sítě. [3] 2.5 Perspektiva využití WGS 84 Koncem 90. let byly péčí TS AČR geodetické polohové základy převedeny ze společného systému ETRS-89 do WGS 84, který je nyní využíván v AČR ke:
•
geodetickému zabezpečení letišť, civilních i vojenských,
•
geodetické lokalizaci prvků VTIS (Vojenský Topografický Informační Systém) a VGIS (Vojenský Geografický informační systém), které jsou podsystémem armádního SIS (Štábní informační systém AČR),
•
zabezpečení jednotek AČR, působících ve svazku IFOR, SFOR,
•
tvorbě mapového standardizovaného díla v zobrazení UTM
2.6 Přechod UTM na WGS 84 Na elipsoid WGS 84 je podle definice kartografického zobrazení UTM nasunut válec tak, aby jeho podélná osa ležela v rovině rovníku; průměr tohoto válce je poněkud menší než je řez elipsoidem WGS 84. Postup:
•
z plochy elipsoidu se vlevo 3° a vpravo 3° od osového poledníku zobrazí na plochu válce obraz zemského povrchu
•
válec se pootočí okolo osy rotace elipsoidu o 6° a opět se vše zobrazí na válec
17
•
tak se pokračuje celkem 60x, až je veškerý povrch Země zobrazen na válec
•
válec se „rozvine“ a je tak získáno 60 šestistupňových pásů na nichž je zobrazena v konformním zobrazení příslušná část povrchu Země, každý z těchto pásů má vlastní souřadnicovou – rovinnou soustavu – osa E je totožná s obrazem rovníku, osa N pak se středovým poledníkem; aby se vyloučilo střídání znamének v kvadrantech, odsune se osa N západním směrem o 500 km – pak jsou všechny souřadnice E kladné. [5]
18
3. UTM – Universal Transverse Mercator Univerzální transverzální Mercatorovo zobrazení je konformní válcové zobrazení v příčné poloze v 6° poledníkových pásech. Bylo vytvořeno roku 1772. Rovnice vznikly přenásobením rovnic zobrazení transverse Mercator modulem m = 0,9996. Tím se docílilo snížení vlivu zkreslení na okrajích pásů na 14 cm/km.
Obr.3: Transverzální válcové zobrazení, [5]
3.1 Charakteristika V roce 1947 bylo převzato americkou armádou pro vojenské mapy velkých měřítek. V současné době se využívá převážně pro potřeby NATO. Jako referenční elipsoid se nejčastěji používá Hayfordův elipsoid (pro NATO) a WGS 84. Od 1.1.2006 je zobrazení UTM na elipsoidu WGS 84 se souřadnicovým systémem WGS 84, používané v NATO, zavedeno též Armádě ČR. Nahradilo principem podobné Gauss – Krügerovo zobrazení na Krasovského elipsoidu se souřadnicovým systémem S - 42, zavedené v bývalé Varšavské smlouvě. 3.2 Popis zobrazení Rovník a poledníky každých 90o od základního poledníku (včetně) se zobrazují jako přímky. Ostatní poledníky jsou komplexní křivky konkávní vzhledem k základnímu poledníku. Ostatní rovnoběžky jsou konkávní vzhledem k nejbližšímu pólu. Zobrazení je symetrické podle rovníku a podle každého poledníku, který je zobrazen jako přímka. Póly se zobrazí jako body na základním poledníku.
19
3.2.1 Značení Referenční elipsoid je rozdělen na 60 očíslovaných pásů (zón), každý má šířku 6°. Zóny se číslují od poledníku 180° východním směrem od 1 do 60. Každá zóna má definován střední poledník, který se zobrazuje jako přímka. Např. zóna 1 se rozkládá od 180°do 174°z.d., střední poledník je 177°z.d.. Kvůli vzrůstajícímu zkreslení při postupu směrem k pólům je UTM definováno od 80°j.š. do 84° s.š.. Pro polární oblasti se používá speciální zobrazení UPS. V UTM je každý poledníkový pás samostatně zobrazen do roviny. Pás má vlastní soustavu pravoúhlých souřadnic vodorovnou osou je obraz rovníku a je kladná na východ, svislou střední poledník daného pásu a má kladný směr na sever. Abychom se vyhnuli záporným hodnotám souřadnice E v 2.a 3. kvadrantu (dle matematického levotočivého číslování kvadrantů) v rámci pásu, je svislá osa rovnoběžně přesunuta na západ o 500 km. Ze stejného důvodu se pro zobrazení jižní polokoule přičítá konstanta 10 000 000 m k souřadnici N. Poloha bodu o zeměpisných souřadnicích [ϕ, λ] je tedy v UTM vyjádřena číslem pásu a dvojicí pravoúhlých souřadnic [E, N]. Rovnoběžkové pásy jsou označeny písmeny A až Z (bez písmen O a I – ta jsou vynechána, aby nedošlo k záměně s čísly), jelikož má každý šestistupňový pás svou soustavu souřadnic, není dvojice hodnot E, N jednoznačným konkrétním vyjádřením polohy. Z tohoto důvodu se k předávání informací o poloze používá tzv. hlásný systém MGRS (viz kapitola 4), což je pouze jiný zápis prvků N a E. 3.3 Průběh zkreslení Zkreslení je v daném bodě ve všech směrech stejné (kružnice zkreslení v originále se zobrazí jako kružnice zkreslení i v mapě). Zkreslení je konstantní na základním poledníku. Dále plošné zkreslení roste směrem od základního poledníku, ale ve stejné vzdálenosti od tohoto poledníku je konstantní. Zkreslení je rovnoměrně rozloženo do plochy daného pásu, na středním poledníku má hodnotu 0,9996. [6]
20
Obr. 4: Rozsah poledníkového pásu UTM, [6]
3.4 Mercator Gerardus Mercator se narodil 5. března 1512 jako Gerhard Kremer ve vlámském městě Rupelmonde (dnes v Belgii, v té době část Holandska). Pouze dvacet let před jeho narozením Kryštof Kolumbus podnikl svoji historickou plavbu do Nového Světa, což mimo jiné vedlo mladého Kremera k myšlenkám o nových zeměpisných objevech. V roce 1530 se přihlásil na Univerzitu v Louvainu a po absolvování se vypracoval do postavení jednoho z hlavních evropských kartografů. V té době bylo obvyklé, že si studovaní mužové polatinštili svá jména. V doslovném překladu je holandské slovo kramer převedeno na merchant, česky nejspíš kupec nebo obchodník. Tedy Gerhard Kremer změnil své jméno na Gerardus Mercator a pod tímto jménem se i proslavil. Mercatorova nadějná kariéra byla ohrožena v roce 1544, kdy byl zatčen pro kacířství. Kvůli svému prosazování protestanství musel následně uprchnout do sousedního Duisburgu (dnešní Německo), kde zůstal až do konce svého života.
Obr.5: Geradus Mercator, [7]
21
Před Mercatorem zdobili kartografové své mapy nereálnými mytologickými figurkami a fiktivními zeměmi. Tyto mapy byly proto spíše uměleckým výtvorem než pravdivou prezentací světa. Mercator byl první, který své mapy stavěl výlučně na naprosto exaktních údajích. Ty získal díky průzkumům, a tak přeměnil kartografii z umění na vědu. Byl také první, u kterého se objevil svazek kolekcí samostatných map. Tento svazek byl nazýván atlas, u příležitosti legendární mytologické postavy, která drží glóbus. Tento umělecký výjev zdobil titulní stranu. Atlas byl publikován ve třech částech, poslední z nich spatřila světlo světa až v roce 1595 po Mercatorově smrti. Stalo se to roku 1568, kdy Mercator sám sobě uložil úkol věnovat svůj čas nové projekci mapy, která by odpovídala námořnickým potřebám a která by přeměnila celosvětovou navigaci z náhodného riskantního snažení na precizní vědu. Ze začátku byl nasměrován dvěma principy: mapa musí být rozprostřena do obdélníkové sítě, kde všechny kružnice představující body téže zeměpisné šířky musí být reprezentovány vodorovnými čarami, které jsou rovnoběžné s rovníkem a stejně dlouhé jako rovník. Druhým principem bylo, že mapa musí být izogonální, jelikož pouze a jen u izogonálních map je zachován skutečný směr mezi libovolnými dvěma body glóbu. Teď budeme mluvit o glóbu. Kružnice zeměpisné šířky zmenšují svoji velikost, když se příslušná zeměpisná šířka zvětšuje, a to až do dosažení nulové délky v kterémkoliv z obou pólů. Ale na Mercatorově mapě jsou tyto kružnice znázorněny vodorovnými čarami o stejné délce jako rovník. Tudíž každá rovnoběžka na mapě je vodorovně (to je z východu na západ) roztažena s jistým koeficientem závislým na zeměpisné šířce dané rovnoběžky.
Obvod kružnice se zeměpisnou šířkou φ je 2πr = 2πRcosφ na glóbu, zatímco na mapě je její délka 2πR . Kružnice je tedy roztažena s koeficientem
2πR = sec ϕ 2πR cos ϕ Aby byla mapa izogonální, s roztažením východo-západních rovnoběžek musí současně existovat i stejné severo-jižní roztažení vzdáleností mezi rovnoběžkami. A tohle severo-jižní roztažení se musí postupně zvětšovat, jak se blížíme k vyšším zeměpisným šířkám. Jinak řečeno, stupně
22
zeměpisné šířky, které jsou na glóbu stejně umístěné podél jednotlivých poledníků, musí být pozvolna zvětšovány na mapě. A to je klíčový princip této mapy.
Obr.6: Mercatorova zeměpisná síť [7]
Nicméně, aby byl tento plán uskutečněn, musí být nejdříve určeny vzdálenosti mezi sousedními rovnoběžkami. Jak přesně tohle Mercator dělal, není známo a stále se o tom vedou debaty v kartografii mezi historiky. Nenechal žádný psaný záznam své metody, až na následující vysvětlující dopis, který byl natištěn na jeho mapě: Při znázorňování světa musíme rozvinout povrch koule do roviny, a to takovým způsobem, aby se polohy míst na všech stranách navzájem shodovaly ve skutečném směru a ve skutečné vzdálenosti... S tímto úmyslem musíme použít nové proporce a nové sestavení poledníků s přihlédnutím k rovnoběžkám... Kvůli tomuto důvodu jsme postupně zvětšovali stupně zeměpisné šířky směrem k jednotlivým pólům přímo úměrně s prodlužováním rovnoběžek do délky rovníku. I toto nejasné vysvětlení nám dává jasně najevo, že Mercator musel dobře pochopit matematické principy, které jsou základem jeho mapy. Když vytvořil síť rovnoběžek a poledníků, zůstalo mu už jen položit kůži na kostru neboli překrýt svoji síť siluetou kontinentů, která byla známa v jeho době. V roce 1569 publikoval svoji mapu světa pod názvem Nový a vylepšený popis zemí ve světě, upraveno a určeno pro užívání navigátorů. Byla to ohromná mapa, která byla kvůli vytištění rozdělena na 21 částí a celkově měřila 54*83 palců. Je to jeden z nejvzácnějších kartografických artefaktů v historii. Jsou známy pouze tři exempláře originálu, které přetrvaly do současnosti. Mercator zemřel 2. prosince 1594 v Duisburgu. Žil dlouhý život, který mu přinesl slávu a bohatství. A přesto jeho proslulý úspěch - mapa, která nosila jeho jméno - nebyla hned přijata námořní komunitou, která nedokázala
23
pochopit nadměrné překroucení tvaru kontinentů. Je fakt, že Mercator nedodal úplné vysvětlení, kterak postupně zvětšoval vzdálenost mezi rovnoběžkami. Tím způsobil zmatek v myslích kartografů, kteří na jeho dílo chtěli navázat. [7]
24
4. MGRS – Military Grid Reference System Hlásné systémy NATO Hlásné systémy představují metodu určování polohy bodů na kartografickém obrazu zemského povrchu, která byla původně vyvinuta pro vojenské účely. Hlásný systém tvoří osnovy stejně vzdálených rovnoběžek, které na kterékoli mapě určují čtvercovou síť. Každý hlásný systém je dále definován klíčem, dle kterého je každému bodu v rámci sítě jednoznačně přiřazen alfanumerický kód. Mezi nejznámější hlásné systémy patří MGRS a GEOREF. Hlásný systém MGRS byl poprvé pospán americkou armádou v roce 1947. Jedná se o alfanumerickou verzi sytému UTM (nepolární oblasti) a UPS (polární oblasti), standardní pro mapy NATO. Umožňuje jednoznačnou identifikaci bodu kdekoli na Zemi. [1] 4.1 Značení Jak už bylo řečeno výše, jelikož má každý šestistupňový pás svou soustavu souřadnic, není dvojice hodnot E, N jednoznačnou identifikací polohy a proto se používá hlásný systém MGRS. Základním prvkem hlásného systému je tzv. zóna (obraz sférického čtyřúhelníku referenčního elipsoidu). Systém rozděluje každý poledníkový pás na 19 vrstev o šířce 8° a 1 vrstvu o šířce 12°, a to od rovnoběžky 80° jižní šířky (v algoritmu uvažujeme –80°) do rovnoběžky 84° severní šířky. Tak je zemský povrch rozdělen na 60 x 20 sférických čtyřúhelníků (ve skutečnosti vzhledem k nepravidelnosti dělení oblasti Špicberků je sférických čtyřúhelníků o 3 méně). Každé zóně je přiřazen hlásným systémem kód. První složka kódu je číslo zobrazeného poledníkového pásu od 1 do 60. Pásy se číslují od obrazu poledníku 180° západní délky, kdy v algoritmu uvažujeme –180°, směrem na východ. Druhá složka kódu je písmeno anglické abecedy C až X (písmena I a O jsou vynechána, aby nedošlo k záměně s číslicemi), které označuje vrstvu (vrstvy se značí od obrazu rovnoběžky 80° jižní šířky směrem na sever). Každá zóna je tedy kódem jednoznačně identifikovatelná.
25
Obr. 7: Souřadný systém MGRS v ČR [1]
V rozměrech zón jsou nepravidelnosti:
•
pás číslo 32 je mezi rovnoběžkami 56°a 64°s.š. (vrstva V) rozšířen o 3°pásu č. 31, tj. na výsledných 9°. Tímto opatřením zůstává jihozápadní část Norska v jednom poledníkovém pásu
•
pásy 33 a 35 jsou mezi rovnoběžkami 72°a 84° kvůli zachování Špicberků ve 2 pásech rozšířeny na 12°
•
tyto nepravidelnosti se kompenzují zvětšením pásů 31 a 37 ze 6°na 9°a vynecháním pásů 32, 34 a 36 mezi těmito rovnoběžkami Hlásnou síť dále tvoří čtverce o straně 100 km. Každý zobrazený
poledníkový pás je rozdělen systémem rovnoběžných čar s obrazem rovníku a příslušného osového poledníku (tj. s osou N a osou E). Obraz poledníkového pásu tedy pokrývá čtvercová síť (obrázek vlevo). V oblasti rovníku, vzhledem k šířce šestistupňového poledníkového pásu (668 km), takto vznikne šest úplných čtverců a na každém okraji pásu jeden neúplný čtverec, který má šířku 34 km. Se zužováním obrazů poledníkových pásů směrem k pólům se snižuje počet úplných čtverců a mění šířka okrajových čtverců. Každý čtverec (úplný i neúplný) se označuje dvěma písmeny, z nichž první písmeno popisuje sloupec (sloupce se značí od obrazu poledníku 180° směrem na východ) čtvercové sítě, ve kterém se nachází příslušný čtverec. Druhé písmeno určuje řádku (řádky se označují od obrazu rovníku směrem na sever a od obrazu rovníku směrem na jih), na níž se příslušný čtverec v rámci čtvercové sítě vyskytuje.
26
Sloupcům (včetně neúplných) jsou přidělena písmena anglické abecedy A až Z (bez I a O). Po písmenu Z se abeceda opakuje. Pro označení vrstev je využito písmen A až V (bez I a O). Po písmenu V se abeceda opakuje. U lichých poledníkových pásů začíná první vrstva od rovníku písmenem A, u sudých písmenem F. [6] Úplný kód polohy objektu v hlásné síti: označení zóny – číslo, písmeno, označení čtverce – písmeno, písmeno, pravoúhlé souřadnice bodu v rámci příslušného čtverce – posloupnost číslic. Každý čtverec určuje lokální soustavu souřadnic LS s počátkem v levém dolním rohu; první polovina posloupnosti číslic udává souřadnici W, tedy vzdálenost bodu od západní svislé strany čtverce, tj. od osy y v soustavy LS, druhá polovina souřadnici S, tedy vzdálenost bodu od jižní vodorovné strany čtverce, tj. od osy x v soustavy LS. Celý údaj (kód) se píše bez mezer a jakýchkoliv interpunkčních znamének následně: V rámci příslušného čtverce je poloha bodů určena posloupností číslic. Tyto číslice udávají pravoúhlé souřadnice lokální souřadnicové soustavy. Podle počtu číslic v souřadnici bodu lze poznat , s jak velkou přesností byla souřadnice určena [6]:
počet číslic 4 6 8 10
příklad 0491 042915 04259152 0425091520
přesnost (m) 1000 100 10 1
příklad kódu [2]:
27
5. Matematický úvod Základním úkolem matematické kartografie
je zobrazení povrchu
referenčního plochy (elipsoidu) do roviny. Vzájemné přiřazení polohy bodů na dvou různých referenčních plochách nazýváme kartografické zobrazení. Kartografická zobrazení rozlišujeme podle použitých referenčních ploch, v našem případě se jedná o zobrazení z elipsoidu do roviny. Kartografická zobrazení lze charakterizovat i z hlediska vlastností kartografických zkreslení. UTM je konformním zobrazením, zachovává tedy úhly. Základní souřadnicovou soustavou na referenčním elipsoidu jsou zeměpisné souřadnice, označované též geodetické zeměpisné souřadnice nebo
pouze
geodetické
souřadnice
(geographic
coordinate
system).
Souřadnice tvoří zeměpisná (geodetická) šířka φ (latitude) a zeměpisná (geodetická) délka λ (longitude). Zeměpisná šířka dosahuje hodnot v rozsahu <-90°, 90°>, často jsou tyto hodnoty označovány i jako jižní zeměpisná šířka (pro hodnoty <-90°, 0°>) a severní zeměpisná šířka (pro hodnoty <0°, 90°>). Zeměpisná délka používaná v běžném životě nabývá hodnot <0°, 360°> s počátkem na základním poledníku s přírůstkem ve směru východním. Čáry s konstantní hodnotou λ, resp. φ jsou nazývány zeměpisné poledníky (meridian), resp. zeměpisné rovnoběžky (parallel). Zeměpisné poledníky
a
rovnoběžky
vytvářejí
na
povrchu
referenčním
elipsoidu
zeměpisnou síť (graticule), která je při klasické tvorbě map důležitým konstrukčním prvkem při zobrazování povrchu elipsoidu do roviny. Zeměpisná síť umožňuje základní orientaci v obsahu map. 5.1 Použitá symbolika: x, y ... rovinné souřadnice poledníkového pásu (UTM) φ, λ ... zeměpisné souřadnice na elipsoidu k ... hodnota zkreslení e ... excentricita M ... meridiánová vzdálenost od rovníku k rovnoběžce o zem. šířce φ a ... velikost velké poloosy elipsoidu N ... příčný poloměr křivosti
28
5.2 Hodnoty vybraných veličin: a ... 6 378 137m k0...0,9996 e ... 0,081 819 190 842 621 5 e´2... 0,006 739 496 742 276 435 5.3 Matematický vztah pro UTM a elipsoid postup je popsán dle [8] ⎡ A3 A5 ⎤ + 5 − 18T + T 2 + 72C − 58e′ 2 x = k 0 N ⎢ A + (1 − T + C ) ⎥ 6 120 ⎦ ⎣
(
)
(1) ⎧ ⎡ A2 A4 A6 ⎤ ⎫ y = k0 ⎨M − M 0 + N ⋅ tgϕ ⎢ + 5 − T + 9C + 4C 2 + 61 − 58T + T 2 + 600C − 330e′2 ⎥⎬ 24 720⎦ ⎭ ⎣2 ⎩ (2)
(
)
(
)
⎡ A2 A4 A6 ⎤ k = k 0 ⎢1 + (1 + C ) + 5 − 4T + 42C + 13C 2 − 28e′ 2 + 61 − 148T + 16T 2 ⎥, 2 24 720 ⎦ ⎣
(
)
(
)
(3)
kde:
•
k0 ... velikost zkreslení na centrálním poledníku daného pásu
e2 • e´ ... je čtverec druhé excentricity, e ′ = 1 − e2 2
2
(4)
•
N ... příčný poloměr křivosti, N =
a 1 − e 2 sin 2 ϕ
(5)
•
T (substituce) ... T = tg 2ϕ (6)
•
C ... C = e′ 2 cos 2 ϕ (7)
•
A ... A = (λ − λ0 ) cos ϕ
... (zeměpisná délka je počítána v radiánech)
(8)
29
•
⎡⎛ e 2 3e 4 5e 6 ⎞ ⎛ 3e 2 3e 4 45e 6 ⎞ M = a ⎢⎜⎜1 − − − − K⎟⎟ϕ − ⎜⎜ + + + K⎟⎟ sin 2ϕ + 4 64 256 32 1024 ⎠ ⎝ 8 ⎠ ⎣⎝ ⎤ ⎛ 15e 4 45e 6 ⎞ ⎛ 35e 6 ⎞ + ⎜⎜ + + K⎟⎟ sin 4ϕ − ⎜⎜ + K⎟⎟ sin 6ϕ + K⎥ ⎝ 256 1024 ⎠ ⎝ 3072 ⎠ ⎦ (9)
M0 = M počítáno pro φ0, šířku protínající centrální poledník λ0 v počátku souřadnicové soustavy poledníkového pásu x,y. Je-li,
ϕ=±
π 2 potom můžeme vynechat všechny rovnice kromě rovnice
(9) , ze kterých vycházíme pro M a M0. Potom x = 0, y = k0(M - M0), k = k0. Rovnice (3) pro k může být rovněž psána jako funkce φ a x:
⎡ x2 ⎤ 2 2 k = k 0 ⎢1 + 1 + e′ cos ϕ ⎥ 2k 02 N 2 ⎦ ⎣
(
)
(10)
K získání souřadnic UTM je přivlastněný „ falešný východ“ přidán k ose x a „falešný sever “ k ose y po výpočtu z rovnic (1) a (2) . Pro inverzní výpočet potom platí: 4 ⎛ ⎛ ϕ1 ⎞ ⎞ ⎡ D 2 2 2 D ⎜ ⎟ + ϕ = ϕ1 − ⎜ N 1tg ⎜⎜ ⎟⎟ ⎟ ⋅ ⎢ − 5 + 3T1 + 10C1 − 4C1 − 9e ′ 24 ⎝ R1 ⎠ ⎠ ⎣ 2 ⎝
(
)
(
+ 61 + 90T1 + 298C1 + 45T12 − 252 e ′ 2 − 3C12
D ⎤ ) 720 ⎥ 6
⎦
(11)
λ = λ0 + [D − (1 + 2T1 + C1 )
D3 D5 ⎤ + 5 − 2C1 + 28T1 − 3C12 + 8e′ 2 + 24T12 ⎥ cos ϕ1 , 6 120 ⎦
(
)
(12)
která vychází z rovnice pro φ1:
⎛ 3e1
ϕ1 = μ + ⎜⎜
⎝ 2
−
⎞ ⎛ 21e12 55e14 ⎞ ⎛ 151e13 ⎞ 27e13 + K⎟⎟ sin 2 μ + ⎜⎜ − + K⎟⎟ sin 4 μ + ⎜⎜ + K⎟⎟ sin 6 μ + 32 32 ⎠ ⎝ 16 ⎠ ⎝ 96 ⎠
⎛ 1097e14 ⎞ + ⎜⎜ − K⎟⎟ sin 8μ + K ⎝ 512 ⎠ (13)
30
kde : •
e1 .. excentricita : e1 =
1 − (1 − e 2 ) 1 + (1 − e 2 )
(14) •
μ ... μ =
M ⎛ ⎞ 3e 4 5e 6 e − − − K⎟⎟ a⎜⎜1 − 4 64 256 ⎝ ⎠ 2
(15) •
M ... M = M 0 +
y , kde M0 se počítá z rovnice (9) pro dané φ0. Index 1 k0
znamená, že v uvedeném vztahu pro danou veličinu zaměníme φza φ1. •
D=
x N1k 0
(16)
•
R1 =
a (1 − e 2 ) 3 (1 − e 2 sin 2 ϕ1 ) 2
(17)
Hlásný systém MGRS je pouze kódovaným zápisem rovinných souřadnic zobrazení UTM.
31
6. Popis programu Program Convert.exe počítá převod souřadnic mezi UTM (N,E), elipsoidem WGS 84 (φ, λ neboli lat, lon) a systémem MGRS. Program pracuje ve 3 modulech:
•
převod z UTM do WGS 84 a MGRS vstup: N, E, zone, layer výstup: lat, lon, kód MGRS
•
převod z WGS 84 do UTM a MGRS vstup: lat, lon výstup: N, E, zone, layer, kód MGRS
•
převod z MGRS do UTM a WGS 84 vstup: kód MGRS výstup: N, E, lat, lon Výběr modulu probíhá v prvním dialogovém okně po spuštění
programu.
Obr.8 : Dialogové okno pro výběr modulu
ad 1) UTM → WGS 84, MGRS vstupní parametry:
•
N … northing, určuje, na které polokouli bod je, nabývá hodnot od 0 do 10 000 000m (na severní polokouli hodnota od rovníku směrem k pólu stoupá, na jižní polokouli je tomu obráceně)
32
•
E … easting, nabývá hodnot od 0 do 668 0000m (počátek souřadnic má hodnotu x = 500 000)
•
Zone … číslo poledníkového pásu, nabývá hodnot od 1 do 60 (rozsah pásu je 6°)
•
Layer … vrstva, je popsána písmeny od C do X , kde I, O se nepoužívají (značení začíná od 80° jižní šířky směrem na sever, rozsah vrstvy je 8°) Postup výpočtu : a) převod na zeměpisné souřadnice:
-
přepočet délky poledníkového oblouku z rovnice (16)
-
výpočet zeměpisných souřadnic dle postupu uvedeného v matematickém úvodu (5.3)
b) převod do MGRS :
-
první část kódu je číslo příslušného poledníkového pásu a označení vrstvy, což se přepíše ze zadaných hodnot
-
ve druhé části jsou dvě písmena označující příslušný 100km čtverec, první písmeno určuje sloupec, druhé řádek
-
poslední údaj uvádí polohu bodu v daném čtverci, počátek souřadnicové soustavy je v levém dolním rohu čtverce, údaj se shoduje se souřadnicí v UTM v řádech desítek tisíc, tedy je-li E = 604 819,71 a N = 1 690 908,32, potom koncový údaj kódu bude 0491890908
Obr.9: Převod z UTM
33
ad 2) WGS → UTM, MGRS vstupní parametry:
•
lat … latitude, zeměpisná šířka, je to úhel mezi rovinou normály v bodě na povrchu elipsoidu a rovinou rovníku, nabývá hodnot od –90°do 90°, rozsah UTM je od –84°do 80°
•
lon … longitude, zeměpisná délka, je to úhel, který svírá rovina místního poledníku, procházejícího určovaným bodem a rovina Greenwichského (nultého) poledníku, nabývá hodnot od –180° do 180° Postup výpočtu: a) převod do UTM:
-
výpočet zóny (poledníkového pásu)
-
výpočet délky poledníkového oblouku od rovníku k dané zeměpisné šířce
-
výpočet zeměpisných souřadnic uvedený v kapitole (5.3)
b) převod do MGRS je prováděn stejným způsobem jako v předchozím modulu
Obr.10: Převod WGS 84ÆUTM, MGRS
ad 3) MGRS → UTM, WGS 84 vstupní parametry:
•
kód MGRS … alfanumerický kód souřadnic UTM, 14-15 znaků (podle počtu číslic v označení poledníkového pásu)
34
postup výpočtu: a) převod do UTM:
-
číslo pásu a označení vrstvy označuje sférický čtyřúhelník (určený vrstvou a poledníkovým pásem), ve kterém se bod nachází
-
převod úhlů na délkový element (z definice metru jako desetimiliontiny zemského kvadrantu, tedy vzdálenosti od pólu k rovníku)
-
určení souřadnic UTM
b) převod do WGS 84 se provádí z určených souřadnic UTM stejným způsobem jako v modulu
Obr.10: Převod MGRSÆUTM, WGS 84
6.1 Microsoft Visual C++ 6.0 Tato verze vývojového prostředí Visual C++ je starší, přesto je možné ji , alespoň částečně, přizpůsobit aktuálnímu stavu vývojových prostředí. Můj výběr software pro tuto práci ovlivnily tři věci: 1) znalost (alespoň základní) jazyka C++ , 2) dostupnost vývojového prostředí (možnosti, které mám doma, v blízkém okolí), 3) fakt, že 90% uživatelů používá produkty společnosti
35
Microsoft a proto je využitelnost výsledného produktu maximální. Vybrané prostředí pro mě představovalo něco nového a naprosto odlišného od dosavadních zkušeností, jelikož programování v tomto prostředí neznamená pouze sestavování algoritmů, ale předpokládá se i určitá znalost operačního systému Windows.
36
Závěr Tato bakalářská práce se zabývá výpočty převodů mezi zeměpisnými souřadnicemi z elipsoidu WGS 84 [φ, λ], rovinnými souřadnicemi zobrazení UTM [N, E] a alfanumerickým systémem MGRS. Převod se týká pouze rozsahu zobrazení UTM, tedy od 80° j.š. do 84°s.š.. Výsledkem celé práce je program pro automatický souřadnicový převod. Program byl vytvářen ve vývojovém prostředí Microsoft Visual C++ 6.0. Vznik funkčního programu vyžadoval rozšíření stávajících znalostí z oblasti programování v jazyce C++, studium nového vývojového prostředí Visual C++ a doplnění znalostí z matematické kartografie. Výběr software byl již naznačen v textové části práce. Představoval pro mě něco zcela nového, využívá ale jazyk C++, který byl součástí výuky mého dosavadního studia na FSv. Při výběru vývojového prostředí jsem se rozhodovala pouze mezi těmi, které jsou pro mne dostupné na domácím PC. Tedy mezi Visual Basic a Visual C++, freeware verze jiných prostředí, dostupných na internetu, mi nebyly před tím známy. Ačkoli C++ patří mezi programovací jazyky se složitější syntaxí, z uvedených dvou voleb mi byl samozřejmě bližší. Práce na finální podobě programu se neobešla bez potíží a nebyla pro mě jako pro laika jednoduchá, přesto věřím, že výsledná podoba programu odpovídá požadavkům zadání. Program převádí souřadnice obousměrně mezi uvedenými systémy. Výběr směru převodu uživatel volí v prvním dialogovém okně.
37
Seznam použité literatury: [1] Veverka B. : Kartografické standardy NATO. Geodetický a kartografický obzor, 45/87, 1999, č.7-8, str.140-147. [2] Kubátová R. : Systém JTSK a WGS 84, jejich charakteristika a vzájemná transformace, Bakalářská práce, ZČU v Plzni [3] Jurkina M.A.,Pick M. : Numerické výpočty ve světovém geodetickém referenčním systému 1984 (WGS 84), Vojenský geografický obzor, 2006, č.1, příloha 2 [4] Geocaching Encyclopaedia, článek WGS 84 , dostupný na http://wiki.geocaching.cz/wiki/WGS-84 [5] Hánek.: Zobrazení UTM, text pro projekt CTU 0513011 (2005), dostupné na http://klobouk.fsv.cvut.cz/~hanek/K154/PDF/UTM.pdf [6] Divišová L.: UTM, semestrální práce z tématické kartografie, ZČU Plzeň, 2000 [7] Smýkalová R.: Mercatorův přínos pro matematickou kartografii, 8.vědecká konference doktorandů a mladých vědeckých pracovníků, 2007, FPV UKF Nitra, dostupné na http://citadel.ukf.sk/konferencia/papers/PDF_Matematika/Smykalova.pdf [8] Snyder, J.P.: Map Projections :A Working Manual, 1987, Geological Survey (U.S.) [9] Buchar P. : Matematická kartografie, skripta, 2005, ČVUT [10] Talhofer V.: Základy matematické kartografie, studijní texty, Univerzita obrany, Fakulta vojenských technologií, katedra vojenské geografie a meteorologie 2007, Brno
[11] Wikipedia,The Free Encyclopedia, article:Universal Transverse Mercator, dostupné na http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system
Seznam příloh: Příloha I: Zdrojový kód programu Convert.exe Příloha II: Zobrazení UTM, značení sítě
38
Příloha I: Zdrojový kód programu Convert.exe Data.cpp #include "stdafx.h" #include "resource.h" #include "param.h" // Global Variables: HWND hWnd; HINSTANCE hInst; // current instance TCHAR szTitle[MAX_LOADSTRING]; // The title bar text TCHAR szWindowClass[MAX_LOADSTRING]; title bar text
// The
// OPENFILE structure and it's pointer LPOPENFILENAME lpofn; OPENFILENAME ofn; // File Name buffer char FileName[MAX_LOADSTRING]; // File structures FILE *in, *out; BOOL fn;
file 2.cpp / file2.cpp : Defines the entry point for the application. // #include "stdafx.h" #include "resource.h" #include "param.h" #include "proto.h" int APIENTRY WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow) { // TODO: Place code here. MSG msg; HACCEL hAccelTable; // Initialize global strings LoadString(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING); LoadString(hInstance, IDC_FILE2, szWindowClass, MAX_LOADSTRING); MyRegisterClass(hInstance); // Perform application initialization: if (!InitInstance (hInstance, nCmdShow)) { return FALSE;
39
} hAccelTable = LoadAccelerators(hInstance, (LPCTSTR)IDC_FILE2); // Main message loop: while (GetMessage(&msg, NULL, 0, 0)) { if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg)) { TranslateMessage(&msg); DispatchMessage(&msg); } } return msg.wParam; }
// // FUNCTION: MyRegisterClass() // // PURPOSE: Registers the window class. // // COMMENTS: // // This function and its usage is only necessary if you want this code // to be compatible with Win32 systems prior to the 'RegisterClassEx' // function that was added to Windows 95. It is important to call this function // so that the application will get 'well formed' small icons associated // with it. // ATOM MyRegisterClass(HINSTANCE hInstance) { WNDCLASSEX wcex; wcex.cbSize = sizeof(WNDCLASSEX); wcex.style wcex.lpfnWndProc wcex.cbClsExtra wcex.cbWndExtra wcex.hInstance wcex.hIcon wcex.hCursor wcex.hbrBackground wcex.lpszMenuName wcex.lpszClassName wcex.hIconSm
= CS_HREDRAW | CS_VREDRAW; = (WNDPROC)WndProc; = 0; = 0; = hInstance; = LoadIcon(hInstance, (LPCTSTR)IDI_FILE2); = LoadCursor(NULL, IDC_ARROW); = (HBRUSH)(COLOR_WINDOW+1); = (LPCSTR)IDC_FILE2; = szWindowClass; = LoadIcon(wcex.hInstance, (LPCTSTR)IDI_SMALL);
return RegisterClassEx(&wcex); } // // FUNCTION: InitInstance(HANDLE, int) // // PURPOSE: Saves instance handle and creates main window // // COMMENTS: // // In this function, we save the instance handle in a global variable and // create and display the main program window.
40
// BOOL InitInstance(HINSTANCE hInstance, int nCmdShow) { hInst = hInstance; // Store instance handle in our global variable hWnd = CreateWindow(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW, CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, NULL, NULL, hInstance, NULL); if (!hWnd) { return FALSE; } // ShowWindow(hWnd, nCmdShow); ShowWindow(hWnd, SW_HIDE); // ShowWindow(hWnd, SW_SHOWNORMAL); UpdateWindow(hWnd); return TRUE; } // // FUNCTION: WndProc(HWND, unsigned, WORD, LONG) // // PURPOSE: Processes messages for the main window. // // WM_COMMAND - process the application menu // WM_PAINT - Paint the main window // WM_DESTROY - post a quit message and return // // LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam) { int wmId, wmEvent; PAINTSTRUCT ps; HDC hdc; TCHAR szHello[MAX_LOADSTRING]; LoadString(hInst, IDS_HELLO, szHello, MAX_LOADSTRING);
switch (message) { case WM_COMMAND: wmId = LOWORD(wParam); wmEvent = HIWORD(wParam); // Parse the menu selections: switch (wmId) { case IDD_DIALOG_MAIN: DialogBox(hInst, (LPCTSTR)IDD_DIALOG_MAIN, hWnd, (DLGPROC)MyMainDialog); break; case IDD_WGS2UTM: DialogBox(hInst, (LPCTSTR)IDD_WGS2UTM, hWnd, (DLGPROC)wgs2utm); break; case IDD_UTM2WGS: DialogBox(hInst, (LPCTSTR)IDD_UTM2WGS, hWnd, (DLGPROC)utm2wgs);
41
break; case IDD_MGRS: DialogBox(hInst, (LPCTSTR)IDD_MGRS, hWnd, (DLGPROC)mgrs); break; case IDD_BATCH: DialogBox(hInst, (LPCTSTR)IDD_BATCH, hWnd, (DLGPROC)BatchDlg); break; case IDM_ABOUT: DialogBox(hInst, (LPCTSTR)IDD_ABOUTBOX, hWnd, (DLGPROC)About); break; case IDM_EXIT: DestroyWindow(hWnd); break; default: return DefWindowProc(hWnd, message, wParam, lParam); } break; case WM_CREATE: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDD_DIALOG_MAIN, NULL); break; case WM_PAINT: hdc = BeginPaint(hWnd, &ps); // TODO: Add any drawing code here... // RECT rt; // GetClientRect(hWnd, &rt); // DrawText(hdc, szHello, strlen(szHello), &rt, DT_CENTER); EndPaint(hWnd, &ps); // SendMessage(hWnd, WM_COMMAND, (WPARAM)IDD_DIALOG_MAIN, NULL); break; case WM_DESTROY: PostQuitMessage(0); break; default: return DefWindowProc(hWnd, message, wParam, lParam); } return 0; } // Mesage handler for about box. LRESULT CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam) { switch (message) { case WM_INITDIALOG: return TRUE; case WM_COMMAND: if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL) { EndDialog(hDlg, LOWORD(wParam)); return TRUE; } break; } return FALSE;
42
}
mgrs.cpp #include "stdafx.h" #include "resource.h" #include <stdio.h> #include <string.h> #include <process.h> #include <math.h> #include
#include <winnls.h> #include #include "proto.h" #define MAX_LOADSTRING 100 extern HWND
hWnd;
int getLatZoneDegree(char l) { char letters[] = { 'A', 'C', 'D', 'E', 'F', 'G', 'H', 'J', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Z' }; int degrees[] = { -90, -84, -72, -64, -56, -48, -40, -32, -24, -16, -8, 0, 8, 16, 24, 32, 40, 48, 56, 64, 72, 84 }; for (int i = 0; i < 22; i++) { if (letters[i] == l) { return degrees[i]; } } return -100; } int getDigraph1Index(char letter) { char digraph1Array[] = { 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'J', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z' }; for (int i = 0; i < 24; i++) { if (digraph1Array[i]==letter) { return i+1; } } return -1; } int getDigraph2Index(char letter) { char digraph2Array[] = { 'V', 'A', 'B', 'C', 'D', 'E', 'F', 'G',
43
'H', 'J', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'U', 'V' }; for (int i = 0; i < 21; i++) { if (digraph2Array[i]==letter) { return i; } } return -1; } LRESULT CALLBACK mgrs(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam) { char degb[20], minb[20], secb[20], degl[20], minl[20], secl[20], mgrs[50]; char text[10]; long zone, eutm, nutm; char latZone, layer, row, col; double a1, a2; long latZoneDegree; double digraph1Index, digraph2Index; double startindexEquator; double zoneCM; double a6[] = { 16, 0, 8 }; long a5; double a7; double a3; double easting, northing; char digraph1, digraph2; char str[20]; bool minus; switch (message) { case WM_INITDIALOG: return TRUE; break; case WM_COMMAND: switch(wParam) { case IDOK: // OK button pressed // algorithmic processing of input data GetDlgItemText(hDlg, IDC_ZONEM, text, 5); zone = atoi(text); GetDlgItemText(hDlg, IDC_LAYERM, text, 5); latZone = (char)text[0]; GetDlgItemText(hDlg, IDC_COLM, text, 5); digraph1 = (char)text[0]; GetDlgItemText(hDlg, IDC_ROWM, text, 5); digraph2 = (char)text[0]; GetDlgItemText(hDlg, IDC_EM, text, 8); easting = atof(text); GetDlgItemText(hDlg, IDC_NM, text, 8); northing = atof(text); latZoneDegree = getLatZoneDegree(latZone);
44
if (latZoneDegree < 0) minus = true; else minus = false; if (minus) { switch(latZoneDegree) { case -84: latZoneDegree = 0; break; case -72: latZoneDegree = 8; break; case -64: latZoneDegree = 16; break; case -56: latZoneDegree = 24; break; case -48: latZoneDegree = 32; break; case -40: latZoneDegree = 40; break; case -32: latZoneDegree = 48; break; case -24: latZoneDegree = 56; break; case -16: latZoneDegree = 64; break; case -8: latZoneDegree = 72; break; case 0: latZoneDegree = 84; break; } } a1 = (double)latZoneDegree * 40000000 / 360.0; if (minus) a2 = 2000000 * ceil(a1 / 2000000.0); else a2 = 2000000 * floor(a1 / 2000000.0); digraph2Index = getDigraph2Index(digraph2); startindexEquator = 1; if (( 1 + zone % 2) == 1) { startindexEquator = 6; } a3 = a2 + (digraph2Index - startindexEquator) * 100000; if (a3 <= 0) { a3 = 10000000 + a3; } northing = a3 + northing; if (northing == 10000000) northing = 0; zoneCM = -183 + 6 * zone; digraph1Index = getDigraph1Index(digraph1); a5 = 1 + zone % 3; a7 = 100000 * (digraph1Index - a6[a5 - 1]); easting = easting + a7; //LatLonToUtm (Bwgs84, Lwgs84, (LPSTR)&north, (LPSTR)&east, (LPSTR)&zone, (LPSTR)&layer, (LPSTR)&mgrs); sprintf(str, "%15.1f", easting); SetDlgItemText(hDlg, IDC_LATM, str); sprintf(str, "%15.1f", northing); SetDlgItemText(hDlg, IDC_LONM, str); UtmToLatLon (zone, latZone, northing, easting, (LPSTR)°b, (LPSTR)&minb, (LPSTR)&secb, (LPSTR)°l, (LPSTR)&minl, (LPSTR)&secl, (LPSTR)&mgrs); SetDlgItemText(hDlg, IDC_STUPNE_BM, (LPSTR)°b); SetDlgItemText(hDlg, IDC_MINUTY_BM, (LPSTR)&minb); SetDlgItemText(hDlg, IDC_VTERINY_BM, (LPSTR)&secb); SetDlgItemText(hDlg, IDC_STUPNE_LM, (LPSTR)°l); SetDlgItemText(hDlg, IDC_MINUTY_LM, (LPSTR)&minl); SetDlgItemText(hDlg, IDC_VTERINY_LM, (LPSTR)&secl);
45
break; case IDCANCEL: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDD_DIALOG_MAIN, NULL); EndDialog(hDlg, LOWORD(wParam)); break; case IDABOUT: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDM_ABOUT, NULL); break; } break; } return FALSE; }
MyMainDialog.cpp #include "stdafx.h" #include "resource.h" #include <stdio.h> #include <string.h> #include <process.h> #include #include <winnls.h> #include #include "proto.h" #define MAX_LOADSTRING 100 extern HWND hWnd; LRESULT CALLBACK MyMainDialog(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam) { switch (message) { case WM_INITDIALOG: return TRUE; case WM_COMMAND: switch(wParam) { case IDCANCEL: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDM_EXIT, NULL); EndDialog(hDlg, LOWORD(wParam)); break; case IDABOUT: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDM_ABOUT, NULL); break; case ID_WGS_TO_UTM: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDD_WGS2UTM, NULL); EndDialog(hDlg, LOWORD(wParam)); break; case ID_UTM_TO_WGS: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDD_UTM2WGS, NULL); EndDialog(hDlg, LOWORD(wParam)); break; case ID_MGRS: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDD_MGRS, NULL); EndDialog(hDlg, LOWORD(wParam));
46
break; case IDBATCH: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDD_BATCH, NULL); EndDialog(hDlg, LOWORD(wParam)); break; } break; } return FALSE; }
StdAfx.cpp // stdafx.cpp : source file that includes just the standard includes // file2.pch will be the pre-compiled header // stdafx.obj will contain the pre-compiled type information #include "stdafx.h" // TODO: reference any additional headers you need in STDAFX.H and not in this file
utm2wgs.cpp #include "stdafx.h" #include "resource.h" #include <stdio.h> #include <string.h> #include <process.h> #include #include <winnls.h> #include #include "proto.h"
#define MAX_LOADSTRING 100 extern HWND
hWnd;
double northing; LRESULT CALLBACK utm2wgs(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam) { char text[20]; int Zone; char layer; double north, east; char degb[20], minb[20], secb[20], degl[20], minl[20], secl[20], mgrs[50]; switch (message)
47
{ case WM_INITDIALOG: return TRUE; case WM_COMMAND: switch(wParam) { case IDCNS: GetDlgItemText(hDlg, IDNORTH, text, 15); if (text[0] == 'N') { SetDlgItemText(hDlg, IDNORTH, (LPSTR)&"S"); northing = -10000000; } else { SetDlgItemText(hDlg, IDNORTH, (LPSTR)&"N"); northing = 0; } break; case IDOK:
// OK button pressed // algorithmic processing of input data GetDlgItemText(hDlg, IDC_NORTHING, text, 15); north = atof(text); GetDlgItemText(hDlg, IDC_EASTING, text, 15); east = atof(text); GetDlgItemText(hDlg, IDC_ZONE, text, 15); Zone = atoi(text); GetDlgItemText(hDlg, IDI_LAYER, text, 15); layer = (char)text[0]; //string s(text);// string sub = s.substr(21); north = north + northing; // north / south decission
UtmToLatLon (Zone, layer, north, east, (LPSTR)°b, (LPSTR)&minb, (LPSTR)&secb, (LPSTR)°l, (LPSTR)&minl, (LPSTR)&secl, (LPSTR)&mgrs); SetDlgItemText(hDlg, IDC_STUPNE_B, (LPSTR)°b); SetDlgItemText(hDlg, IDC_MINUTY_B, (LPSTR)&minb); SetDlgItemText(hDlg, IDC_VTERINY_B, (LPSTR)&secb); SetDlgItemText(hDlg, IDC_STUPNE_L, (LPSTR)°l); SetDlgItemText(hDlg, IDC_MINUTY_L, (LPSTR)&minl); SetDlgItemText(hDlg, IDC_VTERINY_L, (LPSTR)&secl); SetDlgItemText(hDlg, IDC_MGRS, (LPSTR)&mgrs); break; Case IDCANCEL: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDD_DIALOG_MAIN, NULL); EndDialog(hDlg, LOWORD(wParam)); break; case IDABOUT: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDM_ABOUT, NULL); break; } break; } return FALSE; }
wgs2utm.cpp #include "stdafx.h" #include "resource.h"
48
#include <stdio.h> #include <string.h> #include <process.h> #include #include <winnls.h> #include #include "proto.h"
#define MAX_LOADSTRING 100 extern HWND
hWnd;
LRESULT CALLBACK wgs2utm(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam) { double width_deg, width_min; double width_sec, width; double length_deg, length_min; double length_sec, length; double K1=0.863499977, K2=0.504348876, K3=42.52539444, deltaV, pi=3.141592654; double sinw, sinl; char text[10], mgrs[50]; int i_width, i; double Bwgs84, Lwgs84; double degb, minb, secb, degl, minl, secl; char east[30], north[30], zone[30], layer[30]; switch (message) { case WM_INITDIALOG: return TRUE; case WM_COMMAND: switch(wParam) { case IDOK: // OK button pressed // algorithmic processing of input data GetDlgItemText(hDlg, IDC_STUPNE_B, text, 8); degb = atof(text); GetDlgItemText(hDlg, IDC_MINUTY_B, text, 3); minb = atof(text); GetDlgItemText(hDlg, IDC_VTERINY_B, text, 7); secb = atof(text); GetDlgItemText(hDlg, IDC_STUPNE_L, text, 8); degl = atof(text); GetDlgItemText(hDlg, IDC_MINUTY_L, text, 3); minl = atof(text); GetDlgItemText(hDlg, IDC_VTERINY_L, text, 7); secl = atof(text); if (degb < 0) {
49
minb = -minb; secb = -secb; } if (degl < 0) { minl = -minl; secl = -secl; } Bwgs84 = (degb + minb/60 + secb/3600); Lwgs84 = (degl + minl/60 + secl/3600); LatLonToUtm (Bwgs84, Lwgs84, (LPSTR)&north, (LPSTR)&east, (LPSTR)&zone, (LPSTR)&layer, (LPSTR)&mgrs); SetDlgItemText(hDlg, IDC_MGRS, (LPSTR)&mgrs); SetDlgItemText(hDlg, IDC_NORTHING, (LPSTR)&north); SetDlgItemText(hDlg, IDC_EASTING, (LPSTR)&east); SetDlgItemText(hDlg, IDC_ZONE, (LPSTR)&zone); SetDlgItemText(hDlg, IDC_LAYER, (LPSTR)&layer); i = 0; break; case IDCANCEL: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDD_DIALOG_MAIN, NULL); EndDialog(hDlg, LOWORD(wParam)); break; case IDABOUT: PostMessage(hWnd, WM_COMMAND, (WPARAM)IDM_ABOUT, NULL); break; } break; } return FALSE; }
wgs2utm_proc.cpp #include "stdafx.h" #include "resource.h" #include <stdio.h> #include <string.h> #include <process.h> #include #include <winnls.h> #include #include "proto.h" #define MAX_LOADSTRING 100 // Some constants used by these functions. static const double fe = 500000.0; static const double ok = 0.9996; static const double pi = 3.141592654; static const double a = 6378137.0; static const double f = 1 / 298.257223563; static const double deg2rad = pi/180;
50
static const double rad2deg = 180/pi; // An array containing each vertical UTM zone. static char cArray[] = "CDEFGHJKLMNPQRSTUVWX"; double CalculateESquared (double a, double b) { return ((a * a) - (b * b)) / (a * a); } double CalculateE2Squared (double a, double b) { return ((a * a) - (b * b)) / (b * b); } double denom (double es, double sphi) { double sinSphi = sin (sphi); return sqrt (1.0 - es * (sinSphi * sinSphi)); } double sphsr (double a, double es, double sphi) { double dn = denom (es, sphi); return a * (1.0 - es) / (dn * dn * dn); } double sphsn (double a, double es, double sphi) { double sinSphi = sin (sphi); return a / sqrt (1.0 - es * (sinSphi * sinSphi)); } double sphtmd (double ap, double bp, double cp, double dp, double ep, double sphi) { return (ap * sphi) - (bp * sin (2.0 * sphi)) + (cp * sin (4.0 * sphi)) - (dp * sin (6.0 * sphi)) + (ep * sin (8.0 * sphi)); } void LatLonToUtm (double lat, double lon, LPSTR Nutm, LPSTR Eutm, LPSTR _zone, LPSTR _layer, LPSTR mgrs) { double recf; double b; double eSquared; double e2Squared; double tn; double ap; double bp; double cp; double dp; double ep; double olam; double dlam; double s; double c; double t; double eta; double sn; double tmd;
51
double t1, t2, t3, t6, t7; double nfn; double zone, northing, easting; char layer, error[]={"Error"}; zone = (long)((lon + 180) / 6 + 1); if (zone == 61) zone = 1; if (lon > 0 && lon < 9 && lat > 72 && lat < 84) zone = 31; if (lon > 9 && lon < 21 && lat > 72 && lat < 84) zone = 33; if (lon > 21 && lon < 33 && lat > 72 && lat < 84) zone = 35; if (lon > 33 && lon < 42 && lat > 72 && lat < 84) zone = 37; if (lon > 3 && lon < 9 && lat > 56 && lat < 64) zone = 32; if (lat < 84.0 && lat >= 72.0) { // Special case: zone X is 12 degrees from north to south, not 8. layer = cArray[19]; } else { layer = cArray[(int)((lat + 80.0) / 8.0)]; } if (lat >= 84.0 || lat < -80.0) { // Invalid coordinate; the vertical zone is set to the invalid // character. layer = '*'; } double latRad = lat * deg2rad; double lonRad = lon * deg2rad; recf = 1.0 / f; b = a * (recf - 1.0) / recf; eSquared = CalculateESquared (a, b); e2Squared = CalculateE2Squared (a, b); tn = (a - b) / (a + b); ap = a * (1.0 - tn + 5.0 * ((tn * tn) - (tn * tn * tn)) / 4.0 + 81.0 * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 64.0); bp = 3.0 * a * (tn - (tn * tn) + 7.0 * ((tn * tn * tn) - (tn * tn * tn * tn)) / 8.0 + 55.0 * (tn * tn * tn * tn * tn) / 64.0) / 2.0; cp = 15.0 * a * ((tn * tn) - (tn * tn * tn) + 3.0 * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 4.0) / 16.0; dp = 35.0 * a * ((tn * tn * tn) - (tn * tn * tn * tn) + 11.0 * (tn * tn * tn * tn * tn) / 16.0) / 48.0; ep = 315.0 * a * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 512.0; olam = (zone * 6 - 183) * deg2rad; dlam = lonRad - olam; s = sin (latRad); c = cos (latRad); t = s / c; eta = e2Squared * (c * c); sn = sphsn (a, eSquared, latRad); tmd = sphtmd (ap, bp, cp, dp, ep, latRad); t1 = tmd * ok; t2 = sn * s * c * ok / 2.0; t3 = sn * s * (c * c * c) * ok * (5.0 - (t * t) + 9.0 * eta + 4.0 * (eta * eta)) / 24.0; if (latRad < 0.0) nfn = 10000000.0; else nfn = 0; northing = nfn + t1 + (dlam * dlam) * t2 + (dlam * dlam * dlam * dlam) * t3 + (dlam * dlam * dlam * dlam * dlam * dlam) + 0.5; t6 = sn * c * ok; t7 = sn * (c * c * c) * (1.0 - (t * t) + eta) / 6.0; easting = fe + dlam * t6 + (dlam * dlam * dlam) * t7 + 0.5;
52
if (northing >= 9999999.0) northing = 9999999.0; sprintf(Nutm, "%8.1f", northing); sprintf(Eutm, "%8.1f", easting); sprintf(_zone, "%d", (int)zone); sprintf(_layer, "%c", layer); int i = 0; if ((northing < 10000000) && (easting < 1000000)) { mgrs_proc(northing, easting, zone, layer, mgrs); } else { while(error[i]!=0) *mgrs++ = error[i++]; *mgrs = 0; } } void UtmToLatLon (int zone, char layer, double northing, double easting, LPSTR degb, LPSTR minb, LPSTR secb, LPSTR degl, LPSTR minl, LPSTR secl, LPSTR mgrs) { double recf; double b; double eSquared; double e2Squared; double tn; double ap; double bp; double cp; double dp; double ep; double nfn; double tmd; double sr; double sn; double ftphi; double s; double c; double t; double eta; double de; double dlam; double olam; double lat, lon; long deg, min, sec; recf = 1.0 / f; b = a * (recf - 1) / recf; eSquared = CalculateESquared (a, b); e2Squared = CalculateE2Squared (a, b); tn = (a - b) / (a + b); ap = a * (1.0 - tn + 5.0 * ((tn * tn) - (tn * tn * tn)) / 4.0 + 81.0 * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 64.0); bp = 3.0 * a * (tn - (tn * tn) + 7.0 * ((tn * tn * tn) - (tn * tn * tn * tn)) / 8.0 + 55.0 * (tn * tn * tn * tn * tn) / 64.0) / 2.0; cp = 15.0 * a * ((tn * tn) - (tn * tn * tn) + 3.0 * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 4.0) / 16.0; dp = 35.0 * a * ((tn * tn * tn) - (tn * tn * tn * tn) + 11.0 * (tn * tn * tn * tn * tn) / 16.0) / 48.0; ep = 315.0 * a * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 512.0;
53
if ((layer <= 'M' && layer >= 'C') || (layer <= 'm' && layer >= 'c')) { nfn = 10000000.0; } else { nfn = 0; } tmd = (northing - nfn) / ok; sr = sphsr (a, eSquared, 0.0); ftphi = tmd / sr; double t10, t11, t14, t15; for (int i = 0; i < 5; i++) { t10 = sphtmd (ap, bp, cp, dp, ep, ftphi); sr = sphsr (a, eSquared, ftphi); ftphi = ftphi + (tmd - t10) / sr; } sr = sphsr (a, eSquared, ftphi); sn = sphsn (a, eSquared, ftphi); s = sin (ftphi); c = cos (ftphi); t = s / c; eta = e2Squared * (c * c); de = easting - fe; t10 = t / (2.0 * sr * sn * (ok * ok)); t11 = t * (5.0 + 3.0 * (t * t) + eta - 4.0 * (eta * eta) - 9.0 * (t * t) * eta) / (24.0 * sr * (sn * sn * sn) * (ok * ok * ok * ok)); lat = ftphi - (de * de) * t10 + (de * de * de * de) * t11; t14 = 1.0 / (sn * c * ok); t15 = (1.0 + 2.0 * (t * t) + eta) / (6 * (sn * sn * sn) * c * (ok * ok * ok)); dlam = de * t14 - (de * de * de) * t15; olam = (zone * 6 - 183.0) * deg2rad; lon = olam + dlam; deg = (long)(lon*rad2deg); min = (long)((lon*rad2deg - deg)*60); sec = (long)( (lon*rad2deg*3600) - ((deg*3600)+(min*60)) ); if(min < 0) min = - min; if(sec < 0) sec = - sec; sprintf(degl, "%d", deg); sprintf(minl, "%d", min); sprintf(secl, "%d", sec); deg = (long)(lat*rad2deg); min = (long)((lat*rad2deg - deg)*60); sec = (long) ( (lat*rad2deg*3600) - ((deg*3600)+(min*60)) ); if(min < 0) min = - min; if(sec < 0) sec = - sec; sprintf(degb, "%d", deg); sprintf(minb, "%d", min); sprintf(secb, "%d", sec); mgrs_proc(northing, easting, zone, layer, mgrs); } void utmwgs_proc(int Zone, double north, double east, LPSTR degb, LPSTR minb, LPSTR secb, LPSTR degl, LPSTR minl, LPSTR secl, LPSTR mgrs) { double a, b, E2, EE, Eutmred, L0, BR, DB, T, N, B2, B3, L1, L2, L3, A;
54
double Ro; double Eutm, Nutm, B, L; long deg, min, layer; double sec; char mgr[50], clayer; Eutm = east; Nutm = north; if (Nutm < 0) Nutm = Nutm + 10000000; Eutmred = (Eutm - 500000) / 0.9996;
//positive value // redukovana delka
char letter(int index) { char col; switch(index) { case case case case case case case case case case case case case case case case
1: col = 'A'; break; 2: col = 'B'; break; 3: col = 'C'; break; 4: col = 'D'; break; 5: col = 'E'; break; 6: col = 'F'; break; 7: col = 'G'; break; 8: col = 'H'; break; 9: col = 'J'; break; 10: col = 'K'; break; 11: col = 'L'; break; 12: col = 'M'; break; 13: col = 'N'; break; 14: col = 'P'; break; 15: col = 'Q'; break; 16:
55
case case case case case case case case
col = 'R'; break; 17: col = 'S'; break; 18: col = 'T'; break; 19: col = 'U'; break; 20: col = 'V'; break; 21: col = 'W'; break; 22: col = 'X'; break; 23: col = 'Y'; break; 24: col = 'Z'; break;
default: col = 'Z'; break; } return col; } void {
mgrs_proc(double Nutm, double Eutm, long Zone, char clayer, LPSTR mgr)
char col, line, text[20]; int layer[4] = {8,6,4,2}; int windex, bg, j; double mgrsl, mgrsw, ref; LPSTR mg2; mg2 = mgr; mgrsl = Eutm; windex = 0; bg = ((Zone % (24/layer[windex]) - 1) * layer[windex]) + 1; bg = (((Zone-1) % 3) * 8) + 1; col = letter(bg); if (mgrsl <500000) { j = bg + layer[windex]/2; while(mgrsl < 500000) { mgrsl = mgrsl + 100000; j--; } mgrsl = mgrsl - 500000; }
56
else { j = bg + layer[windex]/2; ref = 500000; while(mgrsl >= ref) { ref = ref + 100000; j++; } j--; ref = ref - 100000; mgrsl = mgrsl - ref; } col = letter(j%24); // altitude j = 1; ref = 0; mgrsw = Nutm; if (mgrsw >= 0) { while(mgrsw >= ref) { ref = ref + 100000; j++; } ref = ref - 100000; j--; mgrsw = mgrsw - ref; if ((Zone % 2) == 0) { j = j + 5; } j = j % 20; line = letter(j); } else { mgrsw = -mgrsw; while(mgrsw >= ref) { ref = ref + 100000; j--; if(j == 0 ) j = 20; } ref = ref - 100000; mgrsw = mgrsw - ref; mgrsw = 100000- mgrsw; if ((Zone % 2) == 0) { j = j + 5; } j = j % 20; if(j == 0 ) j = 20; line = letter(j); clayer--; if (clayer=='O') clayer = 'N'; switch(clayer) { case 'X':
57
clayer = 'M'; break; case 'W': clayer = 'L'; break; case 'V': clayer = 'K'; break; case 'U': clayer = 'J'; break; case 'T': clayer = 'H'; break; case 'S': clayer = 'G'; break; case 'R': clayer = 'F'; break; case 'Q': clayer = 'E'; break; case 'P': clayer = 'D'; break; case 'N': clayer = 'C'; break; } } sprintf(text, "%d %c %c %c %05d %05d",Zone, clayer, col, line, (long)mgrsl, (long)mgrsw); j = 0; while(text[j] != 0) *mg2++ = text[j++]; *mg2 = 0; return; }
58
Příloha II : Zobrazení UTM, značení sítě [11]
59
60