TP4 - Levé magnétique sur le Mont-Royal

À rendre le 2023-04-19 - 7% du cours

---

• Contexte
• Objectifs
  • Objectifs généraux
  • Objectifs spécifiques
• Données
• Exemples pour démarrer
  • Librairies
  • Carte de couleur
  • Calcul de la réduction au pôle

Contexte

Le TP4 porte sur la collecte et le traitement des données magnétiques. La première séance de TP est consacrée à la réalisation d’un levé magnétique sur le Mont-Royal, dans le secteur du Lac-aux-Castors. La deuxième séance de TP est consacrée au traitement et à l’interprétation des données magnétiques.

Objectifs

Objectifs généraux

1. Comprendre les enjeux logistiques liés à la réalisation d’un levé magnétique sur le terrain.
2. Comprendre et discuter des différentes étapes du traitement des données magnétiques.

Objectifs spécifiques

1. Réaliser un levé magnétique sur le secteur du Lac-aux-Castors.
2. Calculer la réduction au pôle des données magnétiques.
3. À partir de la réduction au pôle, séparer l’anomalie régionale de l’anomalie résiduelle.
4. Calculer la dérivée verticale de l’anomalie résiduelle.
5. Interpréter les données magnétiques selon vos connaissances de la géologie du Mont-Royal.

Données

Veuillez télécharger les données du levé gravimétrique ici : TP4-WALKMAG-LC.csv. Les coordonnées sont fournies en mètre et l’anomalie de champ total est donnée en nT. L’anomalie de champ total est montrée à la Figure 1 dans le système de coordonnées est WGS 84 / UTM zone 18N (EPSG:32618).

Figure 1. Carte magnétique du secteur du Lac-aux-Castors (Brisebois & Bérubé, 2022).

Les données du levé magnétique ont déjà été corrigées pour la dérive, puis interpolées sur une grille uniforme rectangulaire. Sur la grille fournie, il y a 103 points de mesures dans la direction $x$ et 155 points dans la direction $y$. L’espacement entre les mailles de la grille est de 3 m dans les deux directions.

Exemples pour démarrer

Librairies

import numpy as np  # calculs numériques
import pandas as pd  # gestion des données
import matplotlib.pyplot as plt  # graphiques

Carte de couleur

Attention avec l’utilisation de la fonction .reshape() pour la mise en carte des données lorsque la grille est rectangulaire. Comme le premier indice dans une matrice numpy désigne les rangées de la matrice, le premier indice est celui de la coordonnée $y$. Le deuxième indice est celui de la coordonnée $x$.

# Notez le (155, 103) et non (103, 155) pour 
# obtenir un rectangle avec le nord vers le haut.
plt.pcolormesh(df['X'].values.reshape(155, 103),  
               df['Y'].values.reshape(155, 103), 
               df['g_obs_SI'].values.reshape(155, 103), 
               cmap='magma')

Calcul de la réduction au pôle

Cette fonction permet de calculer la réduction au pôle des données magnétiques. Comme c’est une opération de filtrage dans le domaine spectral, son utilisation est très semblable à celle des fonctions définies dans le TP3 pour filtrer les données gravimétriques.

def reduction_pole(T, x, y, I, D, pad=10):
    """Calcule la réduction au pôle d'une carte d'anomalies magnétiques.
    Args:
        T (1d array): Vecteur de l'anomalie de champ total.
        x (1d array): Vecteur des coordonnées X.
        y (1d array): Vecteur des coordonnées Y.
        I (float): Inclinaison du champ en degrés.
        D (float): Déclinaison du champ en degrés.
        pad (int, optionnel): Nombre de pixels de padding.
    Returns:
        (1d array): L'anomalie réduite au pôle.
    """
    # Préparation des données et padding
    I_rad, D_rad = np.radians(I), np.radians(D)  # conversion en radians
    hw_shape = (len(y.unique()), len(x.unique()))
    T = T.values.reshape(hw_shape)
    T = np.pad(T, (pad, pad), 'reflect')
    # Transformation de Fourier
    T_hat = np.fft.fft2(T)
    k_x = np.fft.fftfreq(n=T.shape[1], d=x.diff().abs().max())  # calcul kx
    k_y = np.fft.fftfreq(n=T.shape[0], d=y.diff().abs().max())  # calcul ky
    KXX, KYY = np.meshgrid(k_x, k_y)  # mettre kx et ky en grille
    # Opération de filtrage
    mu = np.arctan2(KYY, KXX)
    psi = (np.sin(I_rad) + 1j*np.cos(I_rad)*np.cos(D_rad - mu))**-2
    T_hat *= psi  
    # Transformation inverse et retrait du padding
    T_pole = np.fft.ifft2(T_hat).real
    T_pole = T_pole[pad:-pad, pad:-pad].flatten()
    return T_pole

Dans l’exemple d’utilisation suivant, attention de remplacer l’inclinaison (I) et la déclinaison (D) par les bonnes valeurs observées pour le champ magnétique à Montréal :

df['Pole'] = reduction_pole(df['Delta_T'], df['X'], df['Y'], I=None, D=None)

Pour tenir compte du fait que la grille est rectangulaire et non carrée, il y a des différences entre la fonction de réduction au pôle ci-haut et celles du TP3. Pour les autres filtres, comme le prolongement vers le haut, veuillez utiliser le présent gabarit et modifier seulement l’opération de filtrage.