Hlavní navigace

Randonautica v GNU R: kvalitní kulově symetrická metoda

25. 2. 2021
Doba čtení: 5 minut

Sdílet

 Autor: Depositphotos
V dnešním článku si naprogramujeme kvalitní kulově symetrickou metodu, která nebude mít nedostatky jednoduché metody z minulého dílu. Pokud se rozhodnete s GNU R na náhodná místa cestovat, přeji hodně zábavy.

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ů.


Funkce nejsou moc složité.

Úhly:
  • 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á.

Joann94024 at English Wikipedia, CC BY-SA 3.0

Osa X vede od Howlandova Ostrovu (šířka 0, délka +/-180) …


Michal Vogt, CC BY-SA 2.0, Via Wikimedia Commons

… ke Ghaně (šířka 0, délka 0).

  • 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

David Adam Kess CC BY-SA 3.0

Osa Y vede od Galapág (šířka 0, délka –90)…


AbdurrahmanPandia CC BY-SA 4.0, Via Wikimedia Commons

… k Sumatře (šířka 0, délka +90).

  • 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

Daniel Case CC BY-SA 3.0

Osa Z vede od Jižního pólu (šířka –90, délka libovolná)…


Public domain. Via Wikimedia Commons

… 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:

Hacking tip

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!

Autor článku

Karel Kulhavý vystudoval operační systémy, sítě a překladače na MFF UK a je autorem optického pojítka Twibright Ronja a spoluautorem textového a grafického webového prohlížeče Twibright Links.