Fonctions utiles
- Matrice hessienne du champ magnétique dipolaire
- Calculer cosinus directeurs du champ magnétique
- Réduction au pôle des anomalies magnétiques
- Faire la séparation régionale–résiduelle
- Calculer l’intégrand de Talwani
- Importer un array numpy du web
- Importer une image du web
Cette page rassemble des fonctions Python utiles pour le cours et les travaux pratiques.
Matrice hessienne du champ magnétique dipolaire
Cette matrice est le noyau qui permet de calculer le champ magnétique d’un dipole en coordonnées cartésiennes (Section 4.3.4).
import numpy as np
def calculer_J(x, y, z):
# Retourne la matrice J
r = np.sqrt(x*x + y*y + z*z)
J = np.zeros((3, 3))
J[0, 0] = 3*x**2 / r**5 - 1/r**3
J[1, 1] = 3*y**2 / r**5 - 1/r**3
J[2, 2] = 3*z**2 / r**5 - 1/r**3
J[0, 1] = 3*x*y / r**5
J[0, 2] = 3*x*z / r**5
J[1, 0] = 3*x*y / r**5
J[1, 2] = 3*y*z / r**5
J[2, 0] = 3*x*z / r**5
J[2, 1] = 3*y*z / r**5
return J
# Application
xs = np.linspace(-20, 20, 200) # m
y = 0 # m
z = -5 # m
for x in xs:
J = calculer_J(x, y, z)
...
Calculer cosinus directeurs du champ magnétique
Cette fonction permet d’obtenir la direction de l’aimantation induite par un champ magnétique dont l’inclinaison et la déclinaison sont connues.
import numpy as np
def calculer_direction_m(I, D):
# Retourne un vecteur horizontal des
# cosinus directeurs de l'aimantation
I, D = np.radians(I), np.radians(D)
return np.array([np.cos(I)*np.cos(D),
np.cos(I)*np.sin(D),
np.sin(I)])
# Application
I = 60 # degrés
D = 15 # degrés
m = calculer_direction_m(I, D)
Réduction au pôle des anomalies magnétiques
Cette fonction projette les anomalies magnétiques comme si elles avaient été mesurées au pôle nord magnétique (Section 4.4.1).
import numpy as np
def reduction_pole(T, I, D, step):
"""Calcule la réduction au pôle d'une carte d'anomalies magnétiques.
Args:
T (2d array): Anomalie de champ total (m x n).
I (float): Inclinaison du champ en degrés.
D (float): Déclinaison du champ en degrés.
step (tuple): Pas (x, y) entre les stations en mètres.
Returns:
(2d array): L'anomalie réduite au pôle.
"""
I_rad, D_rad = np.radians(I), np.radians(D) # conversion en radians
T_hat = np.fft.fft2(T) # transformation Fourier
k_x = np.fft.fftfreq(n=T.shape[0], d=step[0]) # calcul kx
k_y = np.fft.fftfreq(n=T.shape[1], d=step[1]) # calcul ky
KXX, KYY = np.meshgrid(k_x, k_y) # mettre kx et ky en grille
mu = np.arctan2(KYY, KXX)
psi = (np.sin(I_rad) + 1j*np.cos(I_rad)*np.cos(D_rad - mu))**-2
T_hat *= psi # filtrage
T_pole = np.fft.ifft2(T_hat).real # transformation inverse
return T_pole
# Application
x, y = np.meshgrid(np.linspace(-1, 1, 11), np.linspace(-1, 1, 11))
d = np.sqrt(x*x + y*y)
sigma, mu = 1.0, 0.0
T = np.exp(-( (d-mu)**2 / ( 2.0 * sigma**2 ) ) ) # données synthétiques
I = 80
D = 0
T_pole = reduction_pole(T, I, D, step=(0.2, 0.2))
Faire la séparation régionale–résiduelle
Cette fonction permet de calculer l’anomalie régionale à partir d’une anomalie de Bouguer et d’en extraire la composante résiduelle (Section 3.4).
import numpy as np
def separer_reg_res(gz, h, step):
"""Calcule le prolongement du champ `gz` à une hauteur `h`.
Args:
gz (2d array): Anomalie de Bouguer (m x n).
h (int): Hauteur du prolongement en mètres.
step (tuple): Pas (x, y) entre les stations en mètres.
Returns:
tuple: La régionale et la résiduelle.
"""
gz_hat = np.fft.fft2(gz) # transformation Fourier
k_x = np.fft.fftfreq(n=gz.shape[0], d=step[0]) # calcul kx
k_y = np.fft.fftfreq(n=gz.shape[1], d=step[1]) # calcul ky
KXX, KYY = np.meshgrid(k_x, k_y) # mettre kx et ky en grille
gz_hat *= np.exp(-h * np.sqrt(KXX*KXX + KYY*KYY)) # filtrage
regionale = np.fft.ifft2(gz_hat).real # transformation inverse
residuelle = gz - regionale # séparation rés-rég
return regionale, residuelle
# Application
x, y = np.meshgrid(np.linspace(-1, 1, 11), np.linspace(-1, 1, 11))
d = np.sqrt(x*x + y*y)
sigma, mu = 1.0, 0.0
g = np.exp(-( (d-mu)**2 / ( 2.0 * sigma**2 ) ) ) # données synthétiques
reg, res = separer_reg_res(gz=g, h=10, step=(0.2, 0.2))
Calculer l’intégrand de Talwani
Cette fonction permet de modéliser la réponse gravimétrique d’un polygone arbitraire d’extension infinie (Section 2.4.5).
import numpy as np
def integrand_talwani(x1, x2, z1, z2, density):
"""Suppose que l'intégration se fait en sens horaire.
Args:
x1, x2, z1, z2 (float): coordonnées (x, z) de deux
sommets consécutifs du polygone en sens horaire
density (float): contraste de densité du polygone
Returns:
float: l'intégrand de Talwani
"""
# Calcul des composantes
denom = z2 - z1
alpha = (x2 - x1) / denom
beta = (x1 * z2 - x2 * z1) / denom
factor = beta / (1 + alpha * alpha)
r1sq = (x1 * x1 + z1 * z1) # r1^2
r2sq = (x2 * x2 + z2 * z2) # r2^2
term1 = 0.5 * (np.log(r2sq) - np.log(r1sq)) # pour éviter de calculer ln(r2^2/r1^2)
term2 = np.arctan2(z2, x2) - np.arctan2(z1, x1) # theta2 - theta1
return factor * (term1 - alpha * term2)
# Application
intg_list = [] # liste vide pour mettre l'integrand de chaque arête du polygone
intg_list.append(integrand_talwani(1, 2, -1, -2, density=-2600)) # calcul
gz = np.sum(intg_list) # sommer les integrands pour faire l'integrale de contour
Importer un array numpy du web
Cette fonction permet de charger un vecteur numpy à partir d’une adresse donnée dans les notes de cours.
import requests
from io import BytesIO
import numpy as np
def np_load_from_url(url):
address = BytesIO(requests.get(url).content)
return np.loadtxt(address)
# Application
url = 'https://clberube.org/glq2200/assets/data/gravi_aero.txt'
data = np_load_from_url(url)
Importer une image du web
Cette fonction permet de charger une image à partir d’une adresse donnée dans les notes de cours.
import requests
from io import BytesIO
import numpy as np
from PIL import Image
def load_img_from_url(url):
img = Image.open(BytesIO(requests.get(url).content))
return np.array(img)
# Application
url = 'http://sipi.usc.edu/database/download.php?vol=misc&img=5.2.10'
data = load_img_from_url(url)