Merge pull request #40 from ajithabhks/allow_more_freqencies

Allow more freqencies

Author: ajithabhks

README.md

@@ -1,15 +1,5 @@
 [![PyPI version](https://img.shields.io/pypi/v/sigmt)](https://pypi.org/project/sigmt/)
 
-Note: I’m currently in the process of updating module and variable names throughout the codebase to make them more descriptive and easier to understand.
-
-These updates may temporarily affect stability or introduce minor breaking changes.
-
-If you rely on a stable version for processing or production use, please use the latest release available on PyPI (`pip install sigmt`) or ([GitHub v2.0.4](https://github.com/ajithabhks/sigmt/tree/v2.0.4)) until the refactoring is complete.
-
-Sorry for any inconvenience.
-
----
-
 # SigMT: An open-source python package for magnetotelluric data processing
 
 SigMT is a python package designed for the processing of the raw magnetotelluric (MT) data to obtain the MT impedance and tipper estimates. It works in an automated way, so that manual time series inspection and editing are not required. Mahalanobis based data selection tool is implemented in the package to avoid the manual editing of time series. The final impedance estimation is done using the robust estimation method. Different data selection tools such as coherency threshold, polarization direction are included in this package.

sigmt/__version__.py

@@ -1,4 +1,4 @@
 """
 Returns version of the package
 """
-__version__ = "2.0.4"
+__version__ = "2.0.5"

sigmt/core/band_averaging.py

@@ -217,7 +217,7 @@ def __init__(self,
         self.detrend_time_series()
         self.perform_fft()
         if calibrate_magnetic:
-            self.calibrate_mag()
+            self.calibrate_magnetic()
         del self.calibration_data_magnetic
         del self.time_series
         gc.collect()
@@ -250,7 +250,7 @@ def calibrate_electric(self) -> None:
             dipole_ew = abs(self.calibration_data_electric['ey']['y1']) + abs(self.calibration_data_electric['ey']['y2'])
             self.time_series['ey'] = self.time_series['ey'] / (1 * dipole_ew / 1000)
 
-    def calibrate_mag(self) -> None:
+    def calibrate_magnetic(self) -> None:
         """
         Calibrates the magnetic field channels.
 

sigmt/core/plots.py

@@ -7,6 +7,7 @@
 import xarray as xr
 from matplotlib import pyplot as plt
 
+matplotlib.use("Qt5Agg")
 plt.rcParams['figure.max_open_warning'] = 50
 
 
@@ -55,8 +56,8 @@ def plot_mt_app_res(dataset: xr.Dataset, procinfo: dict) -> None:
     plt.errorbar(ftlist, rho_yx, yerr=err_ryx, ecolor='b', fmt="none")
     plt.xscale('log')
     plt.yscale('log')
-    if max(ftlist) < 15000 and min(ftlist) > 0.001:
-        plt.xlim((15000, 0.001))
+    if max(ftlist) < 15000 and min(ftlist) > 0.0001:
+        plt.xlim((15000, 0.0001))
     else:
         plt.xlim((max(ftlist) + 10, min(ftlist) - 10))
         if min(ftlist) > 0.001:
@@ -75,8 +76,8 @@ def plot_mt_app_res(dataset: xr.Dataset, procinfo: dict) -> None:
     plt.errorbar(ftlist, phase_xy, yerr=err_pxy, ecolor='r', fmt="none")
     plt.errorbar(ftlist, phase_yx, yerr=err_pyx, ecolor='b', fmt="none")
     plt.xscale('log')
-    if max(ftlist) < 15000 and min(ftlist) > 0.001:
-        plt.xlim((15000, 0.001))
+    if max(ftlist) < 15000 and min(ftlist) > 0.0001:
+        plt.xlim((15000, 0.0001))
     else:
         plt.xlim((max(ftlist) + 10, min(ftlist) - 10))
         if min(ftlist) > 0.001:
@@ -125,8 +126,8 @@ def plot_tipper(dataset: xr.Dataset, procinfo: dict) -> None:
     plt.scatter(ftlist, ty_a, c='b', s=10, label='Ty')
     plt.ylim(0, 1)
     plt.xscale('log')
-    if max(ftlist) < 15000 and min(ftlist) > 0.001:
-        plt.xlim((15000, 0.001))
+    if max(ftlist) < 15000 and min(ftlist) > 0.0001:
+        plt.xlim((15000, 0.0001))
     else:
         plt.xlim((max(ftlist) + 10, min(ftlist) - 10))
         if min(ftlist) > 0.001:
@@ -139,8 +140,8 @@ def plot_tipper(dataset: xr.Dataset, procinfo: dict) -> None:
     plt.scatter(ftlist, tx_p, c='r', s=10)
     plt.scatter(ftlist, ty_p, c='b', s=10)
     plt.xscale('log')
-    if max(ftlist) < 15000 and min(ftlist) > 0.001:
-        plt.xlim((15000, 0.001))
+    if max(ftlist) < 15000 and min(ftlist) > 0.0001:
+        plt.xlim((15000, 0.0001))
     else:
         plt.xlim((max(ftlist) + 10, min(ftlist) - 10))
         if min(ftlist) > 0.001:
@@ -175,7 +176,7 @@ def plot_coherency(dataset: xr.Dataset, procinfo: dict) -> None:
     if 'coh_hz' in dataset:
         plt.scatter(frequency_values, dataset['coh_hz'].values, c='g', label='Hz')
     plt.xscale('log')
-    plt.xlim((10000, 0.001))
+    plt.xlim((10000, 0.0001))
     plt.ylim(0, 1)
     # generate tick values for y
     y_ticks = np.arange(0.1, 1.1, 0.1)

sigmt/utils/utils.py

@@ -187,23 +187,27 @@ def get_target_frequency_list(sampling_frequency: float,
 
     fr = parzen_window_radius * frequency_table  # bandwidth of parzen window - oneside from ft
     total_bandwidth = fr * 2  # bandwidth of parzen window - two side from ft
-    # nof spectra in parzen window for each ft
+
+    # No. of spectral bins in the parzen window for each ft
     dof = total_bandwidth / (sampling_frequency / fft_length)
 
     maximum = sampling_frequency / 2
+
+    # Hardcoding maximum frequency to 15000
     if maximum > 15000:
         maximum = 15000
     f_max = max(frequency_table[frequency_table < maximum])  # ft_max, following nyquist
-
     f_max_index = np.where(frequency_table == f_max)[0][0]  # index in ftable
-    min_required = dof[f_max_index] * .01  # min dof required for ft_min
 
-    dof_min = min(dof[dof > min_required])  # finding in dof
-    dof_min_index = np.where(dof == dof_min)[0][0]  # finding index
+    # Finding minimum frequency bin used for a parzen window used.
+    parzen_lower_extend = frequency_table - (frequency_table * parzen_window_radius)
+    # Find the index of the target frequency whose lower Parzen window boundary is
+    # closest to the frequency resolution
+    f_min_index = np.argmin(abs(parzen_lower_extend - (sampling_frequency / fft_length)))
 
-    while dof[dof_min_index] < 10:  # Making sure, atleast 10 spectral lines are averaged
-        dof_min_index = dof_min_index - 1
+    # Looks for minimum dof=10
+    f_min_index = np.min([f_min_index, np.argmin(abs(dof-10))])
 
-    ft_list = frequency_table[f_max_index:dof_min_index + 1]  # Making ft_list
+    ft_list = frequency_table[f_max_index:f_min_index]  # Making ft_list
 
     return ft_list