Add function to create a layer of prisms #186

santisoler · 2020-09-02T21:17:13Z

Add function that creates a layer of prisms as a xr.Dataset which contains
the coordinates of the center of each prism, its top and bottom boundaries (as
xarray coordinates) and each property (as density, magnetization, etc) as data
array. Add PrismsLayer class as a dataset accesor
for xr.Datasets, which allows to define several attributes, properties and
methods for the prism layer, which can be accessed through the prisms_layer
attribute of the xr.Dataset. Add a gravity method that computes
gravitational fields generated by every prism in the layer whose top boundary,
bottom boundary or density are not np.nans.
Add two gallery examples for the prisms layer, and a new South Africa
topography dataset, which is a downsample of ETOPO1.

This PR may serve as a prototype for solving #83, but for prisms. The one for tesseroids could be equivalent to this one, but with some minor checks due to longitude continuity.

Reminders:

Run make format and make check to make sure the code follows the style guide.
Add tests for new features or tests that would have caught the bug that you're fixing.
Add new public functions/methods/classes to doc/api/index.rst and the base __init__.py file for the package.
Write detailed docstrings for all functions/classes/methods. It often helps to design better code if you write the docstrings first.
If adding new functionality, add an example to the docstring, gallery, and/or tutorials.
Add your full name, affiliation, and ORCID (optional) to the AUTHORS.md file (if you haven't already) in case you'd like to be listed as an author on the Zenodo archive of the next release.

santisoler · 2020-09-23T13:22:29Z

I suspect that ETOPO1 is referenced on the geoid (need to check), so we might want to add the geoid height (from the reference ellipsoid) to the topography + bathymetry grid before creating the terrain layer.

santisoler · 2020-09-30T20:45:43Z

@lperozzi Would you like to review this PR?

lperozzi

I will add another points that I've remarked today. If the aim is to extract the topographic correction for the observations stations from the prisms gravity forward, we should be careful with the bord effects.

lperozzi · 2020-10-13T09:38:35Z

examples/forward/prisms_topo_gravity.py

+plt.colorbar(tmp, label="mGal")
+ax.set_aspect("equal")
+plt.tight_layout()
+plt.show()


I do not have any specifical comments. The code is clear and easy to read, however I have a question.
Starting form a land survey, it is also possible to compute the prisms['top'] by gridding the elevation data, using verde for example? Instead to have a DEM or topographic grid for the region?

Sure! You'll probably want to apply some gridding or downsampling to your topographic grid in order to build the prisms: there's no point to use a very high resolution DEM to apply terrain correction, while if you have some topographic height measurements and low resolution DEM you can combine them using a gridder. But those tasks should be carried out by the user.

santisoler · 2020-10-13T20:19:00Z

Thanks for the feedback @lperozzi!

I will add another points that I've remarked today.

Nice! Can't wait to see them!

santisoler · 2020-10-13T20:20:42Z

If the aim is to extract the topographic correction for the observations stations from the prisms gravity forward, we should be careful with the bord effects.

I totally agree, maybe we should change the gallery example so it takes a larger region than the survey area for building the prisms. Thanks for noticing it!

harmonica/forward/prisms_layer.py

santisoler · 2020-11-13T17:40:16Z

Would be nice to use the new make_xarray_grid for creating the prisms layer. We could also wait for fatiando/verde#300 to be merged so we will be able to pass the easting and northing coordinates either as 1d or 2d arrays.
This way we could generate a prisms layer either from a xarray.Dataset or from numpy arrays created through the vd.grid_coordinates without any hassle like slicing 2d arrays.

santisoler · 2021-03-23T20:29:14Z

@leouieda I'm trying to finish this PR. I think it's almost ready, but I'm using the prisms_layer feature for the Transform21 tutorial and I found that it can create some problems when the surface array has nans.

I decided to be conservative and raise and error if any nan is present on surface or reference.

Another option is to have a drop_nans=False optional argument that replaces the nans on the surface for the respective value of the reference, thus generating a prisms without volume.

What do you think?

This reverts commit be344ca.

The method builds the list of prism boundaries and ignores the ones that has nans on their top or bottom boundaries before passing them into the prism_gravity() function. The density of the prisms is caught from the "density" data_var present in the xarray.Dataset by default.

santisoler · 2021-03-29T19:04:25Z

@leouieda I'm trying to finish this PR. I think it's almost ready, but I'm using the prisms_layer feature for the Transform21 tutorial and I found that it can create some problems when the surface array has nans.

I decided to be conservative and raise and error if any nan is present on surface or reference.

Another option is to have a drop_nans=False optional argument that replaces the nans on the surface for the respective value of the reference, thus generating a prisms without volume.

What do you think?

We decided to allow surface array with nans and add a gravity method that computes the gravity field generated by the layer, ignoring the masked points.

Had to disable the proteceted-access on the entire test_prisms_layer.py file, otherwise pylint will keep complaining about it, even if I disable it function wise or even line wise.

harmonica/forward/prisms_layer.py

leouieda · 2021-04-06T15:42:58Z

harmonica/forward/prisms_layer.py

+ nonans_bottom = np.logical_not(np.isnan(self._obj.bottom.values))
+ return np.logical_and(nonans_top, nonans_bottom)
+
+ def get_prisms(self):


Suggested change

def get_prisms(self):

def to_prisms(self):

This fits better with the rest of the xarray methods: http://xarray.pydata.org/en/stable/pandas.html

Could keep this private for now and provide it after if desired.

leouieda · 2021-04-06T16:18:39Z

harmonica/forward/prisms_layer.py

+ # have no nans on their top and bottom
+ nonans_boundaries = self._get_nonans_prisms_mask().ravel()
+ # Get density array
+ density = getattr(self._obj, density_name).values.ravel()


Suggested change

density = getattr(self._obj, density_name).values.ravel()

density = self._obj[density_name].values.ravel()

Works for density names that aren't valid Python variable names.

Write a private method that returns a mask for every prism whose top boundary, bottom boundary or density is nan. Raise a warning if the density array has a nan where neither the top or bottom is nan. Rename the get_prisms() method to _to_prisms() and make it private. Update and improve the tests.

Run the check when defining the prism layer and when computing the spacing of the layer.

santisoler · 2021-04-07T20:05:08Z

I think this is ready to go. I'll merge it even the CIs are failing: the docs build is failing because it can't download the new South Africa topography dataset. If master fails after merge, I'll fix it right away.

santisoler added 8 commits September 2, 2020 18:11

Start adding prisms_layer function and class

a7e3ec1

Fix a few bugs and extend docstrings

3f8e031

Add prisms_layer and PrismsLayer to doc/api/index.rst

3bec013

Remove iteration methods from PrismsLayer

814abb7

Remove private attrs for size and spacing

c171217

Add tests for prisms_layers and PrismsLayer

9eedffb

Fix wrong ordering of prisms on get_prisms()

92c59ac

Add example for computing terrain correction

18816ef

santisoler changed the title ~~WIP Add funciton to create a layer of prisms~~ WIP Add function to create a layer of prisms Sep 21, 2020

Merge branch 'master' into prisms-layer

0289a2f

Merge branch 'master' into prisms-layer

41518a4

santisoler changed the title ~~WIP Add function to create a layer of prisms~~ Add function to create a layer of prisms Oct 1, 2020

lperozzi reviewed Oct 13, 2020

View reviewed changes

santisoler mentioned this pull request Nov 5, 2020

3D modeling and inversion of the density of a mountain #202

Closed

santisoler commented Nov 10, 2020

View reviewed changes

harmonica/forward/prisms_layer.py Outdated Show resolved Hide resolved

santisoler added 3 commits November 11, 2020 09:30

Merge branch 'master' into prisms-layer

4520a62

Add new datasample for South Africa topography

b810809

Start redisigning prisms_layer function

c515281

santisoler added 7 commits December 16, 2020 14:26

Rename DatasetAccessorPrismsLayer on doc/api/index.rst

1708119

Update gallery examples to the new prisms layer implementation

1c8a173

Update example in docstring of prisms_layer

4226fad

Minor correction on docstring

f94d5e8

Minor fix on docstring example

e6be351

Update tests for prisms_layer

03c11d9

Merge branch 'master' into prisms-layer

ef1ca4b

santisoler requested a review from leouieda December 21, 2020 16:21

Improve text on prism_layer example

4fa5790

santisoler added this to the v0.2.0 milestone Feb 10, 2021

santisoler added 6 commits March 19, 2021 12:22

Merge branch 'master' into prisms-layer

c8ca8fd

Use make_xarray_grid when building the prism_layer

5c17f5f

Set data and data_names as None if no properties are passed

d2a7b9c

Raise errors if nans present in surface or reference

be344ca

Run make format

a625b17

Merge branch 'master' into prisms-layer

9f9cf49

santisoler added 3 commits March 26, 2021 16:59

Revert "Raise errors if nans present in surface or reference"

1fb8ce4

This reverts commit be344ca.

Merge branch 'master' into prisms-layer

bf1d688

santisoler added 2 commits March 29, 2021 16:17

Fix pylint complains

af739bc

Disable pylint protected-access globally o test_prisms_layer

6ed8d9a

Had to disable the proteceted-access on the entire test_prisms_layer.py file, otherwise pylint will keep complaining about it, even if I disable it function wise or even line wise.

santisoler commented Apr 5, 2021

View reviewed changes

harmonica/forward/prisms_layer.py Outdated Show resolved Hide resolved

Fix docstring

26a2775

leouieda reviewed Apr 6, 2021

View reviewed changes

santisoler added 7 commits April 6, 2021 15:51

Rename test function

fe93dcb

Check if easting and northing define a regular grid

0582180

Run the check when defining the prism layer and when computing the spacing of the layer.

Fix import order on test_prisms_layer.py

9a35ca7

Merge branch 'master' into prisms-layer

acca8e5

Merge branch 'master' into prisms-layer

9db7652

Improve examples and fix the one for topo effect

cf43e44

santisoler merged commit 69c5a64 into master Apr 7, 2021

santisoler deleted the prisms-layer branch April 7, 2021 20:06

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add function to create a layer of prisms #186

Add function to create a layer of prisms #186

santisoler commented Sep 2, 2020 •

edited

Loading

santisoler commented Sep 23, 2020

santisoler commented Sep 30, 2020

lperozzi left a comment

lperozzi Oct 13, 2020

santisoler Oct 13, 2020

santisoler commented Oct 13, 2020

santisoler commented Oct 13, 2020

santisoler commented Nov 13, 2020 •

edited

Loading

santisoler commented Mar 23, 2021 •

edited

Loading

santisoler commented Mar 29, 2021

leouieda Apr 6, 2021

leouieda Apr 6, 2021

leouieda Apr 6, 2021

santisoler commented Apr 7, 2021

	density = getattr(self._obj, density_name).values.ravel()
	density = self._obj[density_name].values.ravel()

Add function to create a layer of prisms #186

Add function to create a layer of prisms #186

Conversation

santisoler commented Sep 2, 2020 • edited Loading

santisoler commented Sep 23, 2020

santisoler commented Sep 30, 2020

lperozzi left a comment

Choose a reason for hiding this comment

lperozzi Oct 13, 2020

Choose a reason for hiding this comment

santisoler Oct 13, 2020

Choose a reason for hiding this comment

santisoler commented Oct 13, 2020

santisoler commented Oct 13, 2020

santisoler commented Nov 13, 2020 • edited Loading

santisoler commented Mar 23, 2021 • edited Loading

santisoler commented Mar 29, 2021

leouieda Apr 6, 2021

Choose a reason for hiding this comment

leouieda Apr 6, 2021

Choose a reason for hiding this comment

leouieda Apr 6, 2021

Choose a reason for hiding this comment

santisoler commented Apr 7, 2021

santisoler commented Sep 2, 2020 •

edited

Loading

santisoler commented Nov 13, 2020 •

edited

Loading

santisoler commented Mar 23, 2021 •

edited

Loading