Observations covariances, Localisation matrix

The notebook illustrates how to read an EM sensors dataset to prepare it for DA.

Estimated time to run the notebook = 2min

import os
from pathlib import Path

import matplotlib.pyplot as plt
import pandas as pd

from pyCATHY.DA.cathy_DA import DA, dictObs_2pd, localisation
from pyCATHY.DA.observations import read_observations

rootpath = Path(os.getcwd()).parent
fwdEM_path = rootpath / 'EM_forward'        # folder containing the emagpy outputs
SolfwdName = 'my_forward_solution'          # sub-folder of the forward solution
data_err = 0.05                             # instrumental error
selected_indices = [0, 10, 20, 30, 40, 50]  # indices of the EM surveys to read

---

simuWithDA = DA(
    dirName='.',
    prj_name='test_covariance_localisation',
    notebook=True,
)

simuWithDA.run_preprocessor()
simuWithDA.run_processor(IPRT1=3)

🏁 Initiate CATHY object
🍳 gfortran compilation
👟 Run preprocessor
🔄 Update parm file
🔄 Update hap.in file
🔄 Update dem_parameters file
🔄 Update dem_parameters file
🛠  Recompile src files [4s]
🍳 gfortran compilation [11s]
✅ Compilation successful!
👟 Run processor

---

dict_obs = {}  # initiate the dictionary

# EM observations metadata
metadata_EM = {
    'data_type': 'EM',
    'units': r'$mS/m$',        # raw string avoids the invalid '\m' escape
    'instrument': 'GF6L',
    'data_format': 'emagpy',
    'fwdEMnoise': 5,
}

start_date = pd.to_datetime('2023-01-01 08:00:00.00000')
datetime_EM = pd.date_range(start=start_date, periods=60, freq='D')  # 60 daily points
elapsed_seconds = (datetime_EM - start_date).total_seconds()

for i in selected_indices:
    tt = elapsed_seconds[i]
    # .txt file observation generated by emagpy
    filename = fwdEM_path / SolfwdName / f'fwd_ECa{i}.txt'

    read_observations(
        dict_obs,
        filename,
        data_type='EM',
        data_err=data_err,
        show=True,
        tA=tt,
        obs_cov_type='data_err',
        meta=metadata_EM,
        datetime=datetime_EM[i],
    )

data_measure_df = dictObs_2pd(dict_obs)
print(data_measure_df.head())
print(data_measure_df.iloc[0]['data'])

Traceback (most recent call last):
  File "/home/z0272571a@campus.csic.es/Nextcloud/BenCSIC/Codes/BenjMy/pycathy_wrapper/examples/DA/plot_5_covariances.py", line 61, in <module>
    read_observations(
  File "/home/z0272571a@campus.csic.es/Nextcloud/BenCSIC/Codes/BenjMy/pycathy_wrapper/pyCATHY/DA/observations.py", line 188, in read_observations
    df, k = in_meas.read_EM(obs_2_add, data_format)
  File "/home/z0272571a@campus.csic.es/Nextcloud/BenCSIC/Codes/BenjMy/pycathy_wrapper/pyCATHY/importers/sensors_measures.py", line 42, in read_EM
    k = Problem()
NameError: name 'Problem' is not defined

---

nnod = int(simuWithDA.grid3d['nnod'])
grid_coords = simuWithDA.grid3d['mesh3d_nodes'][:nnod, :2]

ncoils = 6
L = 4

# Build localisation matrix (with coil covariance)
local_mat = localisation.build_localization_matrix(
    grid_coords,
    ncoils,
    L,
    with_coil_covariance=True,
)

# Visualise a small block
plt.figure(figsize=(6, 6))
plt.imshow(local_mat[:60, :60], cmap='viridis')
plt.colorbar(label='Localization weight')
plt.title('Localization for Coils at the Same Grid Point')
plt.show()
print("Localization matrix shape:", local_mat.shape)

# One localisation matrix per observation time
simuWithDA.localisationMatrix = [local_mat] * len(data_measure_df.index)