From 0fb8049b18b81440390d1af53ff017ff4280347a Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 13 Mar 2023 17:30:18 -0700 Subject: [PATCH 1/5] Add flag to drop null prisms when converting to pyvista grid --- harmonica/_forward/prism_layer.py | 18 ++++++++++++++++-- harmonica/tests/test_prism_layer.py | 27 ++++++++++++++++++++++++++- 2 files changed, 42 insertions(+), 3 deletions(-) diff --git a/harmonica/_forward/prism_layer.py b/harmonica/_forward/prism_layer.py index 44a8b9e10..ff8a791d6 100644 --- a/harmonica/_forward/prism_layer.py +++ b/harmonica/_forward/prism_layer.py @@ -470,10 +470,16 @@ def get_prism(self, indices): top = self._obj.top.values[indices] return west, east, south, north, bottom, top - def to_pyvista(self): + def to_pyvista(self, drop_null_prisms=False): """ Return a pyvista UnstructuredGrid to plot the PrismLayer + Parameters + ---------- + drop_null_prisms : bool (optional) + If True, prisms with zero volume are not going to be converted. + Default False. + Returns ------- pv_grid : :class:`pyvista.UnstructuredGrid` @@ -481,10 +487,18 @@ def to_pyvista(self): layer as a hexahedron along with their properties. """ prisms = self._to_prisms() + null_prisms = np.zeros_like(prisms[:, 0], dtype=bool) + if drop_null_prisms: + bottom, top = prisms[:, -2], prisms[:, -1] + null_prisms = top == bottom + prisms = prisms[np.logical_not(null_prisms)] + # Define properties properties = None if self._obj.data_vars: properties = { - data_var: np.asarray(self._obj[data_var]).ravel() + data_var: np.asarray(self._obj[data_var]).ravel()[ + np.logical_not(null_prisms) + ] for data_var in self._obj.data_vars } return prism_to_pyvista(prisms, properties=properties) diff --git a/harmonica/tests/test_prism_layer.py b/harmonica/tests/test_prism_layer.py index 134e89fdf..bd9784810 100644 --- a/harmonica/tests/test_prism_layer.py +++ b/harmonica/tests/test_prism_layer.py @@ -419,7 +419,7 @@ def test_to_pyvista(dummy_layer, properties): assert pv_grid.n_points == 20 * 8 # Check coordinates of prisms for i, prism in enumerate(layer.prism_layer._to_prisms()): - npt.assert_allclose(prism, pv_grid.cell_bounds(i)) + npt.assert_allclose(prism, pv_grid.cell[i].bounds) # Check properties of the prisms if properties is None: assert pv_grid.n_arrays == 0 @@ -431,6 +431,31 @@ def test_to_pyvista(dummy_layer, properties): npt.assert_allclose(pv_grid.get_array("density"), layer.density.values.ravel()) +@pytest.mark.skipif(pyvista is None, reason="requires pyvista") +def test_to_pyvista_zero_volume(dummy_layer): + """ + Test the conversion of the prism layer to pyvista.UnstructuredGrid when + some prisms have zero volume + """ + (easting, northing), surface, reference, density = dummy_layer + # Build the layer with or without properties + properties = {"density": density} + layer = prism_layer((easting, northing), surface, reference, properties=properties) + # Assign zero volume to some prisms + layer.top.values[0, 0] = layer.bottom.values[0, 0] + layer.top.values[2, 1] = layer.bottom.values[2, 1] + # Convert the layer to pyvista UnstructuredGrid + pv_grid = layer.prism_layer.to_pyvista(drop_null_prisms=True) + # Check properties of the pyvista grid + expected_n_prisms = 20 - 2 + assert pv_grid.n_cells == expected_n_prisms + assert pv_grid.n_points == expected_n_prisms * 8 + assert pv_grid.n_arrays == 1 + assert pv_grid.array_names == ["density"] + assert pv_grid.get_array("density").ndim == 1 + assert pv_grid.get_array("density").size == expected_n_prisms + + @pytest.mark.skipif(ProgressBar is None, reason="requires numba_progress") @pytest.mark.use_numba def test_progress_bar(dummy_layer): From 1adfae716f874c8ee73a7e419369a20ab4642cd2 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 13 Mar 2023 17:32:19 -0700 Subject: [PATCH 2/5] Improve docstring --- harmonica/_forward/prism_layer.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/harmonica/_forward/prism_layer.py b/harmonica/_forward/prism_layer.py index ff8a791d6..18480c305 100644 --- a/harmonica/_forward/prism_layer.py +++ b/harmonica/_forward/prism_layer.py @@ -477,7 +477,9 @@ def to_pyvista(self, drop_null_prisms=False): Parameters ---------- drop_null_prisms : bool (optional) - If True, prisms with zero volume are not going to be converted. + If True, prisms with zero volume are not going to be included in + the :class:`pyvista.UnstructuredGrid`. + If False, every prism in the layer will be included. Default False. Returns From 065d7bebcbc227e98f83d663548106cd0475c761 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Mar 2023 11:14:48 -0700 Subject: [PATCH 3/5] Discard prisms with nans in their top and bottom boundaries Include prisms that have either a top or bottom boundary equal to nan into the list of null prisms. Update tests to include them as well. --- harmonica/_forward/prism_layer.py | 7 ++++--- harmonica/tests/test_prism_layer.py | 4 +++- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/harmonica/_forward/prism_layer.py b/harmonica/_forward/prism_layer.py index 37db3e428..5b4b488a8 100644 --- a/harmonica/_forward/prism_layer.py +++ b/harmonica/_forward/prism_layer.py @@ -477,8 +477,9 @@ def to_pyvista(self, drop_null_prisms=False): Parameters ---------- drop_null_prisms : bool (optional) - If True, prisms with zero volume are not going to be included in - the :class:`pyvista.UnstructuredGrid`. + If True, prisms with zero volume or with any :class:`numpy.nan` as + their top or bottom boundaries won't be included in the + :class:`pyvista.UnstructuredGrid`. If False, every prism in the layer will be included. Default False. @@ -492,7 +493,7 @@ def to_pyvista(self, drop_null_prisms=False): null_prisms = np.zeros_like(prisms[:, 0], dtype=bool) if drop_null_prisms: bottom, top = prisms[:, -2], prisms[:, -1] - null_prisms = top == bottom + null_prisms = (top == bottom) | (np.isnan(top)) | (np.isnan(bottom)) prisms = prisms[np.logical_not(null_prisms)] # Define properties properties = None diff --git a/harmonica/tests/test_prism_layer.py b/harmonica/tests/test_prism_layer.py index bd9784810..2660d749f 100644 --- a/harmonica/tests/test_prism_layer.py +++ b/harmonica/tests/test_prism_layer.py @@ -444,10 +444,12 @@ def test_to_pyvista_zero_volume(dummy_layer): # Assign zero volume to some prisms layer.top.values[0, 0] = layer.bottom.values[0, 0] layer.top.values[2, 1] = layer.bottom.values[2, 1] + layer.top.values[3, 2] = np.nan + layer.bottom.values[3, 3] = np.nan # Convert the layer to pyvista UnstructuredGrid pv_grid = layer.prism_layer.to_pyvista(drop_null_prisms=True) # Check properties of the pyvista grid - expected_n_prisms = 20 - 2 + expected_n_prisms = 20 - 4 assert pv_grid.n_cells == expected_n_prisms assert pv_grid.n_points == expected_n_prisms * 8 assert pv_grid.n_arrays == 1 From 5c3c5e230225ffe4159949fbd1ae0ef9a988cad3 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Mar 2023 14:33:25 -0700 Subject: [PATCH 4/5] Make drop_null_prisms True by default Update tests functions --- harmonica/_forward/prism_layer.py | 2 +- harmonica/tests/test_prism_layer.py | 16 +++++++++++----- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/harmonica/_forward/prism_layer.py b/harmonica/_forward/prism_layer.py index 5b4b488a8..770d54344 100644 --- a/harmonica/_forward/prism_layer.py +++ b/harmonica/_forward/prism_layer.py @@ -470,7 +470,7 @@ def get_prism(self, indices): top = self._obj.top.values[indices] return west, east, south, north, bottom, top - def to_pyvista(self, drop_null_prisms=False): + def to_pyvista(self, drop_null_prisms=True): """ Return a pyvista UnstructuredGrid to plot the PrismLayer diff --git a/harmonica/tests/test_prism_layer.py b/harmonica/tests/test_prism_layer.py index 2660d749f..d94ae8377 100644 --- a/harmonica/tests/test_prism_layer.py +++ b/harmonica/tests/test_prism_layer.py @@ -401,11 +401,13 @@ def test_prism_layer_gravity_density_nans(field, dummy_layer, prism_layer_with_h @pytest.mark.skipif(pyvista is None, reason="requires pyvista") @pytest.mark.parametrize("properties", (False, True)) -def test_to_pyvista(dummy_layer, properties): +@pytest.mark.parametrize("drop_null_prisms", (False, True)) +def test_to_pyvista(dummy_layer, properties, drop_null_prisms): """ Test the conversion of the prism layer to pyvista.UnstructuredGrid """ (easting, northing), surface, reference, density = dummy_layer + reference -= 1 # lower the reference so we don't get prisms with zero volume # Build the layer with or without properties if properties: properties = {"density": density} @@ -413,7 +415,7 @@ def test_to_pyvista(dummy_layer, properties): properties = None layer = prism_layer((easting, northing), surface, reference, properties=properties) # Convert the layer to pyvista UnstructuredGrid - pv_grid = layer.prism_layer.to_pyvista() + pv_grid = layer.prism_layer.to_pyvista(drop_null_prisms=drop_null_prisms) # Check properties of the pyvista grid assert pv_grid.n_cells == 20 assert pv_grid.n_points == 20 * 8 @@ -432,12 +434,14 @@ def test_to_pyvista(dummy_layer, properties): @pytest.mark.skipif(pyvista is None, reason="requires pyvista") -def test_to_pyvista_zero_volume(dummy_layer): +@pytest.mark.parametrize("drop_null_prisms", (False, True)) +def test_to_pyvista_drop_null_prisms(dummy_layer, drop_null_prisms): """ Test the conversion of the prism layer to pyvista.UnstructuredGrid when some prisms have zero volume """ (easting, northing), surface, reference, density = dummy_layer + reference -= 1 # lower the reference so we don't get prisms with zero volume # Build the layer with or without properties properties = {"density": density} layer = prism_layer((easting, northing), surface, reference, properties=properties) @@ -447,9 +451,11 @@ def test_to_pyvista_zero_volume(dummy_layer): layer.top.values[3, 2] = np.nan layer.bottom.values[3, 3] = np.nan # Convert the layer to pyvista UnstructuredGrid - pv_grid = layer.prism_layer.to_pyvista(drop_null_prisms=True) + pv_grid = layer.prism_layer.to_pyvista(drop_null_prisms=drop_null_prisms) # Check properties of the pyvista grid - expected_n_prisms = 20 - 4 + expected_n_prisms = 20 + if drop_null_prisms: + expected_n_prisms -= 4 assert pv_grid.n_cells == expected_n_prisms assert pv_grid.n_points == expected_n_prisms * 8 assert pv_grid.n_arrays == 1 From 33e55840b70bc86bb903ef4b7d96bf4bc4e42718 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Mar 2023 14:34:25 -0700 Subject: [PATCH 5/5] Update docstrings --- harmonica/_forward/prism_layer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/_forward/prism_layer.py b/harmonica/_forward/prism_layer.py index 770d54344..8cecf9cae 100644 --- a/harmonica/_forward/prism_layer.py +++ b/harmonica/_forward/prism_layer.py @@ -481,7 +481,7 @@ def to_pyvista(self, drop_null_prisms=True): their top or bottom boundaries won't be included in the :class:`pyvista.UnstructuredGrid`. If False, every prism in the layer will be included. - Default False. + Default True. Returns -------