Fonctions utiles

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)

Retour en haut

© 2020–2024 Charles L. Bérubé
Polytechnique Montréal