Full sEIT data processing example

Full processing of one timestep of the data from Weigand and Kemna 2017 (Biogeosciences).

imports

import os
import subprocess

import reda
from reda.utils.fix_sign_with_K import fix_sign_with_K
import reda.utils.geometric_factors as geom_facs
import numpy as np
import reda.importers.eit_fzj as eit_fzj

data import

seit = reda.sEIT()
seit.import_eit_fzj(
    'data/bnk_raps_20130408_1715_03_einzel.mat',
    'data/configs.dat'
)
print(seit.data[['a', 'b', 'm', 'n']].iloc[0:10])
# In order to keep things fast, we only keep the following frequencies
seit.keep_frequencies((1, 10, 100, 1000, 10000))

Constructing four-point measurements
Summary:
                  a             b  ...     frequency          rpha
count  67375.000000  67375.000000  ...  67375.000000  67375.000000
mean      17.907532     20.831169  ...   2577.580109    473.723981
std       10.333167     10.677015  ...   8258.716055   2729.735973
min        1.000000      2.000000  ...      0.462963  -3141.592649
25%       10.000000     12.000000  ...     29.411765  -3131.399200
50%       17.000000     22.000000  ...    122.222220      5.510219
75%       27.000000     30.000000  ...    700.000000   3139.297772
max       37.000000     38.000000  ...  45000.000000   3141.592652

[8 rows x 7 columns]
    a   b   m   n
17  1   2  20  21
47  1  27  30  28
48  1  27  31  29
49  1  27  32  28
50  1  27  32  30
73  2   3  20  21
86  2   3  32   1
93  2   4   5   8
94  2   4   6   9
95  2   4   7  10

compute geometric factors and correct for signs/phase shifts by pi

settings = {
    'rho': 100,
    'elem': 'data/elem.dat',
    'elec': 'data/elec.dat',
    'sink_node': '6467',
    '2D': True,
}
k = geom_facs.compute_K_numerical(seit.data, settings)
seit.data = geom_facs.apply_K(seit.data, k)

# input('nr 2, press enter to continue')
fix_sign_with_K(seit.data)

SETTINGS
{'rho': 100, 'elem': 'data/elem.dat', 'elec': 'data/elec.dat', 'sink_node': '6467', '2D': True}
2D modeling

apply correction factors for 2D rhizotron tank

corr_facs_nor = np.loadtxt('data/corr_fac_avg_nor.dat')
corr_facs_rec = np.loadtxt('data/corr_fac_avg_rec.dat')
corr_facs = np.vstack((corr_facs_nor, corr_facs_rec))
seit.data, cfacs = eit_fzj.apply_correction_factors(seit.data, corr_facs)

apply data filters

seit.filter('r < 0')
seit.filter('rho_a < 15 or rho_a > 35')
seit.filter('rpha < - 40 or rpha > 3')
seit.filter('rphadiff < -5 or rphadiff > 5')
seit.filter('k > 400')
seit.filter('rho_a < 0')
seit.filter('a == 12 or b == 12 or m == 12 or n == 12')
seit.filter('a == 13 or b == 13 or m == 13 or n == 13')

# import IPython
# IPython.embed()
seit.print_data_journal()
seit.filter_incomplete_spectra(flimit=300, percAccept=85)
seit.print_data_journal()
seit.print_log()

--- Data Journal Start ---
2024-11-14 10:34:05.382368
A filter was applied with query "r < 0". In total 75 records were removed
A filter was applied with query "rho_a < 15 or rho_a > 35". In total 1963 records were removed
A filter was applied with query "rpha < - 40 or rpha > 3". In total 1570 records were removed
A filter was applied with query "rphadiff < -5 or rphadiff > 5". In total 2016 records were removed
A filter was applied with query "k > 400". In total 557 records were removed
A filter was applied with query "rho_a < 0". In total 0 records were removed
A filter was applied with query "a == 12 or b == 12 or m == 12 or n == 12". In total 141 records were removed
A filter was applied with query "a == 13 or b == 13 or m == 13 or n == 13". In total 130 records were removed
--- Data Journal End ---


--- Data Journal Start ---
2024-11-14 10:34:05.594825
A filter was applied with query "r < 0". In total 75 records were removed
A filter was applied with query "rho_a < 15 or rho_a > 35". In total 1963 records were removed
A filter was applied with query "rpha < - 40 or rpha > 3". In total 1570 records were removed
A filter was applied with query "rphadiff < -5 or rphadiff > 5". In total 2016 records were removed
A filter was applied with query "k > 400". In total 557 records were removed
A filter was applied with query "rho_a < 0". In total 0 records were removed
A filter was applied with query "a == 12 or b == 12 or m == 12 or n == 12". In total 141 records were removed
A filter was applied with query "a == 13 or b == 13 or m == 13 or n == 13". In total 130 records were removed
--- Data Journal End ---

2024-11-14 10:34:05,353 - reda.main.logger - INFO - Data sized changed from 9625 to 9550
2024-11-14 10:34:05,357 - reda.main.logger - INFO - Data sized changed from 9550 to 7587
2024-11-14 10:34:05,361 - reda.main.logger - INFO - Data sized changed from 7587 to 6017
2024-11-14 10:34:05,365 - reda.main.logger - INFO - Data sized changed from 6017 to 4001
2024-11-14 10:34:05,369 - reda.main.logger - INFO - Data sized changed from 4001 to 3444
2024-11-14 10:34:05,373 - reda.main.logger - INFO - Data sized changed from 3444 to 3444
2024-11-14 10:34:05,377 - reda.main.logger - INFO - Data sized changed from 3444 to 3303
2024-11-14 10:34:05,382 - reda.main.logger - INFO - Data sized changed from 3303 to 3173

import crtomo
grid = crtomo.crt_grid('data/elem.dat', 'data/elec.dat')

seitinv = seit.export_to_crtomo_seit_manager(grid, norrec='nor')
seitinv.crtomo_cfg['robust_inv'] = 'F'
seitinv.crtomo_cfg['mag_abs'] = 0.012
seitinv.crtomo_cfg['mag_rel'] = 0.5
seitinv.crtomo_cfg['hom_bg'] = 'T'
seitinv.crtomo_cfg['d2_5'] = 0
seitinv.crtomo_cfg['fic_sink'] = 'T'
seitinv.crtomo_cfg['fic_sink_node'] = 6467
seitinv.apply_crtomo_cfg()

This grid was sorted using CutMcK. The nodes were resorted!
Rectangular grid found
1.0 (295, 28) 8260
10.0 (295, 28) 8260
100.0 (295, 28) 8260
1000.0 (190, 28) 5320
10000.0 (102, 28) 2856

now run the inversion we do this the “old” style using the command td_run_all_local, which is also included in the crtomo_tools

# only invert if the sipdir does not already exist
if not os.path.isdir('sipdir'):
    # save to a sip-directory
    seitinv.save_to_eitdir('sipdir')

    os.chdir('sipdir')
    subprocess.call('td_run_all_local -t 1 -n 2', shell=True)
    os.chdir('..')

Now plot the results at this point all plot scripts for sEIT results are located in old, deprecated, packages, and thus we need to write a new plotting tools. In the mean time you can enter each tomodir in the sipdir/invmod subdirectory and plot it using the single-frequency plot command “td_plot”

# TODO

import crtomo
import numpy as np
sinv = crtomo.eitMan(seitdir='sipdir/')
# sinv.extract_points('rpha', np.atleast_2d(np.array((-0.5, 13))))
sinv.extract_points('rpha', np.atleast_2d(np.array((20, -5))))

This grid was sorted using CutMcK. The nodes were resorted!
Rectangular grid found
reading voltages
Reading inversion results
is robust False
Info: res_m.diag not found: sipdir//invmod/00_1.000000//inv/res_m.diag
/home/runner/work/crtomo_tools/crtomo_tools/examples/04_full_dataprocessing
Info: ata.diag not found: sipdir//invmod/00_1.000000//inv/ata.diag
/home/runner/work/crtomo_tools/crtomo_tools/examples/04_full_dataprocessing
reading voltages
Reading inversion results
is robust False
Info: res_m.diag not found: sipdir//invmod/01_10.000000//inv/res_m.diag
/home/runner/work/crtomo_tools/crtomo_tools/examples/04_full_dataprocessing
Info: ata.diag not found: sipdir//invmod/01_10.000000//inv/ata.diag
/home/runner/work/crtomo_tools/crtomo_tools/examples/04_full_dataprocessing
reading voltages
Reading inversion results
is robust False
Info: res_m.diag not found: sipdir//invmod/02_100.000000//inv/res_m.diag
/home/runner/work/crtomo_tools/crtomo_tools/examples/04_full_dataprocessing
Info: ata.diag not found: sipdir//invmod/02_100.000000//inv/ata.diag
/home/runner/work/crtomo_tools/crtomo_tools/examples/04_full_dataprocessing
reading voltages
Reading inversion results
is robust False
Info: res_m.diag not found: sipdir//invmod/03_1000.000000//inv/res_m.diag
/home/runner/work/crtomo_tools/crtomo_tools/examples/04_full_dataprocessing
Info: ata.diag not found: sipdir//invmod/03_1000.000000//inv/ata.diag
/home/runner/work/crtomo_tools/crtomo_tools/examples/04_full_dataprocessing
reading voltages
Reading inversion results
is robust False
Info: res_m.diag not found: sipdir//invmod/04_10000.000000//inv/res_m.diag
/home/runner/work/crtomo_tools/crtomo_tools/examples/04_full_dataprocessing
Info: ata.diag not found: sipdir//invmod/04_10000.000000//inv/ata.diag
/home/runner/work/crtomo_tools/crtomo_tools/examples/04_full_dataprocessing