Processing a project

The quick reading functionality of resistics focuses on analysis of single continuous recordings. When there are multiple recordings at a site or multiple sites, it can be more convenient to use a resistics project. This is generally easier to manage and use, especially when doing remote reference or intersite processing.

The data in this example has been provided for use by the SAMTEX consortium. For more information, please refer to [Jones2009]. Additional details about the data can be found at https://www.mtnet.info/data/kap03/kap03.html.

from pathlib import Path
import seedir as sd
import plotly
import resistics.letsgo as letsgo

Let’s remind ourselves of the project contents and then load the project.

project_path = Path("..", "..", "data", "project", "kap03")
sd.seedir(str(project_path), style="emoji")
resenv = letsgo.load(project_path)

Out:

📁 kap03/
├─📁 images/
├─📁 time/
│ ├─📁 kap163/
│ │ └─📁 meas01/
│ │   ├─📄 data.npy
│ │   └─📄 metadata.json
│ ├─📁 kap160/
│ │ └─📁 meas01/
│ │   ├─📄 data.npy
│ │   └─📄 metadata.json
│ └─📁 kap172/
│   └─📁 meas01/
│     ├─📄 data.npy
│     └─📄 metadata.json
├─📄 resistics.json
├─📁 results/
├─📁 calibrate/
├─📁 evals/
├─📁 masks/
├─📁 spectra/
└─📁 features/

Inspect the current configuration. As no custom configuration has been specified, this will be the default configuration.

resenv.config.summary()

Out:

{
    'name': 'default',
    'time_readers': [
        {
            'name': 'TimeReaderAscii',
            'apply_scalings': True,
            'extension': '.txt',
            'delimiter': None,
            'n_header': 0
        },
        {
            'name': 'TimeReaderNumpy',
            'apply_scalings': True,
            'extension': '.npy'
        }
    ],
    'time_processors': [
        {'name': 'InterpolateNans'},
        {'name': 'RemoveMean'}
    ],
    'dec_setup': {
        'name': 'DecimationSetup',
        'n_levels': 8,
        'per_level': 5,
        'min_samples': 256,
        'div_factor': 2,
        'eval_freqs': None
    },
    'decimator': {
        'name': 'Decimator',
        'resample': True,
        'max_single_factor': 3
    },
    'win_setup': {
        'name': 'WindowSetup',
        'min_size': 128,
        'min_olap': 32,
        'win_factor': 4,
        'olap_proportion': 0.25,
        'min_n_wins': 5,
        'win_sizes': None,
        'olap_sizes': None
    },
    'windower': {'name': 'Windower'},
    'fourier': {
        'name': 'FourierTransform',
        'win_fnc': ['kaiser', 14],
        'detrend': 'linear',
        'workers': -2
    },
    'spectra_processors': [],
    'evals': {'name': 'EvaluationFreqs'},
    'sensor_calibrator': {
        'name': 'SensorCalibrator',
        'chans': None,
        'readers': [
            {
                'name': 'SensorCalibrationJSON',
                'extension': '.json',
                'file_str': 'IC_$sensor$extension'
            }
        ]
    },
    'tf': {
        'name': 'ImpedanceTensor',
        'variation': 'default',
        'out_chans': ['Ex', 'Ey'],
        'in_chans': ['Hx', 'Hy'],
        'cross_chans': ['Hx', 'Hy'],
        'n_out': 2,
        'n_in': 2,
        'n_cross': 2
    },
    'regression_preparer': {'name': 'RegressionPreparerGathered'},
    'solver': {
        'name': 'SolverScikitTheilSen',
        'fit_intercept': False,
        'normalize': False,
        'n_jobs': -2,
        'max_subpopulation': 2000,
        'n_subsamples': None
    }
}

And it’s always useful to know what transfer function will be calculated out.

print(resenv.config.tf)

Out:

| Ex | = | Ex_Hx Ex_Hy | | Hx |
| Ey |   | Ey_Hx Ey_Hy | | Hy |

Now let’s run single site processing on a site and then look at the directory structure again. Begin by transforming to frequency domain and reducing to the evaluation frequencies. Note that whilst there is only a single measurement for this site, the below is written to work when there are more measurements.

site = resenv.proj["kap160"]
for meas in site:
    letsgo.process_time_to_evals(resenv, "kap160", meas.name)
sd.seedir(str(project_path), style="emoji")

Out:

📁 kap03/
├─📁 images/
├─📁 time/
│ ├─📁 kap163/
│ │ └─📁 meas01/
│ │   ├─📄 data.npy
│ │   └─📄 metadata.json
│ ├─📁 kap160/
│ │ └─📁 meas01/
│ │   ├─📄 data.npy
│ │   └─📄 metadata.json
│ └─📁 kap172/
│   └─📁 meas01/
│     ├─📄 data.npy
│     └─📄 metadata.json
├─📄 resistics.json
├─📁 results/
├─📁 calibrate/
├─📁 evals/
│ └─📁 kap160/
│   └─📁 default/
│     └─📁 meas01/
│       ├─📄 data.npz
│       └─📄 metadata.json
├─📁 masks/
├─📁 spectra/
└─📁 features/

Now let’s run single site processing on a site and then look at the directory structure again. To run the transfer function calculation, the sampling frequency to process needs to be specified. In this case, it’s 0.2 Hz.

letsgo.process_evals_to_tf(resenv, 0.2, "kap160")
sd.seedir(str(project_path), style="emoji")

Out:

  0%|          | 0/20 [00:00<?, ?it/s]
 15%|#5        | 3/20 [00:00<00:00, 28.86it/s]
 45%|####5     | 9/20 [00:00<00:00, 45.49it/s]
100%|##########| 20/20 [00:00<00:00, 88.57it/s]

  0%|          | 0/20 [00:00<?, ?it/s]
  5%|5         | 1/20 [00:01<00:29,  1.56s/it]
 10%|#         | 2/20 [00:03<00:28,  1.56s/it]
 15%|#5        | 3/20 [00:04<00:26,  1.56s/it]
 20%|##        | 4/20 [00:06<00:24,  1.56s/it]
 25%|##5       | 5/20 [00:07<00:23,  1.58s/it]
 30%|###       | 6/20 [00:08<00:16,  1.20s/it]
 35%|###5      | 7/20 [00:08<00:12,  1.04it/s]
 40%|####      | 8/20 [00:09<00:09,  1.25it/s]
 45%|####5     | 9/20 [00:09<00:07,  1.44it/s]
 50%|#####     | 10/20 [00:10<00:06,  1.61it/s]
 55%|#####5    | 11/20 [00:10<00:04,  2.10it/s]
 60%|######    | 12/20 [00:10<00:02,  2.68it/s]
 65%|######5   | 13/20 [00:10<00:02,  3.30it/s]
 70%|#######   | 14/20 [00:10<00:01,  3.95it/s]
 75%|#######5  | 15/20 [00:10<00:01,  4.56it/s]
 80%|########  | 16/20 [00:10<00:00,  5.32it/s]
 85%|########5 | 17/20 [00:11<00:00,  6.08it/s]
 90%|######### | 18/20 [00:11<00:00,  6.75it/s]
 95%|#########5| 19/20 [00:11<00:00,  7.31it/s]
100%|##########| 20/20 [00:11<00:00,  7.74it/s]
100%|##########| 20/20 [00:11<00:00,  1.75it/s]
📁 kap03/
├─📁 images/
├─📁 time/
│ ├─📁 kap163/
│ │ └─📁 meas01/
│ │   ├─📄 data.npy
│ │   └─📄 metadata.json
│ ├─📁 kap160/
│ │ └─📁 meas01/
│ │   ├─📄 data.npy
│ │   └─📄 metadata.json
│ └─📁 kap172/
│   └─📁 meas01/
│     ├─📄 data.npy
│     └─📄 metadata.json
├─📄 resistics.json
├─📁 results/
│ └─📁 kap160/
│   └─📁 default/
│     └─📄 0_200000_impedancetensor_default.json
├─📁 calibrate/
├─📁 evals/
│ └─📁 kap160/
│   └─📁 default/
│     └─📁 meas01/
│       ├─📄 data.npz
│       └─📄 metadata.json
├─📁 masks/
├─📁 spectra/
└─📁 features/

Get the transfer function

soln = letsgo.get_solution(
    resenv,
    "kap160",
    resenv.config.name,
    0.2,
    resenv.config.tf.name,
    resenv.config.tf.variation,
)
fig = soln.tf.plot(
    soln.freqs,
    soln.components,
    to_plot=["ExHy", "EyHx"],
    x_lim=[1, 5],
    res_lim=[1, 4],
    legend="128",
    symbol="circle",
)
fig.update_layout(height=900)
plotly.io.show(fig)