Ke zopakování: úhlové souřadnice na povrchu koule převedeme na 3D kartézské v prostoru. Pak k nim přičteme náhodný posuv s 3D normálním rozdělením. To se snadno vyrobí ze 3 nezávislých 1D normálních rozložení. Pravoúhlé souřadnice převedeme zpět do kulových. Souřadnici vzdálenosti od středu Země nepoužijeme, čímž se místo promítne kolmo na povrch Země.
2D normální rozdělení. Směrodatná odchylka je 1 km, stupnice jsou v km. 3D normální rozdělení vypadá na obrazovce zcela identicky. Je to neostrý kulovitý 3D obláček, jehož průmět do 2D je 2D normální rozdělení. Graf jsem vykreslil R příkazem plot(rnorm(4000),rnorm(4000),pch=20)
Nejdříve jsem hledal funkce na převod mezi kartézskými a kulovými souřadnicemi v GNU R. V CRAN balíku sphereplot je funkce sph2car. Nastává ale typické nabobtnalé tzv. dependency hell (komiks): CRAN balík sphereplot při instalaci vyžaduje CRAN balík rgl, což je interaktivní vizualizace 3D grafů pomocí OpenGL, a ten vyžaduje dalších 119 CRAN balíků, které zabírají na disku 339 MB. To už si ty funkce radši naprogramuju sám a celý program bude zabírat 31 řádků a 567 bajtů.
- voszm, vzdálenost od středu Země v metrech. Použijeme průměrný poloměr Země 6 371 000 m.
- zss, zeměpisná šířka ve stupních. 0 je rovník, +90 je severní pól, –90 je jižní pól.
- zds, zeměpisná délka ve stupních. 0 je Greenwich, +90 je Sumatra, –90 jsou Galapágy, –180 a 180 jsou identické, Howlandův Ostrov
Kartézské souřadnice
- xm, od středu Země v metrech, má zeměpisnou délku zds=0/+/-180, + směřuje zhruba k Africe, – zhruba k Howlandově ostrovu v Tichém oceánu, kde v roce 1937 zmizela slavná letkyně Amelia Earhartová.
- ym, od středu Země v metrech, má zeměpisnou délku zds=+/-90, + směřuje zhruba k Sumatře, – zhruba ke Galapágám
- zm, od středu Země v metrech, má zeměpisnou šířku zss=+/-90, + směřuje k Severnímu pólu, – k Jižnímu pólu
Osa Z vede od Jižního pólu (šířka –90, délka libovolná)…
… k Severnímu pólu (šířka 90, délka libovolná).
Konstanty a matematické funkce:
- pi, 3,1415926535897932384626433… Nepoužijeme vestavěnou konstantu GNU R pi, ale číslo přímo. Vestavěná konstanta-nekonstanta se dá totiž přepisovat a když si s GNU R bude někdo hrát, konstantu přepíše nepřesnou hodnotou a zůstane mu v prostředí uložena, výpočet mu poté i po potenciálně dlouhé době bude probíhat špatně, aniž by si toho všiml. float (single) má mantisu 24 bitů odpovídající 1,19 metrům na polovičním obvodu Země. double má mantisu 53 bitů odpovídající 2,2 nanometrům. GNU R v single přesnosti vůbec počítat neumí. GNU R příkaz help(double):
_R has no single precision data type. All real numbers are stored
in double precision format_.
- sin, cos berou argument v radiánech. 180 stupňů je π radiánů.
- arcsin je inverzní k sin. Vrací argument v radiánech.
- snr, stupně na radiány. Vynásobíme-li stupně konstantou snr, vyjde stejný úhel v radiánech.
- rns, radiány na stupně. Vynásobíme-li radiány konstantou rns, vyjde stejný úhel ve stupních.
xm a ym jsem původně převedl na zeměpisnou délku pomocí funkce arctan. Pro půlku kruhu ale dává špatné hodnoty. Proto použijeme funkci GNU R atan2:
Trigonometric Functions [...] Usage: [...] atan2(y, x) [...] The arc-tangent of two arguments ‘atan2(y, x)’ returns the angle between the x-axis and the vector from the origin to (x, y), i.e., for positive arguments ‘atan2(y, x) == atan(y/x)’. Angles are in radians, not degrees[...]
Jak bude funkce atan2 reagovat, pokud se náhodou strefíme do pólu a jak x tak y budou 0? Bude hrozit dělení nulou, a navíc výsledek bude mnohoznačný. Vrátí atan2 NaN? Naštěstí nevrátí. Je nám jedno, co vrátí za zeměpisnou délku, protože na tu se na pólu nehledí, ale chci aby nějaké číslo vrátila:
> atan2(0,0) [1] 0
Do souboru randonaut_3d.R si napíšeme nový program, tentokrát si dáme směrodatnou odchylku 5 tisíc kilometrů. Spustíme příkazemR --no-save < randonaut_3d.R v shellu:
home=c(50.0917174, 14.4055056)
som=5e6
zss=home[1]
zds=home[2]
pi=3.1415926535897932384626433
snr=pi/180
zdr=zds*snr
zsr=zss*snr
voszm=6371000
xm=voszm*cos(zdr)*cos(zsr)
ym=voszm*sin(zdr)*cos(zsr)
zm=voszm*sin(zsr)
xm=xm+rnorm(1,,som)
ym=ym+rnorm(1,,som)
zm=zm+rnorm(1,,som)
rns=180/pi
voszm=sqrt(xm*xm+ym*ym+zm*zm)
zdr=atan2(ym,xm)
zsr=asin(zm/voszm)
zds=zdr*rns
zss=zsr*rns
sprintf("(%.7f, %.7f)" , zss, zds)
sprintf("https://www.google.com/maps/search/?api=1&query=%.7f,%.7f" , zss, zds)
sprintf("https://openstreetmap.org/?mlat=%.7f&mlon=%.7f" , zss, zds)
Ukázkové výstupy, tentokrát daleko od Prahy:
[1] "(51.9016328, -8.9872039)" [1] https://www.google.com/maps/search/?api=1&query=51.9016328,-8.9872039 [1] https://openstreetmap.org/?mlat=51.9016328&mlon=-8.9872039 [1] "(46.7804985, -2.1882695)" [1] https://www.google.com/maps/search/?api=1&query=46.7804985,-2.1882695 [1] https://openstreetmap.org/?mlat=46.7804985&mlon=-2.1882695 [1] "(48.8423209, -23.8639521)" [1] https://www.google.com/maps/search/?api=1&query=48.8423209,-23.8639521 [1] https://openstreetmap.org/?mlat=48.8423209&mlon=-23.8639521
Pokud chceme generovat body na Zemi úplně rovnoměrně rozložené, v programu před přičítáním náhodných hodnot všechny 3 souřadnice xm, ym a zm přepíšeme nulami – bod se tak posune do středu Země a (nenulová) hodnota proměnné som nebude mít na výsledek vliv:
xm=ym=zm=0 xm=xm+rnorm(1,,som) ym=ym+rnorm(1,,som) zm=zm+rnorm(1,,som)
Následují ukázkové výstupy. 71 % povrchu Země je pokryto oceánem, takže je normální, že většina bodů bude ve vodě:
[1] "(12.8138031, -102.9801809)" [1] https://www.google.com/maps/search/?api=1&query=12.8138031,-102.9801809 [1] https://openstreetmap.org/?mlat=12.8138031&mlon=-102.9801809 [1] "(5.9652654, 105.2035492)" [1] https://www.google.com/maps/search/?api=1&query=5.9652654,105.2035492 [1] https://openstreetmap.org/?mlat=5.9652654&mlon=105.2035492 [1] "(-46.8074627, -142.2256476)" [1] https://www.google.com/maps/search/?api=1&query=-46.8074627,-142.2256476 [1] https://openstreetmap.org/?mlat=-46.8074627&mlon=-142.2256476
Pokud se rozhodnete s GNU R na náhodná místa cestovat, přeju hodně zábavy a nějaký ten bizarní nález!
ČLÁNKY DO MAILU