various tweaks to pass tests and improve usability

CHANGES.md

@@ -1,12 +1,14 @@
 # Changelog
 
-## 0.5.2 — coming February 2022
+## 0.5.2 — 18 February 2022
 
-- **Breaking change:** Mode is now 'same' by default on reflectivity functions. If you were assuming `mode='valid'` you should change your code.
+- **Breaking change:** Mode is now `'same'` by default on reflectivity functions. If you were assuming `mode='valid'` you should change your code.
 - `reflection.reflectivity()` now works properly on 2D panels of Vp, Vs, and rho.
+- `reflection.acoustic_reflectivity()` now optionally accepts *only* Vp (or impedance) or Rho. You should pass both if you have them, but if you only pass one, reflectivity will be computed from that alone.
 - `reflection.convolve()` now works properly on 2D and 3D reflectivity series, and even works with 2D wavelet banks (provided the bank is shorter in time than the model, which it usually will be).
 - `reflection.elastic_impedance()` now works properly on 2D panels of Vp, Vs, and rho.
 - We have started adding better documentation; check it out at [code.agilescientific.com/bruges](https://code.agilescientific.com/bruges).
+- `models.wedge()` now behaves as expected if you pass a NumPy array of `int` values to `strat`.
 
 ## 0.5.1 — 17 February 2022
 

bruges/attribute/energy.py

@@ -25,9 +25,8 @@ def energy(traces, duration, dt=1):
     Returns:
         ndarray: An array the same dimensions as the input array.
     """
-    data = traces.astype(np.float).reshape(-1, traces.shape[-1])
+    # data = traces.astype(np.float).reshape(-1, traces.shape[-1])
     n_samples = int(duration / dt)
     window = np.ones(n_samples) / n_samples
-    energy = convolve(data**2, window)
-    return energy.reshape(traces.shape)
+    return convolve(traces**2, window)
 

bruges/filters/convolve.py

@@ -41,8 +41,3 @@ def convolve(reflectivity, wavelet, axis=-1):
     pos = wavelet.ndim - 1
 
     return np.moveaxis(syn.reshape(outshape), pos, axis)
-
-
-def apply_along_axis(func_1d, arr, *args, **kwargs):
-    mapobj = map(lambda tr: func_1d(tr, *args, **kwargs), arr)
-    return np.array(list(mapobj))

bruges/filters/test/wavelets_test.py

@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 import unittest
 import numpy as np
 from numpy import array
@@ -64,4 +63,4 @@ def test_ormsby(self):
 
 if __name__ == '__main__':
     suite = unittest.TestLoader().loadTestsFromTestCase(WaveletTest)
-    unittest.TextTestRunner(verbosity=2).run(suite)
+    unittest.TextTestRunner(verbosity=2).run(suite)

bruges/models/test/models_test.py

@@ -44,6 +44,11 @@ def test_wedge(self):
         self.assertTrue(np.allclose(base, b))
         self.assertTrue(ref == 6)
 
+    def test_wedge_numpy(self):
+        strat = np.array([2000, 2100, 2300])
+        w, *_ = wedge(depth=10, width=7, strat=strat)
+        self.assertTrue(w.shape == (10, 7))
+
     def test_netgross(self):
         w, top, *_ = wedge(depth=10, width=7, breadth=3, strat=(10, (20, 30), 40))
         self.assertTrue(w.sum() == 6003)

bruges/models/wedge.py

@@ -35,9 +35,9 @@ def get_strat(strat,
     orders = {'nearest': 0, 'linear': 1, 'quadratic': 2, 'cubic': 3}
     order = orders.get(kind, 0)
 
-    if isinstance(strat, int) and (order == 0):
+    if np.issubdtype(type(strat), int) and (order == 0):
         out = np.repeat([strat], thickness)
-    elif isinstance(strat, float) and (order == 0):
+    elif np.issubdtype(type(strat), float) and (order == 0):
         out = np.repeat([strat], thickness)
     elif isinstance(strat, tuple) and (order == 0):
         out = np.repeat(strat, int(round(thickness/len(strat))))

bruges/reflection/reflection.py

@@ -216,17 +216,18 @@ def vectorize(func):
     @wraps(func)
     def wrapper(vp1, vs1, rho1, vp2, vs2, rho2, theta1=0, **kwargs):
 
-        new_shape = [-1] + vp1.ndim * [1]
-        theta1 = theta1.reshape(*new_shape)
-        if (np.nan_to_num(theta1) > np.pi/2.).any():
-            raise ValueError("Incidence angle theta1 must be less than 90 deg.")
-
         vp1 = np.asanyarray(vp1, dtype=float)
         vs1 = np.asanyarray(vs1, dtype=float) + 1e-12  # Prevent singular matrix.
         rho1 = np.asanyarray(rho1, dtype=float)
         vp2 = np.asanyarray(vp2, dtype=float)
         vs2 = np.asanyarray(vs2, dtype=float) + 1e-12  # Prevent singular matrix.
         rho2 = np.asanyarray(rho2, dtype=float)
+
+        new_shape = [-1] + vp1.ndim * [1]
+        theta1 = theta1.reshape(*new_shape)
+        if (np.nan_to_num(theta1) > np.pi/2.).any():
+            raise ValueError("Incidence angle theta1 must be less than 90 deg.")
+
         return func(vp1, vs1, rho1, vp2, vs2, rho2, theta1, **kwargs)
     return wrapper
 

bruges/rockphysics/elastic.py

@@ -40,17 +40,17 @@ def elastic_impedance(vp, vs, rho, theta1,
     Returns:
         ndarray: The elastic impedance log at the specficied angle or angles.
     """
-    new_shape = [-1] + vp.ndim * [1]
-    theta1 = np.radians(theta1).reshape(*new_shape)
-    if (np.nan_to_num(theta1) > np.pi/2.).any():
-        raise ValueError("Incidence angle theta1 must be less than 90 deg.")
-
     alpha = np.asanyarray(vp, dtype=float)
     beta = np.asanyarray(vs, dtype=float)
     rho = np.asanyarray(rho, dtype=float)
     op = np.sin if use_sin else np.tan
     k = np.mean(beta**2.0 / alpha**2.0) if k is None else k
 
+    new_shape = [-1] + alpha.ndim * [1]
+    theta1 = np.radians(theta1).reshape(*new_shape)
+    if (np.nan_to_num(theta1) > np.pi/2.).any():
+        raise ValueError("Incidence angle theta1 must be less than 90 deg.")
+
     a = 1 + op(theta1)**2.0
     b = -8 * k * np.sin(theta1)**2.0
     c = 1 - 4 * k * np.sin(theta1)**2.0

bruges/rockphysics/test/elastic_impedance_test.py

@@ -55,8 +55,8 @@ def test_ei_log(self):
         x3 = elastic_impedance(arr_vp, arr_vs, arr_rho, theta1=theta2)
 
         self.assertEqual(x1.shape, (ones_vp.size,))
-        self.assertEqual(x2.shape, (ones_vp.size, theta2.size))
-        self.assertEqual(x3.shape, (arr_vp.size, theta2.size))
+        self.assertEqual(x2.shape, (theta2.size, ones_vp.size))
+        self.assertEqual(x3.shape, (theta2.size, arr_vp.size))
 
     def test_ei_nan(self):
         """

bruges/util/util.py

@@ -263,7 +263,7 @@ def error_flag(pred, actual, dev=1.0, method=1):
     return flag
 
 
-def apply_along_axis(func_1d, arr, kernel, **kwargs):
+def apply_along_axis(func_1d, arr, *args, **kwargs):
     """
     Apply 1D function across 2D slice as efficiently as possible.
 
@@ -277,7 +277,7 @@ def apply_along_axis(func_1d, arr, kernel, **kwargs):
     Example
     >>> apply_along_axes(np.convolve, reflectivity_2d, wavelet, mode='same') 
     """
-    mapobj = map(lambda tr: func_1d(tr, kernel, **kwargs), arr)
+    mapobj = map(lambda tr: func_1d(tr, *args, **kwargs), arr)
     return np.array(list(mapobj))