Merge branch 'dev' of https://github.com/gimli-org/gimli into dev

.github/workflows/main.yml

@@ -39,7 +39,7 @@ jobs:
       OPENBLAS_CORETYPE: "ARMV8" # current bug in OpenBlas core detection
     steps:
       - name: Run pg.test()
-        run: python3 -c "import pygimli; pygimli.test(show=False)"
+        run: python3 -c "import pygimli; pygimli.test(show=False, abort=True)"
         working-directory: source
       - name: Install as development package
         working-directory: source

doc/examples/3_ert/plot_02_ert_field_data.py

@@ -9,8 +9,11 @@
 - 2D inversion with topography
 - geometric factor generation
 - topography effect
+
+The data is the profile 11 already shown by Günther et al. (2006, Fig. 11).
 """
-# sphinx_gallery_thumbnail_number = 5
+# sphinx_gallery_thumbnail_number = 7
+import matplotlib.pyplot as plt
 import pygimli as pg
 from pygimli.physics import ert
 
@@ -21,36 +24,46 @@
 data = pg.getExampleData('ert/slagdump.ohm', verbose=True)
 print(data)
 
+###############################################################################
+# Let us first have a look at the topography contained in the data
+plt.plot(pg.x(data), pg.z(data), 'x-')
+
 ###############################################################################
 # The data file does not contain geometric factors (token field 'k'),
 # so we create them based on the given topography.
+k0 = ert.createGeometricFactors(data)  # the analytical one
 data['k'] = ert.createGeometricFactors(data, numerical=True)
 
-###############################################################################
-# We initialize the ERTManager for further steps and eventually inversion.
-mgr = ert.ERTManager(sr=False)
-
 ###############################################################################
 # It might be interesting to see the topography effect, i.e the ratio between
-# the numerically computed geometry factor and the analytical formula
-k0 = ert.createGeometricFactors(data)
-_ = ert.showData(data, vals=k0/data['k'], label='Topography effect')
+# the numerically computed geometry factor and the analytical formula after
+# Rücker et al. (2006). We display it using a colormap with neutral white.
+_ = ert.showData(data, vals=k0/ data['k'], label='Topography effect',
+                 cMin=2/3, cMax=3/2, logScale=True, cMap="bwr")
 
 ###############################################################################
-# The data container has no apparent resistivities (token field 'rhoa') yet.
-# We can let the Manager fix this later for us (as we now have the 'k' field),
-# or we do it manually.
-mgr.checkData(data)
-print(data)
+# We can now compute the apparent resistivity and display it, once with the
+# wrong analytical formula and once with the numerical values in data['k']
+data['rhoa'] = data['r'] * data['k']
+kw = dict(cMin=6, cMax=33)
+fig, ax = plt.subplots(ncols=2)
+data.show(data['r']*k0, ax=ax[0], **kw);
+data.show(ax=ax[1], **kw)
+ax[0].set_title('Uncorrected')
+ax[1].set_title('Corrected');
 
 ###############################################################################
 # The data container does not necessarily contain data errors data errors
 # (token field 'err'), requiring us to enter data errors. We can let the
 # manager guess some defaults for us automaticly or set them manually
-data['err'] = ert.estimateError(data, relativeError=0.03, absoluteUError=5e-5)
-# or manually:
-# data['err'] = data_errors  # somehow
-_ = ert.show(data, data['err']*100)
+data.estimateError(relativeError=0.03, absoluteUError=5e-5)
+# which internally calls
+# data['err'] = ert.estimateError(data, ...)  # can also set manually
+_ = data.show(data['err']*100, label='error estimate (%)')
+
+###############################################################################
+# We initialize the ERTManager for further steps and eventually inversion.
+mgr = ert.ERTManager(data)
 
 ###############################################################################
 # Now the data have all necessary fields ('rhoa', 'err' and 'k') so we can run
@@ -59,8 +72,7 @@
 #
 mod = mgr.invert(data, lam=10, verbose=True,
                  paraDX=0.3, paraMaxCellSize=10, paraDepth=20, quality=33.6)
-
-_ = mgr.showResult()
+ax, cb = mgr.showResult()
 
 ###############################################################################
 # We can view the resulting model in the usual way.
@@ -69,4 +81,4 @@
 
 ###############################################################################
 # Or just plot the model only using your own options.
-_ = mgr.showResult(mod, cMin=5, cMax=30, cMap="Spectral_r", logScale=True)
+ax, cb = mgr.showResult(mod, cMin=5, cMax=30, cMap="Spectral_r", logScale=True)

pygimli/physics/ert/ertManager.py

@@ -390,9 +390,13 @@ def showModel(self, model=None, ax=None, elecs=True, **kwargs):
             _, ax = pg.plt.subplots()
 
         kwargs.setdefault("coverage", self.coverage())
-        ax, cBar = self.fop.drawModel(ax, model, **kwargs)
+        color = kwargs.pop("color", "magenta")
+        ax, cBar = self.fop.drawModel(ax, model, **kwargs)        
         if elecs:
-            pg.viewer.mpl.drawSensors(ax, self.fop.data.sensors())
+            if isinstance(elecs, str):
+                color = elecs
+
+            pg.viewer.mpl.drawSensors(ax, self.fop.data.sensors(), color=color)
 
         return ax, cBar
 

pygimli/physics/ert/importData.py

@@ -334,7 +334,7 @@ def importAsciiColumns(filename, verbose=False, return_header=False):
 
         if content[0].startswith('Filename'):  # ABEM lead-in
             for n, line in enumerate(content):
-                if line.find("MeasID") > 0:
+                if line.find("MeasID") >= 0:
                     break
 
             for i in range(n):
@@ -506,3 +506,7 @@ def readAsDictionary(content, token=None, sep=None):  # obsolote due to numpy?
                 data[token[j]][i] = v
 
     return data
+
+
+if __name__ == "__main__":
+    pass

pygimli/viewer/showmesh.py

@@ -545,9 +545,8 @@ def _drawField(ax, mesh, data, kwargs):  # like view.mpl.drawField?
     if axisLabels == True and mesh.dim() == 2:
 
         try:
-            pg.viewer.mpl.adjustWorldAxes(ax,
-                                       useDepth=min(mesh.boundaryMarkers()) < 0,
-                                       xl=xl, yl=yl)
+            useDepth = min(mesh.boundaryMarkers()) < 0 and max(pg.y(mesh)) <= 0
+            pg.viewer.mpl.adjustWorldAxes(ax, useDepth=useDepth, xl=xl, yl=yl)
         except BaseException:
             pass
     else: