Merge pull request #100 from MTgeophysics/updates

Updates to Simpeg inversions

Author: kujaku11

.github/workflows/testing_in_conda.yml

@@ -96,6 +96,15 @@ jobs:
         python -c "import pandas as pd; print(f'pandas version: {pd.__version__}')"
         conda list
 
+    - name: Run Simpeg1D tests (non-Windows)
+      if: runner.os != 'Windows'
+      shell: bash
+      run: |
+        source ~/.profile
+        conda init --all
+        conda activate mtpy-v2-test
+        pytest -rA -q tests/modeling/simpeg/test_simpeg_1d_inversion_recipe.py
+
     - name: Run Tests
       shell: bash
       run: |

README.md

@@ -7,7 +7,7 @@
 [![Documentation Status](https://readthedocs.org/projects/mtpy-v2/badge/?version=latest)](https://mtpy-v2.readthedocs.io/en/latest/?badge=latest)
 [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/MTgeophysics/mtpy-v2/main)
 
-## Version 2.1.1
+## Version 2.1.2
 
 # Description
 

docs/source/conf.py

@@ -67,9 +67,9 @@
 #
 
 # The short X.Y version.
-version = "2.1.1"
+version = "2.1.2"
 # The full version, including alpha/beta/rc tags.
-release = "2.1.1"
+release = "2.1.2"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

mtpy/__init__.py

@@ -18,7 +18,7 @@
 from mtpy.imaging.mtcolors import MT_CMAP_DICT, register_cmaps
 
 
-__version__ = "2.1.1"
+__version__ = "2.1.2"
 __all__ = ["MT", "MTData", "MTCollection"]
 
 # =============================================================================

mtpy/core/mt.py

@@ -986,7 +986,8 @@ def to_dataframe(
         cols : list of str, optional
             Column names to include, by default None
         impedance_units : str, optional
-            Impedance units, "mt" [mV/km/nT] or "ohm" [Ohms], by default "mt"
+            Impedance output units in the dataframe,
+            "mt" [mV/km/nT] or "ohm" [Ohms], by default "mt"
 
         Returns
         -------
@@ -1017,8 +1018,10 @@ def to_dataframe(
         mt_df.dataframe.loc[:, "period"] = self.period
         if self.has_impedance():
             z_object = self.Z
-            z_object.units = impedance_units
-            mt_df.from_z_object(z_object)
+            # Keep Z internal storage in mt units while exposing dataframe values
+            # in requested output units.
+            z_object.output_units = impedance_units
+            mt_df.from_z_object(z_object, units=impedance_units)
         if self.has_tipper():
             mt_df.from_t_object(self.Tipper)
 
@@ -1035,7 +1038,8 @@ def from_dataframe(
         mt_df : MTDataFrame or DataFrame-like
             Dataframe containing MT data for a single station
         impedance_units : str, optional
-            Impedance units, "mt" [mV/km/nT] or "ohm" [Ohms], by default "mt"
+            Impedance units represented in the dataframe,
+            "mt" [mV/km/nT] or "ohm" [Ohms], by default "mt"
 
         Raises
         ------
@@ -1077,8 +1081,9 @@ def from_dataframe(
 
         self.tf_id = self.station
 
-        z_obj = mt_df.to_z_object()
-        z_obj.units = impedance_units
+        # Construct with dataframe units so incoming arrays are interpreted
+        # correctly, and keep output units aligned with MT impedance units.
+        z_obj = mt_df.to_z_object(units=impedance_units)
         self.Z = z_obj
         self.Tipper = mt_df.to_t_object()
 
@@ -1542,7 +1547,9 @@ def to_occam1d(
 
         return occam_data
 
-    def to_simpeg_1d(self, mode: str = "det", **kwargs: Any) -> Simpeg1D:
+    def to_simpeg_1d(
+        self, mode: str = "det", resistivity_error=10, phase_error=2.5, **kwargs: Any
+    ) -> Simpeg1D:
         """
         Run a 1D inversion using SimPEG.
 
@@ -1577,7 +1584,13 @@ def to_simpeg_1d(self, mode: str = "det", **kwargs: Any) -> Simpeg1D:
         if not self.Z._has_tf_model_error():
             self.compute_model_z_errors()
             self.logger.info("Using default errors for impedance")
-        simpeg_1d = Simpeg1D(self.to_dataframe(), mode=mode, **kwargs)
+        simpeg_1d = Simpeg1D(
+            self.to_dataframe(impedance_units="mt"),
+            mode=mode,
+            resistivity_error=resistivity_error,
+            phase_error=phase_error,
+            **kwargs,
+        )
         simpeg_1d.run_fixed_layer_inversion(**kwargs)
         simpeg_1d.plot_model_fitting(fig_num=1)
         simpeg_1d.plot_response(fig_num=2, **kwargs)

mtpy/core/mt_data.py

@@ -1252,6 +1252,12 @@ def rotate_dask(
     ) -> "MTData":
         """Rotate stations using dask-delayed execution.
 
+        **Note**: Rotation is defined as a rotation of the coordinate reference frame
+        of the transfer function, not a rotation of the physical measurement. So,
+        for example, a 90 degree rotation of a station with NED coordinates would
+        swap the North and East components of the transfer function and change the
+        coordinate reference frame to ENU.
+
         Parameters
         ----------
         rotation_angle : float or ndarray
@@ -1720,9 +1726,10 @@ def _station_dataset_to_dataframe(
                     z_error=tf_error,
                     frequency=1.0 / period,
                     z_model_error=tf_model_error,
-                    units=impedance_units,
+                    input_units="mt",
+                    output_units=impedance_units,
                 )
-                station_df.from_z_object(z_object)
+                station_df.from_z_object(z_object, units=impedance_units)
 
             t_output = self._pick_channel_labels(output_labels, ["hz", "z"], 1)
             t_inputs = self._pick_channel_labels(
@@ -3318,6 +3325,16 @@ def rotate(
     ) -> "MTData" | None:
         """Rotate all station transfer functions.
 
+        **Note**: This method is not intended for general coordinate
+        rotation of station locations.  It is designed to rotate the
+        impedance and tipper channels of each station dataset
+        according to the specified angle and coordinate reference frame.
+        The station location coordinates are not modified by this method.
+        Rotation is off the coordinate system therefore a positive
+        clockwise rotation of 10 degrees will rotate the coordinate
+        system 10 degrees and the estimated strike angle will be 10
+        degrees less than the strike angle before rotation.
+
         Parameters
         ----------
         rotation_angle : float or ndarray

mtpy/core/mt_dataframe.py

@@ -1049,34 +1049,39 @@ def _fill_data(
                     column,
                 ] = data_array[:, index["ii"], index["jj"]]
 
-    def from_z_object(self, z_object: Z) -> None:
+    def from_z_object(self, z_object: Z, units="mt") -> None:
         """
         Fill dataframe from a Z impedance object.
 
         Parameters
         ----------
         z_object : Z
             Z object containing impedance tensor data
+        units : str, optional
+            Units to write impedance values into the dataframe,
+            by default "mt"
 
         Notes
         -----
         Populates impedance, resistivity, phase, and phase tensor columns
 
         """
+        z_export = z_object.copy()
+        z_export.output_units = units
 
         self.dataframe.loc[self.dataframe.station == self.station, "period"] = (
-            z_object.period
+            z_export.period
         )
 
         # should make a copy of the phase tensor otherwise it gets calculated
         # multiple times and becomes a time sink.
-        pt_object = z_object.phase_tensor.copy()
+        pt_object = z_export.phase_tensor.copy()
 
         for error in ["", "_error", "_model_error"]:
-            if getattr(z_object, f"_has_tf{error}")():
+            if getattr(z_export, f"_has_tf{error}")():
                 for key in ["z", "res", "phase"]:
                     obj_key = self._key_dict[key]
-                    data_array = getattr(z_object, f"{obj_key}{error}").copy()
+                    data_array = getattr(z_export, f"{obj_key}{error}").copy()
                     for comp in ["xx", "xy", "yx", "yy"]:
                         index = self._get_index(comp)
                         # if key in ["pt"]:
@@ -1150,7 +1155,9 @@ def to_z_object(self, units: str = "mt") -> Z:
         Parameters
         ----------
         units : str, optional
-            Impedance units ('mt' or 'ohm'), by default 'mt'
+            Units represented by impedance values in the dataframe
+            and returned by the resulting Z object ('mt' or 'ohm'),
+            by default 'mt'
 
         Returns
         -------
@@ -1191,7 +1198,14 @@ def to_z_object(self, units: str = "mt") -> Z:
                 self.dataframe.station == self.station, f"z_{comp}_model_error"
             ]
 
-        z_object = Z(z, z_err, self.frequency, z_model_err, units=units)
+        z_object = Z(
+            z,
+            z_err,
+            self.frequency,
+            z_model_err,
+            input_units=units,
+            output_units=units,
+        )
 
         if (z == 0).all():
             for comp in ["xx", "xy", "yx", "yy"]:
@@ -1224,6 +1238,8 @@ def to_z_object(self, units: str = "mt") -> Z:
 
             if not (res == 0).all():
                 if not (phase == 0).all():
+                    # Resistivity/phase-to-Z computation returns mt-based impedance.
+                    z_object.input_units = "mt"
                     z_object.set_resistivity_phase(
                         res,
                         phase,

mtpy/core/transfer_function/base.py

@@ -610,6 +610,14 @@ def rotate(
             positive angle. In this coordinate system the rotation matrix is
             the inverse of the conventional rotation matrix.
 
+        **Note**: Rotation is a rotation of the coordinate reference
+        frame of the transfer function, not a rotation of the physical measurement.
+        So, for example, if you have a transfer function that is referenced to
+        a coordinate system that is rotated 30 degrees clockwise from north,
+        then you would rotate by -30 degrees to get the transfer function in
+        a coordinate system that is aligned with north. This is because the
+        rotation is of the reference frame, not the physical measurement.
+
         Parameters
         ----------
         alpha : float or int or str or array-like

mtpy/core/transfer_function/tf_helpers.py

@@ -6,19 +6,31 @@
 
 import mtpy.utils.calculator as MTcc
 
+MT_TO_OHM_FACTOR = 1.0 / np.pi * np.sqrt(5.0 / 8.0) * 10**3.5
+IMPEDANCE_UNITS = {"mt": 1, "ohm": MT_TO_OHM_FACTOR}
+
 
 def compute_resistivity(
-    z_mt: np.ndarray | None, frequency: np.ndarray
+    z_mt: np.ndarray | None, frequency: np.ndarray, z_units: str = "mt"
 ) -> np.ndarray | None:
     """Compute apparent resistivity from impedance in mt units."""
     if z_mt is None:
         return None
 
-    return np.apply_along_axis(
-        lambda values: np.abs(values) ** 2 / frequency * 0.2,
-        0,
-        z_mt,
-    )
+    if z_units == "mt":
+        return np.apply_along_axis(
+            lambda values: np.abs(values) ** 2 / frequency * 0.2,
+            0,
+            z_mt,
+        )
+    elif z_units == "ohm":
+        return np.apply_along_axis(
+            lambda values: np.abs(values / IMPEDANCE_UNITS["ohm"]) ** 2
+            * frequency
+            * 0.2,
+            0,
+            z_mt,
+        )
 
 
 def compute_phase(z_mt: np.ndarray | None) -> np.ndarray | None:
@@ -33,29 +45,43 @@ def compute_resistivity_error(
     z_mt: np.ndarray | None,
     z_error_mt: np.ndarray | None,
     frequency: np.ndarray,
+    z_units: str = "mt",
 ) -> np.ndarray | None:
     """Compute apparent resistivity error from impedance and standard deviation."""
     if z_mt is None or z_error_mt is None:
         return None
 
     with np.errstate(divide="ignore", invalid="ignore"):
-        return np.apply_along_axis(
-            lambda values: values / frequency * 0.2,
-            0,
-            2 * z_error_mt * np.abs(z_mt),
-        )
+        if z_units == "mt":
+            return np.apply_along_axis(
+                lambda values: values / frequency * 0.2,
+                0,
+                2 * z_error_mt * np.abs(z_mt),
+            )
+        elif z_units == "ohm":
+            return np.apply_along_axis(
+                lambda values: values / frequency * 0.2,
+                0,
+                2 * z_error_mt * np.abs(z_mt / IMPEDANCE_UNITS["ohm"]),
+            )
 
 
 def compute_phase_error(
     z_mt: np.ndarray | None,
     z_error: np.ndarray | None,
+    z_units: str = "mt",
 ) -> np.ndarray | None:
     """Compute impedance phase error in degrees from complex impedance error."""
     if z_mt is None or z_error is None:
         return None
 
     with np.errstate(divide="ignore", invalid="ignore"):
-        return np.degrees(np.arctan(z_error / np.abs(z_mt)))
+        if z_units == "mt":
+            return np.degrees(np.arctan(z_error / np.abs(z_mt)))
+        elif z_units == "ohm":
+            return np.degrees(
+                np.arctan(z_error / np.abs(z_mt / IMPEDANCE_UNITS["ohm"]))
+            )
 
 
 def compute_impedance_error(