Add support for n-dimensional coordinates on BlockReduce

[santisoler]

Include additional coordinates when reducing data and observation points.

Fixes #195

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

[santisoler]

@leouieda I think this solves the issue, but please take a look at it because I'm not so familiar with Verde source code (it's highly probable that I missed something).

[leouieda]

Hi @santisoler this looks really nice! I made a few comments on the wording of the docs but the code looks good to me. I had to figure out how the original code worked again (of course I didn't remember).

One thing I'm wondering is if we should make this optional and turned off by default. This changes the behavior of the class and might break existing code that's relying on BlockReduce to drop coordinates like lon, lat = BlockReduce(...).filter(coordinates, data).

This could be done with a flag drop_coords=True (or something like that) to BlockReduce. The code could stay the same but drop the coordinates from the tuple before doing the reductions if it's true.

[santisoler]

Thanks @leouieda!

I see the problem by returning n coordinates when only two were expected before this PR.
Although I'm not a big fan of adding so many flags that configure this little things. But I know Verde is already on > 1.0.0 and we should ensure backward compatibility.

I leave the decision to you. Let me know what you decide and I'll make any change you propose.

[leouieda]

Yeah, this would have been good as a default if I had the foresight. But for now, we shouldn't break the API. We can put a warning later on if we decide to change the default to True for a 2.0.0 release.

[santisoler]

@leouieda I gave this a second thought. The problematic line you mentioned is not quite right. BlockReduce.filter() currently returns two tuples: blocked_coordinates and blocked_data.
The common usage is:

coordinates, data = BlockReduce(...).filter(coordinates, data)

So the number of output variables is not changed on this PR, therefore it's not breaking backward compatibility.

Nevertheless it could create problems if the user use it like this:

(lon, lat), data = BlockReduce(...).filter(coordinates, data)

But I don't think this is the recommended way according to the documentation.

I'm on a grey area here (mainly due to my lack of experience with supporting backward compatibility).
What do you think? Should we still add the drop_coords flag?

[leouieda]

I'm on a grey area here (mainly due to my lack of experience with supporting backward compatibility). What do you think? Should we still add the drop_coords flag?

Yep, I always forget that filter returns 3 things always. But it the docstring says that it's only the first 2 coordinates. So even if there are no users relying on this, I don't want to break the promises made in the docstring. I don't think it's a usability problem if we have the drop_coords argument.

Actually, I probably would never use BlockReduce on gravity and magnetic data this way. These are really sensitive to the vertical dimension so averaging out the height seems like a big introduction of errors. Instead, I would the BlockMean coordinates to set the point masses and fit the entire dataset.

But it's good to have the option of not dropping coordinates since the average height, etc, might be useful for other things.

[leouieda]

Looks good to me! Thanks, Santi!