diff --git a/CHANGES.md b/CHANGES.md index f175e3b..7cf8a44 100644 --- a/CHANGES.md +++ b/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 diff --git a/bruges/attribute/energy.py b/bruges/attribute/energy.py index 7c47389..360b578 100644 --- a/bruges/attribute/energy.py +++ b/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) \ No newline at end of file diff --git a/bruges/filters/convolve.py b/bruges/filters/convolve.py index 99872fd..835017b 100644 --- a/bruges/filters/convolve.py +++ b/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)) diff --git a/bruges/filters/test/wavelets_test.py b/bruges/filters/test/wavelets_test.py index f7039ef..d2fdc57 100644 --- a/bruges/filters/test/wavelets_test.py +++ b/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) \ No newline at end of file + unittest.TextTestRunner(verbosity=2).run(suite) diff --git a/bruges/models/test/models_test.py b/bruges/models/test/models_test.py index 5a4c5f3..1b1fbb4 100644 --- a/bruges/models/test/models_test.py +++ b/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) diff --git a/bruges/models/wedge.py b/bruges/models/wedge.py index df7fb3d..867ad2e 100644 --- a/bruges/models/wedge.py +++ b/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)))) diff --git a/bruges/reflection/reflection.py b/bruges/reflection/reflection.py index eda44bc..3ac4f0c 100644 --- a/bruges/reflection/reflection.py +++ b/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 diff --git a/bruges/rockphysics/elastic.py b/bruges/rockphysics/elastic.py index e2fe2eb..f1d135a 100644 --- a/bruges/rockphysics/elastic.py +++ b/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 diff --git a/bruges/rockphysics/test/elastic_impedance_test.py b/bruges/rockphysics/test/elastic_impedance_test.py index e187daf..2666e4e 100644 --- a/bruges/rockphysics/test/elastic_impedance_test.py +++ b/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): """ diff --git a/bruges/util/util.py b/bruges/util/util.py index 7b102c2..aff9379 100644 --- a/bruges/util/util.py +++ b/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))