fixed auto extent and made some other minor changes

gempy/data_management.py

@@ -134,15 +134,6 @@ def __init__(self,
         self.interfaces['isFault'] = self.interfaces['isFault'].astype('bool')
         self.orientations['isFault'] = self.orientations['isFault'].astype('bool')
 
-    def get_auto_extent(self):
-        X = np.concatenate((self.orientations['X'], self.interfaces['X']))
-        Y = np.concatenate((self.orientations['Y'], self.interfaces['Y']))
-        Z = np.concatenate((self.orientations['Z'], self.interfaces['Z']))
-        x_min, x_max, y_min, y_max, z_min, z_max = X.min(), X.max(), Y.min(), Y.max(), Z.min(), Z.max()
-        # some extra space at borders of model
-        xspace, yspace, zspace = (x_max-x_min) / 10, (y_max-y_min) / 10, (z_max-z_min) / 10
-        return np.array([x_min-xspace, x_max+xspace, y_min-yspace, y_max+yspace, z_min-zspace, z_max+zspace]).astype(int)
-
     def set_basement(self):
 
         try:
@@ -429,7 +420,10 @@ def import_data_csv(self, path_i, path_o, **kwargs):
             interfaces_read = interfaces_read.assign(**dict.fromkeys(c[~np.in1d(c, interfaces_read.columns)], False))
             self.set_interfaces(interfaces_read, append=True)
             #self.interfaces[interfaces_read.columns] = interfaces_read[interfaces_read.columns]
-            #gagag
+
+        # set auto extent if it is not defined yet
+        if self.extent is None:
+            self.extent = self.get_auto_extent()
         self.update_df()
 
     def modify_interface(self, index, **kwargs):
@@ -1122,6 +1116,37 @@ def _read_vox(self, path):
         return block_geomodeller
 
 
+    def get_auto_extent(self):
+        """
+        calculates the extent based in the minumum and maximum coordinates in the input data.
+        Returns: np.array([x_min, x_max, y_min, y_max, z_min, z_max])
+        """
+
+        def ru(x, n=100):
+            """rounds up to the nearest n of x+40"""
+            return x+100 if x % n < 40 else x+n - x % n
+
+        def rd(x, n=100, stop0=None):
+            """rounds down to the nearest n of x+40, if stop0 zero is the minimum"""
+            if x % n < 40:
+                return x-100 if x > 40 else 0 if stop0 else ru(x)-n
+            elif x < 0:
+                return 0 if stop0 else ru(x)-n
+            else:
+                return x-x % n
+
+        X = np.concatenate((self.orientations['X'], self.interfaces['X']))
+        Y = np.concatenate((self.orientations['Y'], self.interfaces['Y']))
+        Z = np.concatenate((self.orientations['Z'], self.interfaces['Z']))
+
+        ext = np.array([X.min(), X.max(), Y.min(), Y.max(), Z.min(), Z.max()])
+
+        if not any(x < 0 for x in ext): #avoids rounding below 0 if there is no negative number in the input data
+            return np.array([rd(ext[0], stop0=1), ru(ext[1]), rd(ext[2], stop0=1), ru(ext[3]), rd(ext[4], stop0=1), ru(ext[5])])
+        else:
+            return np.array([rd(ext[0]), ru(ext[1]), rd(ext[2]), ru(ext[3]), rd(ext[4]), ru(ext[5])])
+
+
 def get_orientation(normal):
     """Get orientation (dip, azimuth, polarity ) for points in all point set"""
     #    if "normal" not in dir(self):

gempy/gempy_front.py

@@ -182,7 +182,7 @@ def create_data(extent=None, resolution=(50, 50, 50), **kwargs):
 def create_from_geomodeller_xml(fp, resolution=(50, 50, 50), return_xml=False, **kwargs):
     """
     Creates InputData object from a GeoModeller xml file. Automatically extracts and sets model extent, interface
-    and orientation data as well as the stratigraphic pile.
+    and orientation data (as well as the stratigraphic pile, soon).
 
     Args:
         fp (str): Filepath for the GeoModeller xml file to be read.
@@ -194,6 +194,7 @@ def create_from_geomodeller_xml(fp, resolution=(50, 50, 50), return_xml=False, *
     Returns:
         gp.data_management.InputData
     """
+    #Todo finish building straigraphic pile
     gmx = _ReadGeoModellerXML(fp)  # instantiate parser class with filepath of xml
 
     # instantiate InputData object with extent and resolution

gempy/plotting/plot.py

@@ -621,6 +621,7 @@ def plot_data(geo_data, direction="y", data_type = 'all', series="all", legend_f
     plot = PlotData2D(geo_data)
 
     # TODO saving options
+
     return plot.plot_data(direction=direction, data_type=data_type, series=series, legend_font_size=legend_font_size, **kwargs)
 
 
@@ -712,18 +713,20 @@ def plot_topology(geo_data, G, centroids, direction="y"):
     PlotData2D.plot_topo_g(geo_data, G, centroids, direction=direction)
 
 
-def plot_stereonet(geo_data, series_only=False, litho=None, planes=True, poles=True, single_plots=False, show_density=False):
+def plot_stereonet(geo_data, litho=None, series_only=False, planes=True, poles=True, single_plots=False, show_density=False, legend=True):
     '''
-    Plot an equal-area projection of the orientations dataframe using mplstereonet.
+    Plots an equal-area projection of the orientations dataframe using mplstereonet.
+    Only works after assigning the series for the right color assignment.
 
     Args:
         geo_data (gempy.DataManagement.InputData): Input data of the model
-        series_only: To select whether a stereonet is plotted per series or per formation
-        litho: selection of formation or series names, as list. If None, all are plotted
-        planes: If True, azimuth and dip are plotted as great circles
-        poles: If True, pole points (plane normal vectors) of azimuth and dip are plotted
-        single_plots: If True, each formation is plotted in a single stereonet
-        show_density: If True, density contour plot around the pole points is shown
+        litho (list): selection of formation or series names. If None, all are plotted
+        series_only (bool): to decide if the data is plotted per series or per formation
+        planes (bool): plots azimuth and dip as great circles
+        poles (bool): plots pole points (plane normal vectors) of azimuth and dip
+        single_plots (bool): plots each formation in a single stereonet
+        show_density (bool): shows density contour plot around the pole points
+        legend (bool): shows legend
 
     Returns:
         None
@@ -769,10 +772,10 @@ def plot_stereonet(geo_data, series_only=False, litho=None, planes=True, poles=T
         if poles:
             ax.pole(df_sub['azimuth'] - 90, df_sub['dip'], marker='o', markersize=7,
                     markerfacecolor=cmap(df_sub['formation_number'].values[0]),
-                    markeredgewidth=1.1, markeredgecolor='gray', label=formation+': '+'pole point')
+                    markeredgewidth=1.1, markeredgecolor='gray', label=formation)#+': '+'pole point')
         if planes:
             ax.plane(df_sub['azimuth'] - 90, df_sub['dip'], color=cmap(df_sub['formation_number'].values[0]),
-                     linewidth=1.5, label=formation+': '+'azimuth/dip')
+                     linewidth=1.5, label=formation)
         if show_density:
             if single_plots:
                 ax.density_contourf(df_sub['azimuth'] - 90, df_sub['dip'],
@@ -782,8 +785,9 @@ def plot_stereonet(geo_data, series_only=False, litho=None, planes=True, poles=T
                                     alpha=.5)
 
         fig.subplots_adjust(top=0.8)
-        handles, labels = ax.get_legend_handles_labels()
-        by_label = OrderedDict(zip(labels, handles))
-        ax.legend(by_label.values(), by_label.keys(), bbox_to_anchor=(1.9, 1.1))
+        if legend:
+            handles, labels = ax.get_legend_handles_labels()
+            by_label = OrderedDict(zip(labels, handles))
+            ax.legend(by_label.values(), by_label.keys(), bbox_to_anchor=(1.9, 1.1))
         ax.grid(True, color='black', alpha=0.25)
-
+    return fig

gempy/posterior_analysis_simple.py

@@ -14,16 +14,19 @@
 import matplotlib.pyplot as plt
 
 class Posterior:
-    def __init__(self, dbname, verbose=False, entropy=False, interp_data = None):
+
+
+    def __init__(self, dbname, verbose=False, entropy=False, interp_data=None):
         """
         Posterior database analysis for GemPy-pymc2 hdf5 databases.
         Args:
             dbname (str): Path of the hdf5 database.
-            pymc_model_f (str, optional): name of the model output function used (default: "gempy_model).
-            pymc_topo_f (str, optional): name of the topology output function used (default: "gempy_topo).
-            topology (bool, optional):  if a topology trace should be loaded from the database (default: False).
+            entropy (bool): if true, all postmodels are calculated (may take some time!) to visualize entropy
             verbose (bool, optional): Verbosity switch.
         """
+        if entropy:
+            warnings.warn('All post models are calculated. Based on the model complexity and the number of iterations, '
+                          'this may take some time!')
         # TODO: Add a method to set the lith_block and fault_block
         self.interp_data = interp_data
         self.verbose = verbose
@@ -37,10 +40,6 @@ def __init__(self, dbname, verbose=False, entropy=False, interp_data = None):
         # load input data
         self.input_data = self.db.input_data.gettrace()
 
-        self.lith_prob = None
-        self.ie = None
-        self.ie_total = None
-
         if entropy is True:
             self.lbs, self.fbs = self.all_post_models()
 
@@ -52,19 +51,21 @@ def __init__(self, dbname, verbose=False, entropy=False, interp_data = None):
 
             self.ie_total = self.calculate_ie_total()
 
-    def plot_lith_entropy(self, resolution):
+    def plot_lith_entropy(self, extent, resolution):
         '''plots information entropy in middle of block model in y-direction'''
         y = int(resolution[1] / 2)
         ie_reshaped = self.lb_ie.reshape(resolution)
         plt.figure(figsize=(15, 5))
         fig, ax = plt.subplots(1)
         plt.rc('xtick', labelsize=18)
         plt.rc('ytick', labelsize=18)
-        plt.imshow(ie_reshaped[:, y, :].T, origin="lower", cmap="viridis")
+        plt.imshow(ie_reshaped[:, y, :].T, origin="lower", cmap="viridis", extent=[extent[0], extent[1], extent[4], extent[5]])
+        plt.xlabel('X')
+        plt.ylabel('Y')
         plt.colorbar()
         return fig
 
-    def plot_fault_entropy(self, resolution):
+    def plot_fault_entropy(self, extent, resolution):
         '''plots information entropy in middle of block model in y-direction'''
         y = int(resolution[1] / 2)
         # print(y, resolution)
@@ -74,7 +75,9 @@ def plot_fault_entropy(self, resolution):
         fig, ax = plt.subplots(1)
         plt.rc('xtick', labelsize=18)
         plt.rc('ytick', labelsize=18)
-        plt.imshow(ie_reshaped[:, y, :].T, origin="lower", cmap="viridis")
+        plt.imshow(ie_reshaped[:, y, :].T, origin="lower", cmap="viridis", extent=[extent[0], extent[1], extent[4], extent[5]])
+        plt.xlabel('X')
+        plt.ylabel('Y')
         plt.colorbar()
         return fig
 
@@ -115,7 +118,6 @@ def all_post_models(self):
         return lbs, fbs
 
     def compute_prob(self, lith_blocks):
-        """Blocks must be just the lith blocks!"""
         lith_id = np.unique(lith_blocks)
         # print(len(lith_id))
         # lith_count = np.zeros_like(lith_blocks[0:len(lith_id)])

gempy/utils/fishdist.py

@@ -6,13 +6,13 @@
 from mpl_toolkits.mplot3d import proj3d
 import mplstereonet
 
-'''This solution for vMF sampling originally appeared here: 
+'''Parts of this solution for vMF sampling originally appeared here: 
 - https://github.com/pymc-devs/pymc3/issues/2458 and
 - https://github.com/jasonlaska/spherecluster/blob/master/spherecluster/util.py'''
 
 
 class vMF():
-    '''draws and visualizes samples from a van mises fisher distribution based on mean, concentration and number of
+    '''draws and visualizes samples from a von mises fisher distribution based on mean, concentration and number of
     samples for stochastic simulations with orientation uncertainty'''
 
     def __init__(self, mu, kappa, num_samples):
@@ -126,7 +126,7 @@ def plot_vMF_2D(self, poles=True):
             ax.pole(mean[0] - 90, mean[1], color='r', markersize=6, label='mean')
         ax.grid()
         ax.density_contourf(points_sph[:, 0] - 90, points_sph[:, 1], measurement='poles', cmap='inferno', alpha=0.7)
-        ax.set_title('kappa = '+str(self.kappa),y=1.1)
+        ax.set_title('kappa = '+str(self.kappa), y=1.2)
         return fig
 
     def xyz_to_spherical_coordinates(self, gamma1):