J’ai pris la carte ETOPO1 de la NOAA qui donne la topographie de la surface de la Terre en fonction de la longitude et de la latitude (plus exactement celle qui prend la surface des glaces et pas du bedrock, étant donné que pour l’instant, si le doigt tombe en Arctique, on n’est pas dans l’océan même si c’est le cas topographiquement). Il y a un point par minute d’arc, donc une fois décompressé, c’est un fichier de 4Go quand même.
Puis j’ai considéré une densité de proba du point choisi indépendante de la longitude et forte vers l’équateur. Libre à vous de compliquer ça par la suite, autant commencer simple. Ma densité s’exprime $f(\theta)=\dfrac{3}{8\pi}\cos^2\theta$. Le facteur est là pour assurer que l’intégrale sur la sphère $\int f\mathrm d\Omega=\int_{\theta =-\pi/2}^{\pi/2}\int_{\phi = 0}^{2\pi} f\cos\theta\mathrm d\theta\mathrm d\varphi=1$.
La probabilité recherchée est alors simplement $\int o\times f\mathrm d\Omega$ où $o$ est une fonction de $\theta$ et $\phi$ valant $1$ si le point est dans l’océan et $0$ sinon. Comme critère, j’ai simplement pris qu’une altitude négative signifie que l’on est dans l’océan. C’est pas tout à fait vrai pour plein de raisons, mais c’est pas mal pour une première approximation.
J’ai écrit un code en Fortran (Python est à la ramasse sur un fichier aussi gros…) qui fait ce calcul :
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29 | program ocean
implicit none
integer:: unt, i
integer, parameter:: ndata = 233312401
integer, parameter:: output = 1000000
real:: lon, lat, topo
real:: theta, phi, proba, ddpint, delta
real, parameter:: pi = 4 * atan(1.)
real, parameter:: one_deg = pi / 180
real, parameter:: norm = 3 / (8 * pi)
real, parameter:: dtheta_dphi = (one_deg / 60) ** 2
proba = 0
ddpint = 0
open(newunit=unt, file='ETOPO1_Ice_g_int.xyz/data', status='old')
do i = 1, ndata
read(unt, *) lon, lat, topo
theta = lat * one_deg
delta = norm * cos(theta)**3 * dtheta_dphi
proba = proba + delta * merge(1, 0, topo <= 0)
ddpint = ddpint + delta
if (mod(i, output) == 0) print *, proba, ddpint, proba/ddpint
end do
close(unt)
print *, 'Finished:'
print *, proba, ddpint, proba/ddpint
end program ocean
|
Là-dedans, proba
représente ce que l’on cherche. ddpint
l’intégrale de la densité de proba (pour vérifier que cela fait bien $1$ numériquement) et je calcule aussi le quotient entre les deux proba/ddpint
pour se donner une idée d’à quel point l’erreur sur ddpint
affecte proba
.
Sur ma machine, après compilation avec gfortran -fdefault-real-8 -O3 probaOcean.f90 -o ocean
, j’obtiens pour ces trois variables les valeurs suivantes :
| 0.73560247442935078 1.0000462961701353 0.73556842042861248
|
on voit déjà que ddpint
est très proche de $1$, surtout compte tenu du nombre important de points sur lesquels l’intégration est effectuée (plus de 233 millions de points). La probabilité de tomber dans un océan est légèrement plus forte qu’avec le calcul basé sur une distribution uniforme de LudoBike, ce qui est logique puisqu’il y de fortes masses continentales près des pôles qui sont ici défavorisées par ma densité.
Pour rigoler, j’ai aussi calculé la probabilité de tomber sur une zone de haute montagne avec comme critère une topo supérieure à 2000 m. Les résultats sont pas folichons, avec une proba inférieure à 2%.
| 1.7484948079127604E-002 1.0000462961701353 1.7484138630470897E-002
|