WIP Add tesseroids layer

[santisoler]

Add function to create a teseroids layer as a Dataset, which contains the longitudinal and latitudinal coordinates of the center point of each tesseroid, along with a few data variables containing the top and bottom coordinates and the density of each tesseroid.

The top and bottom parameters can be passed in either a geocentric spherical or a geodetic coordinate system.
The region, spacing and shape are handled in a similar way as Verde's grid_coordinates does.

Fixes #83

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.
• Write detailed docstrings for all functions/methods.
• If adding new functionality, add an example to the docstring, gallery, and/or tutorials.

[leouieda]

Hi Santi, thanks for taking a stab at this!

1. Is there a reason for choosing to store the center of the tesseroid instead of the boundary coordinates?
2. An advantage for storing the boundaries is that it enables non-uniform meshes (with padding in the edges, for example).
3. Now that I think about, I don't think you can have coordinates that don't match the data shape (which would happen with storing boudaries). But maybe we could the boundaries as auxiliary coordinates or data? When generating actual tesseroids it would be best to have the boundaries. The cell size would also be OK.
4. I see the application you have in mind for the region_centers argument but it might case problems if you forget that you set it to True and then the region doesn't match the data. Could we do the same with verde.pad_region before passing it to this function?

There are also a few comments inline.

This is looking pretty good! We can use something like this to make things like topography_to_tesseroids etc.

[leouieda]

Is there a reason for repeating all of these lines instead of using verde.grid_coordinates?

[santisoler]

I wanted to take advantage of xarray coordinates: we don't need to make np.meshgrid(longitude, latitude), so I though we can save a little bit of memory by omitting it.

Nevertheless I was just drafting, so we can still think how to reduce code repetition.

[leouieda]

cell_size is probably a better term than spacing for meshes

[santisoler]

I kept spacing because the adjust argument of verde.coordinates.spacing_to_shape must be either shape or spacing.
I don't quite like it for the tesseroids, so we could change the argument name and internally pass spacing to spacing_to_shape when adjust="spacing".

PS: I want the legal values for adjust and the cell_size/spacing argument are in agreement.

[leouieda]

Can't these be arrays as well?

[leouieda]

What about the latitudes? It would make more sense for everything to be in geodetic, not just top and bottom. It's OK if the grid is not evenly spaced in spherical latitude if you ask for geodetic coordinates.

[santisoler]

I left the latitudes unchanged because I want to have evenly spaced grids on spherical coordinates.
Should we handle this differently?

Maybe the docstring is not clear enough, but the tesseroid_layer will always return the xarray in spherical coordinates. So when I pass coordinates="geodetic" it means that the top and bottom are in geodetic coordinates, but the returned xarray will be in spherical.

[leouieda]

What if instead of this, we take projection = None or function. In this case, projection would be any function that transforms coordinates into a spherical geocentric system.

[santisoler]

I like that!
With this in mind, maybe we could create a private function to create layers or meshes and then call it for each specific geometry (tesseorid, prisms, etc).
So, when using prisms, the projection could be one from pyproj, while on the tesseroid case, it could be geodetic_to_spherical.

[santisoler]

Hi Santi, thanks for taking a stab at this!

1. Is there a reason for choosing to store the center of the tesseroid instead of the boundary coordinates?

2. An advantage for storing the boundaries is that it enables non-uniform meshes (with padding in the edges, for example).

3. Now that I think about, I don't think you can have coordinates that don't match the data shape (which would happen with storing boudaries). But maybe we could the boundaries as auxiliary coordinates or data? When generating actual tesseroids it would be best to have the boundaries. The cell size would also be OK.

The main reason to store tesseroid centers is, as you pointed out, because I wanted to have data arrays with one value for each tesseroid.

I assume that by non-uniform meshes you mean not equally spaced, right?
We can still generate non-uniform meshes under this schema: you'll get the tesseroid boundary from the next neighbor:

east = longitude[i] + (longitude[i + 1] - longitude[i]) / 2
west = longitude[i] - (longitude[i] - longitude[i - 1]) / 2

That way we also ensure that there are no gaps inside the layer.

Storing the cell size is also a good option.
But we should be careful because if we or a user accidentally mess with this, could generate gaps or overlap between tesseroids.

In summary, I would store only independent information about the geometry of the tesseroids.
So, if we don't want gaps or overlapping, storing their centers is enough because the boundaries can be get from them by the no-gaps-no-overlapping assumption.

4. I see the application you have in mind for the `region_centers` argument but it might case problems if you forget that you set it to `True` and then the region doesn't match the data. Could we do the same with `verde.pad_region` before passing it to this function?

I'm a little bit confused about this. What do you mean by the region doesn't match the data?

Now you commented the region_centers argument, I don't know if it's absolutely needed or if it will cause more headaches than the ones that will solve.
I'm thinking that an obvious application of tesseroid_layer is to generate a tesseroid layer that mimics some relief data (topography, bathymetry, moho depth, ice bed, etc). Therefore the function should create a tesseroid for each grid point, locating its center on the coordinate of the grid point (i.e. region_center=True). This also prevents having computation points on tesseroids vertices in case that we use a computation grid located on the same longitude, latitude coordinates as the relief (e.g. for inversions that simplify the jacobian matrix).
I cannot see an application for the region_center=False case.

There are also a few comments inline.

This is looking pretty good! We can use something like this to make things like topography_to_tesseroids etc.

Thanks!

[santisoler]

I'm closing this PR. We should base our tesseroid layer on the prism layer already implemented in Harmonica. I think it worth start writing the tesseroid implementation from scratch.