Merge pull request #25 from MTgeophysics/updates

- Adding `to/from_json` to `MTLocation`
- Fix `plot_mt_response` bugs
- Changed keyword in `MTData.from_modem` from `file_type` to `survey`

Author: kujaku11

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 2.0.6
+current_version = 2.0.7
 files = setup.py mtpy/__init__.py README.md docs/source/conf.py
 commit = True
 tag = True

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.0.6
+## Version 2.0.7
 
 # Description
 

docs/source/conf.py

@@ -67,9 +67,9 @@
 #
 
 # The short X.Y version.
-version = "2.0.6"
+version = "2.0.7"
 # The full version, including alpha/beta/rc tags.
-release = "2.0.6"
+release = "2.0.7"
 
 # 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.0.6"
+__version__ = "2.0.7"
 __all__ = ["MT", "MTData", "MTCollection"]
 
 # =============================================================================

mtpy/core/mt_data.py

@@ -889,7 +889,7 @@ def to_modem(self, data_filename=None, **kwargs):
 
         return modem_data
 
-    def from_modem_data(self, data_filename, file_type="data", **kwargs):
+    def from_modem_data(self, data_filename, survey="data", **kwargs):
         """
 
         :param data_filename: DESCRIPTION
@@ -907,7 +907,7 @@ def from_modem_data(self, data_filename, file_type="data", **kwargs):
             "'from_modem_data' will be deprecated in future versions, use 'from_modem'"
         )
 
-        self.from_modem(data_filename, file_type=file_type, **kwargs)
+        self.from_modem(data_filename, survey=survey, **kwargs)
 
     def from_modem(self, data_filename, survey="data", **kwargs):
         """
@@ -923,6 +923,7 @@ def from_modem(self, data_filename, survey="data", **kwargs):
         """
         modem_data = Data(**kwargs)
         mdf = modem_data.read_data_file(data_filename)
+
         # set survey name to something useful
         mdf.dataframe["survey"] = survey
 

mtpy/core/mt_location.py

@@ -14,6 +14,7 @@
 from pyproj import CRS
 import numpy as np
 from loguru import logger
+import json
 
 from mtpy.utils.gis_tools import project_point
 
@@ -30,7 +31,6 @@ class MTLocation:
     """
 
     def __init__(self, survey_metadata=None, **kwargs):
-
         self.logger = logger
         if survey_metadata is None:
             self._survey_metadata = self._initiate_metadata()
@@ -256,11 +256,7 @@ def utm_crs(self, value):
         new_crs = CRS.from_user_input(value)
         if value != self._utm_crs:
             # reproject easting, northing to new zone
-            if (
-                self._utm_crs is not None
-                and self.east != 0
-                and self.north != 0
-            ):
+            if self._utm_crs is not None and self.east != 0 and self.north != 0:
                 self._east, self._north = project_point(
                     self.east, self.north, self._utm_crs, new_crs
                 )
@@ -483,3 +479,52 @@ def get_elevation_from_national_map(self):
             self.logger.warning(
                 "Could not get elevation data, not setting elevation"
             )
+
+    def to_json(self, filename):
+        """
+        write out information to a json file
+
+        :param filename: DESCRIPTION
+        :type filename: TYPE
+        :return: DESCRIPTION
+        :rtype: TYPE
+
+        """
+
+        js_dict = {
+            "latitude": self.latitude,
+            "longitude": self.longitude,
+            "elevation": self.elevation,
+            "north": self.north,
+            "east": self.east,
+            "model_east": self.model_east,
+            "model_north": self.model_north,
+            "model_elevation": self.model_elevation,
+            "datum_crs": self.datum_crs.to_json_dict(),
+            "utm_crs": self.utm_crs.to_json_dict(),
+        }
+        with open(filename, "w") as fid:
+            json.dump(js_dict, fid)
+
+    def from_json(self, filename):
+        """
+        read in json formatted location
+
+        :param filename: DESCRIPTION
+        :type filename: TYPE
+        :return: DESCRIPTION
+        :rtype: TYPE
+
+        """
+
+        with open(filename, "r") as fid:
+            js_dict = json.load(fid)
+
+        for key, value in js_dict.items():
+            if key in ["datum_crs", "utm_crs"]:
+                if isinstance(value, dict):
+                    setattr(self, key, CRS.from_json_dict(value))
+                else:
+                    setattr(self, key, CRS.from_user_input(value))
+            else:
+                setattr(self, key, value)

mtpy/imaging/mtplot_tools/plot_settings.py

@@ -17,6 +17,7 @@
 
 import mtpy.imaging.mtcolors as mtcl
 
+
 # =============================================================================
 # ==============================================================================
 # Plot settings
@@ -273,8 +274,8 @@ def set_phase_limits(self, phase, mode="od"):
                     ph_min = round(ph_min / 5) * 5
                 ph_max = max(
                     [
-                        np.nanmax(phase[:, 0, 1]),
-                        np.nanmax(phase[:, 1, 0] + 180),
+                        np.nanmax(phase[nz_xy, 0, 1]),
+                        np.nanmax(phase[nz_yx, 1, 0] + 180),
                     ]
                 )
                 if ph_max < 91:
@@ -373,7 +374,6 @@ def font_dict(self):
         return {"size": self.font_size + 2, "weight": self.font_weight}
 
     def make_pt_cb(self, ax):
-
         cmap = mtcl.cmapdict[self.ellipse_cmap]
         if "seg" in self.ellipse_cmap:
             # normalize the colors

mtpy/imaging/plot_mt_response.py

@@ -149,7 +149,6 @@ def rotation_angle(self, theta_r):
         only a single value is allowed
         """
         if not theta_r == 0:
-
             self.Z.rotate(theta_r)
             self.Tipper.rotate(theta_r)
             self.pt.rotate(theta_r)
@@ -249,7 +248,6 @@ def _setup_subplots(self):
         return label_coords
 
     def _plot_resistivity(self, axr, period, z_obj, mode="od"):
-
         if mode == "od":
             comps = ["xy", "yx"]
             props = [
@@ -262,8 +260,6 @@ def _plot_resistivity(self, axr, period, z_obj, mode="od"):
                 self.xy_error_bar_properties,
                 self.yx_error_bar_properties,
             ]
-        res_limits = self.set_resistivity_limits(z_obj.resistivity, mode=mode)
-        x_limits = self.set_period_limits(period)
 
         eb_list = []
         label_list = []
@@ -282,8 +278,8 @@ def _plot_resistivity(self, axr, period, z_obj, mode="od"):
 
         axr.set_yscale("log", nonpositive="clip")
         axr.set_xscale("log", nonpositive="clip")
-        axr.set_xlim(x_limits)
-        axr.set_ylim(res_limits)
+        axr.set_xlim(self.x_limits)
+        axr.set_ylim(self.res_limits)
         axr.grid(
             True, alpha=0.25, which="both", color=(0.25, 0.25, 0.25), lw=0.25
         )
@@ -318,8 +314,8 @@ def _plot_phase(self, axp, period, z_obj, mode="od", index=0):
                 self.xy_error_bar_properties,
                 self.yx_error_bar_properties,
             ]
-        phase_limits = self.set_phase_limits(z_obj.phase, mode=mode)
 
+        phase_limits = self.set_phase_limits(z_obj.phase, mode=mode)
         for comp, prop in zip(comps, props):
             if comp == "yx":
                 plot_phase(
@@ -385,9 +381,12 @@ def _plot_determinant(self):
             borderpad=0.02,
         )
 
-        self.axr.set_ylim(
-            self.set_resistivity_limits(self.Z.resistivity, mode="det")
-        )
+        if self.res_limits is None:
+            self.axr.set_ylim(
+                self.set_resistivity_limits(self.Z.resistivity, mode="det")
+            )
+        else:
+            self.axr.set_ylim(self.res_limits)
 
     def _plot_tipper(self):
         if self.Tipper is not None:
@@ -406,7 +405,8 @@ def _plot_tipper(self):
                 self.axt.set_xlim(
                     np.log10(self.x_limits[0]), np.log10(self.x_limits[1])
                 )
-
+                if self.tipper_limits is not None:
+                    self.axt.set_ylim(self.tipper_limits)
                 self.axt.yaxis.set_major_locator(MultipleLocator(0.2))
                 self.axt.yaxis.set_minor_locator(MultipleLocator(0.1))
                 self.axt.set_xlabel("Period (s)", fontdict=self.font_dict)
@@ -420,11 +420,13 @@ def _plot_tipper(self):
     def _plot_pt(self):
         # ----plot phase tensor ellipse---------------------------------------
         if self.plot_pt:
-
             color_array = self.get_pt_color_array(self.pt)
 
             # -------------plot ellipses-----------------------------------
-            self.cbax, self.cbpt, = plot_pt_lateral(
+            (
+                self.cbax,
+                self.cbpt,
+            ) = plot_pt_lateral(
                 self.axpt,
                 self.pt,
                 color_array,

mtpy/modeling/structured_mesh_3d.py

@@ -572,8 +572,7 @@ def make_mesh(self, verbose=True):
         for s_north in sorted(self.station_locations.model_north):
             try:
                 node_index = np.where(
-                    abs(s_north - self.grid_north)
-                    < 0.02 * self.cell_size_north
+                    abs(s_north - self.grid_north) < 0.02 * self.cell_size_north
                 )[0][0]
                 if s_north - self.grid_north[node_index] > 0:
                     self.grid_north[node_index] -= 0.02 * self.cell_size_north
@@ -585,7 +584,7 @@ def make_mesh(self, verbose=True):
         if self.z_mesh_method == "custom":
             if self.grid_z is None:
                 self.z_mesh_method = "default"
-                self._logger.warn(
+                self._logger.warning(
                     "No grid_z provided, creating new z mesh using default method"
                 )
 
@@ -602,9 +601,7 @@ def make_mesh(self, verbose=True):
             )
 
         # compute grid center
-        center_east = np.round(
-            self.grid_east.min() - self.grid_east.mean(), -1
-        )
+        center_east = np.round(self.grid_east.min() - self.grid_east.mean(), -1)
         center_north = np.round(
             self.grid_north.min() - self.grid_north.mean(), -1
         )
@@ -692,7 +689,7 @@ def add_layers_to_mesh(
             layer_thickness = layer_thickness[-1]
 
             if n_add_layers != len(add_layers):
-                self._logger.warn(
+                self._logger.warning(
                     "Updating number of layers to reflect the length of the layer thickness array"
                 )
             n_add_layers = len(add_layers)
@@ -826,8 +823,7 @@ def to_modem(self, model_fn=None, **kwargs):
         if self.res_scale.lower() == "loge":
             write_res_model = np.log(self.res_model[::-1, :, :])
         elif (
-            self.res_scale.lower() == "log"
-            or self.res_scale.lower() == "log10"
+            self.res_scale.lower() == "log" or self.res_scale.lower() == "log10"
         ):
             write_res_model = np.log10(self.res_model[::-1, :, :])
         elif self.res_scale.lower() == "linear":
@@ -945,9 +941,7 @@ def from_modem(self, model_fn=None):
         self.nodes_east = np.array(
             [float(nn) for nn in ilines[3].strip().split()]
         )
-        self.nodes_z = np.array(
-            [float(nn) for nn in ilines[4].strip().split()]
-        )
+        self.nodes_z = np.array([float(nn) for nn in ilines[4].strip().split()])
 
         self.res_model = np.zeros((n_north, n_east, n_z))
 
@@ -1631,10 +1625,10 @@ def add_topography_to_model(
             )
 
         if self.n_air_layers is None or self.n_air_layers == 0:
-            self._logger.warn(
+            self._logger.warning(
                 "No air layers specified, so will not add air/topography !!!"
             )
-            self._logger.warn(
+            self._logger.warning(
                 "Only bathymetry will be added below according to the topofile: sea-water low resistivity!!!"
             )
 
@@ -1679,9 +1673,7 @@ def add_topography_to_model(
                 # adjust level to topography min
                 if max_elev is not None:
                     self.grid_z -= max_elev
-                    ztops = np.where(
-                        self.surface_dict["topography"] > max_elev
-                    )
+                    ztops = np.where(self.surface_dict["topography"] > max_elev)
                     self.surface_dict["topography"][ztops] = max_elev
                 else:
                     self.grid_z -= topo_core.max()
@@ -2770,7 +2762,7 @@ def from_ws3dinv(self, model_fn):
             lagrange = float(info[8])
         except IndexError:
             self._logger.warning("Did not get Lagrange Multiplier")
-            lagrange = 1.
+            lagrange = 1.0
 
         # get lengths of things
         n_north, n_east, n_z, n_res = np.array(
@@ -2792,7 +2784,7 @@ def from_ws3dinv(self, model_fn):
                 self._nodes_north[count_n] = float(north_node)
                 count_n += 1
             line_index += 1
-        self.grid_north = np.insert(np.cumsum(self.nodes_north),0,0)
+        self.grid_north = np.insert(np.cumsum(self.nodes_north), 0, 0)
 
         count_e = 0  # number of east nodes found
         while count_e < n_east:
@@ -2801,7 +2793,7 @@ def from_ws3dinv(self, model_fn):
                 self._nodes_east[count_e] = float(east_node)
                 count_e += 1
             line_index += 1
-        self.grid_east = np.insert(np.cumsum(self.nodes_east),0,0)
+        self.grid_east = np.insert(np.cumsum(self.nodes_east), 0, 0)
 
         count_z = 0  # number of vertical nodes
         zdep = 0

setup.py

@@ -65,7 +65,7 @@
     test_suite="tests",
     tests_require=test_requirements,
     url="https://github.com/MTgeophysics/mtpy-v2",
-    version="2.0.6",
+    version="2.0.7",
     zip_safe=False,
     package_data={"": []},
 )