Link Search Menu Expand Document

Fonctions utiles


Cette page est une collection de fonctions Python pertinentes 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

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

import numpy as np

def reduction_pole(T, I, D, step):
    """Calcule le prolongement du champ `gz` à une hauteur `h`.
    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

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

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

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

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-2021 Charles L. Bérubé
Polytechnique Montréal