Remote reference statistics

Remote reference statistics are a means of understanding how adding a remote reference site to the processing of magnetotelluric data is changing the solution. There are currently five remote reference statistics.

  • “RR_coherence”

  • “RR_coherenceEqn”

  • “RR_absvalEqn”

  • “RR_transferFunction”

  • “RR_resPhase”

More information on each of these is given in remote reference statistics.

Important

Remote reference statistics are always saved with the local site.

Calculating and viewing remote reference statistics

Begin by loading the project and calculate the remote statistics using the calculateRemoteStatistics() of module statistics. The calculateRemoteStatistics() method will calculate the remote reference statistics for sites=[“M6”] and at sampling frequencies sampleFreqs=[128] using the remote site “Remote”. In this instance, all the remote statistics are being calculated.

 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
from datapaths import remotePath, remoteImages
from resistics.project.io import loadProject
from resistics.project.statistics import (
    calculateRemoteStatistics,
    viewStatistic,
    viewStatisticHistogram,
    viewStatisticDensityplot,
)
from resistics.project.mask import newMaskData, calculateMask
from resistics.project.transfunc import processProject, viewImpedance
from resistics.common.plot import plotOptionsStandard, getPaperFonts

plotOptions = plotOptionsStandard(plotfonts=getPaperFonts())
proj = loadProject(remotePath)

calculateRemoteStatistics(
    proj,
    "Remote",
    sites=["M6"],
    sampleFreqs=[128],
    remotestats=[
        "RR_coherence",
        "RR_coherenceEqn",
        "RR_absvalEqn",
        "RR_transferFunction",
        "RR_resPhase",
    ],
)

To get an idea how the remote reference is changing the impedance tensor values, the impedance tensor estimates for each window for the single site and remote reference processing can now be plotted. To do this, the viewStatisticDensityplot() of module statistics is used. This method plots 2-D histograms of statistic data cross pairs. In the following code snippet, density plots are being saved for the first evaluation frequency of decimation levels 0, 1, 2 and 3. Limits for x and y axes are also being set explicitly.

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
lims = {0: 200, 1: 120, 2: 50, 3: 30}
for declevel in range(0, 4):
    fig = viewStatisticDensityplot(
        proj,
        "M6",
        128,
        "RR_transferFunction",
        declevel=declevel,
        crossplots=[["ExHyRealRR", "ExHyImagRR"], ["EyHxRealRR", "EyHxImagRR"]],
        xlim=[-lims[declevel], lims[declevel]],
        ylim=[-lims[declevel], lims[declevel]],
        plotoptions=plotOptions,
        show=False,
    )
    fig.savefig(remoteImages / "densityPlot_remoteRef_{}.png".format(declevel))
    fig = viewStatisticDensityplot(
        proj,
        "M6",
        128,
        "transferFunction",
        declevel=declevel,
        crossplots=[["ExHyReal", "ExHyImag"], ["EyHxReal", "EyHxImag"]],
        xlim=[-lims[declevel], lims[declevel]],
        ylim=[-lims[declevel], lims[declevel]],
        plotoptions=plotOptions,
        show=False,
    )
    fig.savefig(remoteImages / "densityPlot_singleSite_{}.png".format(declevel))

For single site statistics, the window-by-window impedance tensor estimates are calculated in the “transferFunction” statistic. This has components: ExHxReal, ExHxImag, ExHyReal, ExHyImag, EyHxReal, EyHxImag, EyHyReal, EyHyImag. For the density plot, the cross pairs ExHyReal vs. ExHyImag and EyHxReal vs. EyHxImag are plotted. This is essentially the Ex-Hy and Ey-Hx polarisations of the single site impedance tensor estimates plotted in the complex plane.

alternate text

Density plot of single site, window-by-window, impedance tensor Ex-Hy and Ey-Hx estimates for decimation level 0, evaluation frequency 32 Hz or period 0.03125 seconds.

For remote reference statistics, the window-by-window impedance tensor estimates are calculated in the “RR_transferFunction” statistic. This has components: ExHxRealRR, ExHxImagRR, ExHyRealRR, ExHyImagRR, EyHxRealRR, EyHxImagRR, EyHyRealRR, EyHyImagRR. For the density plot, the cross pairs ExHyRealRR vs. ExHyImagRR and EyHxRealRR vs. EyHxImagRR are plotted. This is essentially the Ex-Hy and Ey-Hx polarisations of the remote reference impedance tensor estimates plotted in the complex plane.

alternate text

Density plot of remote reference, window-by-window, impedance tensor Ex-Hy and Ey-Hx estimates for decimation level 0, evaluation frequency 32 Hz or period 0.03125 seconds.

The same plots for decimation level 2 are shown below.

alternate text

Density plot of single site, window-by-window, impedance tensor Ex-Hy and Ey-Hx estimates for decimation level 2, evaluation frequency 0.25 Hz or period 4 seconds.

alternate text

Density plot of remote reference, window-by-window, impedance tensor Ex-Hy and Ey-Hx estimates for decimation level 2, evaluation frequency 0.25 Hz or period 4 seconds.

The plots demonstrate how adding a remote reference is having an impact on the impedance tensor estimates. Particularly the plots for decimation level 2 (period 4 seconds) in the dead band show that by adding the remote reference, not only are the impedance tensor estimates changing, but the spread is becoming larger. Below are the single site impedance tensor estimates for Site M6 and the standard remote reference imepdance tensor estimate for Site M6. Visible in the impedance tensor plots is that the single site phase is smooth (though probably biased) in the dead band, whilst for the same periods, the remote reference estimates show significantly more noise. This is due to the greater spread in window-by-window impedance tensor estimates apparent in the remote reference processing.

alternate text

Single site impedance tensor estimates for Site M6 at 128 Hz

alternate text

Standard remote reference impedance tensor estimate for Site M6.

Masking with remote reference statistics

Now, the remote reference statistics can be used to generate a mask. The easiest one to use is the remote reference equivalent of coherence. The statistic for this is “RR_coherenceEqn”. Below, a mask is generated for “RR_coherenceEqn” requiring a minimum coherence of 0.8 for the pairs “ExHyR-HyHyR” and “EyHxR-HxHxR”.

59
60
61
62
63
64
65
66
67
68
69
# create a mask that uses some of the remote reference statistics that were calculated
# get a mask data object and specify the sampling frequency to mask (128Hz)
maskData = newMaskData(proj, 128)
maskData.setStats(["RR_coherenceEqn"])
maskData.addConstraint(
    "RR_coherenceEqn", {"ExHyR-HyHyR": [0.8, 1.0], "EyHxR-HxHxR": [0.8, 1.0]}
)
# finally, lets give maskData a name, which will relate to the output file
maskData.maskName = "rr_cohEqn_80_100"
calculateMask(proj, maskData, sites=["M6"])
maskData.printInfo()

Important

Because remote reference statistics are always associated with the local site, any masks generated from them are also saved to the local site.

Now, the Site M6 can be remote reference processed using the new mask and the impedance tensor estimate plotted.

71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
# process with masks
processProject(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    remotesite="Remote",
    masks={"M6": "rr_cohEqn_80_100"},
    postpend="rr_cohEqn_80_100",
)

from resistics.common.plot import plotOptionsTransferFunction

plotOptionsTF = plotOptionsTransferFunction(plotfonts=getPaperFonts())
figs = viewImpedance(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    postpend="rr_cohEqn_80_100",
    oneplot=False,
    plotoptions=plotOptionsTF,
    save=False,
    show=False,
)
figs[0].savefig(remoteImages / "remoteReferenceM6_128_RR_coh.png")
alternate text

Remote reference impedance tensor estimate for Site M6 using a mask created from remote reference statistics.

It is also possible to see how the mask has affected the density plot by providing a mask into the viewStatisticDensityplot() of module statistics. An example is shown below.

 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
# see how the masks changed the density plot
lims = {0: 200, 1: 120, 2: 50, 3: 30}
for declevel in range(0, 4):
    fig = viewStatisticDensityplot(
        proj,
        "M6",
        128,
        "RR_transferFunction",
        declevel=declevel,
        crossplots=[["ExHyRealRR", "ExHyImagRR"], ["EyHxRealRR", "EyHxImagRR"]],
        maskname="rr_cohEqn_80_100",
        xlim=[-lims[declevel], lims[declevel]],
        ylim=[-lims[declevel], lims[declevel]],
        plotoptions=plotOptions,
        show=False,
    )
    fig.savefig(remoteImages / "densityPlot_remoteRef_{}_mask.png".format(declevel))
alternate text

The effect of masking on the density plot for remote reference impedance tensor estimates for Site M6 and decimation level 2, evaluation frequency 0.25 Hz and period 4 seconds.

Masks and date or time constraints

Recall, however, that masks can be combined with other masks as well as with datetime constraints. In the following example, the remote reference processing is performed again, but with the addition of a datetime constraint which includes only data recorded during the night.

114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
# process with masks and datetime constraints
processProject(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    remotesite="Remote",
    masks={"M6": "rr_cohEqn_80_100"},
    datetimes=[{"type": "time", "start": "20:00:00", "stop": "07:00:00"}],
    postpend="rr_cohEqn_80_100_night",
)

figs = viewImpedance(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    postpend="rr_cohEqn_80_100_night",
    oneplot=False,
    plotoptions=plotOptionsTF,    
    save=False,
    show=False,
)
figs[0].savefig(remoteImages / "remoteReferenceM6_128_RR_coh_datetime.png")
alternate text

Remote reference impedance tensor estimate for Site M6 using a mask created from remote reference statistics and with additional time constraints.

The new impedance tensor esimate appears better at the low decimation levels where there are plently of windows and worse at the high decimation levels where it is reducing the number of windows which can be included in the processing. The remote reference processing can be performed again, but this time applying the time constraints only to decimation levels 0 and 1.

137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
# process with masks and datetime constraints, but only for decimation levels 0 and 1
processProject(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    remotesite="Remote",
    masks={"M6": "rr_cohEqn_80_100"},
    datetimes=[
        {"type": "time", "start": "20:00:00", "stop": "07:00:00", "levels": [0, 1]}
    ],
    postpend="rr_cohEqn_80_100_night2",
)

figs = viewImpedance(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    postpend="rr_cohEqn_80_100_night2",
    oneplot=False,
    plotoptions=plotOptionsTF,    
    save=False,
    show=False,
)
figs[0].savefig(remoteImages / "remoteReferenceM6_128_RR_coh_datetime_01.png")
alternate text

Remote reference impedance tensor estimate for Site M6 using a mask created from remote reference statistics and with additional time constraints only applied to decimation levels 0 and 1.

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
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
from datapaths import remotePath, remoteImages
from resistics.project.io import loadProject
from resistics.project.statistics import (
    calculateRemoteStatistics,
    viewStatistic,
    viewStatisticHistogram,
    viewStatisticDensityplot,
)
from resistics.project.mask import newMaskData, calculateMask
from resistics.project.transfunc import processProject, viewImpedance
from resistics.common.plot import plotOptionsStandard, getPaperFonts

plotOptions = plotOptionsStandard(plotfonts=getPaperFonts())
proj = loadProject(remotePath)

calculateRemoteStatistics(
    proj,
    "Remote",
    sites=["M6"],
    sampleFreqs=[128],
    remotestats=[
        "RR_coherence",
        "RR_coherenceEqn",
        "RR_absvalEqn",
        "RR_transferFunction",
        "RR_resPhase",
    ],
)

lims = {0: 200, 1: 120, 2: 50, 3: 30}
for declevel in range(0, 4):
    fig = viewStatisticDensityplot(
        proj,
        "M6",
        128,
        "RR_transferFunction",
        declevel=declevel,
        crossplots=[["ExHyRealRR", "ExHyImagRR"], ["EyHxRealRR", "EyHxImagRR"]],
        xlim=[-lims[declevel], lims[declevel]],
        ylim=[-lims[declevel], lims[declevel]],
        plotoptions=plotOptions,
        show=False,
    )
    fig.savefig(remoteImages / "densityPlot_remoteRef_{}.png".format(declevel))
    fig = viewStatisticDensityplot(
        proj,
        "M6",
        128,
        "transferFunction",
        declevel=declevel,
        crossplots=[["ExHyReal", "ExHyImag"], ["EyHxReal", "EyHxImag"]],
        xlim=[-lims[declevel], lims[declevel]],
        ylim=[-lims[declevel], lims[declevel]],
        plotoptions=plotOptions,
        show=False,
    )
    fig.savefig(remoteImages / "densityPlot_singleSite_{}.png".format(declevel))

# create a mask that uses some of the remote reference statistics that were calculated
# get a mask data object and specify the sampling frequency to mask (128Hz)
maskData = newMaskData(proj, 128)
maskData.setStats(["RR_coherenceEqn"])
maskData.addConstraint(
    "RR_coherenceEqn", {"ExHyR-HyHyR": [0.8, 1.0], "EyHxR-HxHxR": [0.8, 1.0]}
)
# finally, lets give maskData a name, which will relate to the output file
maskData.maskName = "rr_cohEqn_80_100"
calculateMask(proj, maskData, sites=["M6"])
maskData.printInfo()

# process with masks
processProject(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    remotesite="Remote",
    masks={"M6": "rr_cohEqn_80_100"},
    postpend="rr_cohEqn_80_100",
)

from resistics.common.plot import plotOptionsTransferFunction

plotOptionsTF = plotOptionsTransferFunction(plotfonts=getPaperFonts())
figs = viewImpedance(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    postpend="rr_cohEqn_80_100",
    oneplot=False,
    plotoptions=plotOptionsTF,
    save=False,
    show=False,
)
figs[0].savefig(remoteImages / "remoteReferenceM6_128_RR_coh.png")

# see how the masks changed the density plot
lims = {0: 200, 1: 120, 2: 50, 3: 30}
for declevel in range(0, 4):
    fig = viewStatisticDensityplot(
        proj,
        "M6",
        128,
        "RR_transferFunction",
        declevel=declevel,
        crossplots=[["ExHyRealRR", "ExHyImagRR"], ["EyHxRealRR", "EyHxImagRR"]],
        maskname="rr_cohEqn_80_100",
        xlim=[-lims[declevel], lims[declevel]],
        ylim=[-lims[declevel], lims[declevel]],
        plotoptions=plotOptions,
        show=False,
    )
    fig.savefig(remoteImages / "densityPlot_remoteRef_{}_mask.png".format(declevel))

# process with masks and datetime constraints
processProject(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    remotesite="Remote",
    masks={"M6": "rr_cohEqn_80_100"},
    datetimes=[{"type": "time", "start": "20:00:00", "stop": "07:00:00"}],
    postpend="rr_cohEqn_80_100_night",
)

figs = viewImpedance(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    postpend="rr_cohEqn_80_100_night",
    oneplot=False,
    plotoptions=plotOptionsTF,    
    save=False,
    show=False,
)
figs[0].savefig(remoteImages / "remoteReferenceM6_128_RR_coh_datetime.png")

# process with masks and datetime constraints, but only for decimation levels 0 and 1
processProject(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    remotesite="Remote",
    masks={"M6": "rr_cohEqn_80_100"},
    datetimes=[
        {"type": "time", "start": "20:00:00", "stop": "07:00:00", "levels": [0, 1]}
    ],
    postpend="rr_cohEqn_80_100_night2",
)

figs = viewImpedance(
    proj,
    sites=["M6"],
    sampleFreqs=[128],
    postpend="rr_cohEqn_80_100_night2",
    oneplot=False,
    plotoptions=plotOptionsTF,    
    save=False,
    show=False,
)
figs[0].savefig(remoteImages / "remoteReferenceM6_128_RR_coh_datetime_01.png")