Multiple spectra

It is often beneficial to compare the difference between various spectra calculation parameters. To make this easier, resistics supports the calculation of multiple spectra for each timeseries data folder. This is achieved through the specdir (short for spectra directory) option.

Note

This section describes the specdir option and how multiple spectra are saved and supported. However, the easiest way to specify the specdir option is in configuration files. For more information, please see Using configuration files.

By default, the specdir is called spectra and spectra data is saved in the following location:

exampleProject
├── calData
├── timeData
│   └── site1
|       |── dataFolder1
│       |── dataFolder2
|       |──     .
|       |──     .
|       |──     .
|       └── dataFolderN
├── specData
│   └── site1
|       |── dataFolder1
|       |   └── spectra
│       |── dataFolder2
|       |   └── spectra
|       |──     .
|       |──     .
|       |──     .
|       └── dataFolderN
|           └── spectra
├── statData
├── maskData
├── transFuncData
├── images
└── mtProj.prj

However, by specifying the specdir option in calculateSpectra(), a new set of spectra can be calculated. In the Viewing spectra section, it was clear that the time series spectra had a peak at 50 Hz, which is a common powerline frequency. Therefore, an example is shown below where a second set of spectra are calculated with a notch filter applied at 50 Hz and specdir is specified as specdir = “notch”.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
from datapaths import projectPath, imagePath
from resistics.project.io import loadProject

projData = loadProject(projectPath)

# calculate another set of spectra for the 128 Hz data with notching at 50Hz and 16.667Hz
from resistics.project.spectra import calculateSpectra

calculateSpectra(projData, sampleFreqs=[128], notch=[50], specdir="notch")
projData.refresh()

The new set of spectra data with specdir = “notch” are saved in the following way:

exampleProject
├── calData
├── timeData
│   └── site1
|       |── dataFolder1
│       |── dataFolder2
|       |──     .
|       |──     .
|       |──     .
|       └── dataFolderN
├── specData
│   └── site1
|       |── dataFolder1
|       |   |── notch
|       |   └── spectra
|       |
│       |── dataFolder2
|       |   |── notch
|       |   └── spectra
|       |──     .
|       |──     .
|       |──     .
|       └── dataFolderN
|           |── notch
|           └── spectra
├── statData
├── maskData
├── transFuncData
├── images
└── mtProj.prj

To view the new spectra data, specdir needs to be specified in the calls to viewSpectra(), viewSpectraSection() and viewSpectraStack(). An example is provided below using viewSpectraSection().

12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
# view the spectra
from resistics.common.plot import plotOptionsSpec, getPaperFonts
from resistics.project.spectra import viewSpectra, viewSpectraSection

plotOptions = plotOptionsSpec(plotfonts=getPaperFonts())
fig = viewSpectra(
    projData,
    "site1",
    "meas_2012-02-10_11-30-00",
    specdir="notch",
    plotoptions=plotOptions,
    show=False,
    save=False,
)
fig.savefig(imagePath / "multspec_viewspec_notch_spec")

In the plots below, the default spectra data and the notched spectra data are shown.

alternate text

Default spectra data

alternate text

Notched spectra data

A spectra section can also be plotted in a similar to way the example shown in Viewing Spectra.

28
29
30
31
32
33
34
35
36
37
fig = viewSpectraSection(
    projData,
    "site1",
    "meas_2012-02-10_11-30-00",
    specdir="notch",
    plotoptions=plotOptions,
    show=False,
    save=False,
)
fig.savefig(imagePath / "multspec_viewspec_notch_section")
alternate text

Notched spectra section

For the sake of comarpison, the default (no notch) spectra section.

alternate text

Default spectra section

Further, inspecting the comments of the notched spectra data, the application of the notch filter and the frequency is recorded.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
Unscaled data 2012-02-10 11:30:00 to 2012-02-11 23:03:43.992188 read in from measurement E:\magnetotellurics\code\resisticsdata\tutorial\tutorialProject\timeData\site1\meas_2012-02-10_11-30-00, samples 0 to 16387071
Sampling frequency 128.0
Scaling channel Ex with scalar -0.000112591 to give mV
Dividing channel Ex by electrode distance 0.08 km to give mV/km
Scaling channel Ey with scalar -0.000112512 to give mV
Dividing channel Ey by electrode distance 0.091 km to give mV/km
Scaling channel Hx with scalar -0.00360174 to give mV
Scaling channel Hy with scalar -0.00360038 to give mV
Scaling channel Hz with scalar -0.0036007 to give mV
Remove zeros: False, remove nans: False, remove average: True
---------------------------------------------------
Calculating project spectra
Using default configuration
Channel Ex not calibrated
Channel Ey not calibrated
Channel Hx calibrated with calibration data from file E:\magnetotellurics\code\resisticsdata\tutorial\tutorialProject\calData\Hx_MFS06365.TXT
Channel Hy calibrated with calibration data from file E:\magnetotellurics\code\resisticsdata\tutorial\tutorialProject\calData\Hy_MFS06357.TXT
Channel Hz calibrated with calibration data from file E:\magnetotellurics\code\resisticsdata\tutorial\tutorialProject\calData\Hz_MFS06307.TXT
Notch filter applied at 50 Hz
Decimating with 7 levels and 7 frequencies per level
Evaluation frequencies for this level 32.0, 22.62741699796952, 16.0, 11.31370849898476, 8.0, 5.65685424949238, 4.0
Windowing with window size 512 samples and overlap 128 samples
Time data decimated from 128.0 Hz to 16.0 Hz, new start time 2012-02-10 11:30:00, new end time 2012-02-11 23:03:43.937500
Evaluation frequencies for this level 2.82842712474619, 2.0, 1.414213562373095, 1.0, 0.7071067811865475, 0.5, 0.35355339059327373
Windowing with window size 512 samples and overlap 128 samples
Time data decimated from 16.0 Hz to 2.0 Hz, new start time 2012-02-10 11:30:00, new end time 2012-02-11 23:03:43.500000
Evaluation frequencies for this level 0.25, 0.17677669529663687, 0.125, 0.08838834764831843, 0.0625, 0.044194173824159216, 0.03125
Windowing with window size 512 samples and overlap 128 samples
Time data decimated from 2.0 Hz to 0.25 Hz, new start time 2012-02-10 11:30:00, new end time 2012-02-11 23:03:40
Time data decimated from 0.25 Hz to 0.125 Hz, new start time 2012-02-10 11:30:00, new end time 2012-02-11 23:03:36
Evaluation frequencies for this level 0.022097086912079608, 0.015625, 0.011048543456039804, 0.0078125, 0.005524271728019902, 0.00390625, 0.002762135864009951
Windowing with window size 512 samples and overlap 128 samples
Time data decimated from 0.125 Hz to 0.015625 Hz, new start time 2012-02-10 11:30:00, new end time 2012-02-11 23:03:20
Evaluation frequencies for this level 0.001953125, 0.0013810679320049755, 0.0009765625, 0.0006905339660024878, 0.00048828125, 0.0003452669830012439, 0.000244140625
Windowing with window size 512 samples and overlap 128 samples
Spectra data written out to E:\magnetotellurics\code\resisticsdata\tutorial\tutorialProject\specData\site1\meas_2012-02-10_11-30-00\notch on 2019-10-05 20:10:06.845930 using resistics 0.0.6.dev2
---------------------------------------------------

The next step is to process the notched spectra data. To do this, specdir has to be specified in the call to processProject() as shown below:

39
40
41
42
# process the new set of spectra
from resistics.project.transfunc import processProject

processProject(projData, sites=["site1"], specdir="notch")

Finally, the transfer function data for the notched spectra data can be viewed by once more specifying the specdir option in the call to viewImpedance()

44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
# plot the transfer functions, again with specifying the relevant specdir
from resistics.project.transfunc import viewImpedance

figs = viewImpedance(
    projData,
    sites=["site1"],
    sampleFreqs=[128],
    oneplot=False,
    specdir="notch",
    save=False,
    show=False,
)
figs[0].savefig(imagePath / "multspec_viewimp_notch")

# and compare to the original
figs = viewImpedance(
    projData,
    sites=["site1"],
    sampleFreqs=[128],
    oneplot=False,
    specdir="spectra",
    save=False,
    show=False,
)
figs[0].savefig(imagePath / "multspec_viewimp_nonotch")

The impedance for the default spectra calculated parameters and then for the notch parameters are shown below.

alternate text

Default impedance tensor estimation

alternate text

Impedance tensor estimation using the notched spectra data

As mentioned in the note at the top of this page, the simplest way to use the specdir option is through configuration files.

Important

When specifying specdir in a configuration file, the specdir keyword no longer has to be passed to the various method calls as it is picked up directly from the confiugration file. This makes everything much simpler and removes a source of mistakes.

More information about configuration files is provided in the Configuration files section. An example of using configuration files is provided in the Using configuration files section.

Complete example script

For the purposes of clarity, the complete example script is provided below.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
from datapaths import projectPath, imagePath
from resistics.project.io import loadProject

projData = loadProject(projectPath)

# calculate another set of spectra for the 128 Hz data with notching at 50Hz and 16.667Hz
from resistics.project.spectra import calculateSpectra

calculateSpectra(projData, sampleFreqs=[128], notch=[50], specdir="notch")
projData.refresh()

# view the spectra
from resistics.common.plot import plotOptionsSpec, getPaperFonts
from resistics.project.spectra import viewSpectra, viewSpectraSection

plotOptions = plotOptionsSpec(plotfonts=getPaperFonts())
fig = viewSpectra(
    projData,
    "site1",
    "meas_2012-02-10_11-30-00",
    specdir="notch",
    plotoptions=plotOptions,
    show=False,
    save=False,
)
fig.savefig(imagePath / "multspec_viewspec_notch_spec")

fig = viewSpectraSection(
    projData,
    "site1",
    "meas_2012-02-10_11-30-00",
    specdir="notch",
    plotoptions=plotOptions,
    show=False,
    save=False,
)
fig.savefig(imagePath / "multspec_viewspec_notch_section")

# process the new set of spectra
from resistics.project.transfunc import processProject

processProject(projData, sites=["site1"], specdir="notch")

# plot the transfer functions, again with specifying the relevant specdir
from resistics.project.transfunc import viewImpedance

figs = viewImpedance(
    projData,
    sites=["site1"],
    sampleFreqs=[128],
    oneplot=False,
    specdir="notch",
    save=False,
    show=False,
)
figs[0].savefig(imagePath / "multspec_viewimp_notch")

# and compare to the original
figs = viewImpedance(
    projData,
    sites=["site1"],
    sampleFreqs=[128],
    oneplot=False,
    specdir="spectra",
    save=False,
    show=False,
)
figs[0].savefig(imagePath / "multspec_viewimp_nonotch")