Merge overlapping grids

[ckohnke]

Description of the desired feature:

Grid stitching methods to handle grid discrepancies at overlapping grid points. This may already be possible with the tools in Verde and an example may be beneficial if it is.

Using airborne magnetic data as an example:

• After processing, leveling, and gridding individual flights of a magnetic survey there is usually a difference between grids in the region where they overlap.
• Theoretically, these overlapping points are at the same time (diurnal removed), and space (equivalent source, upward continuation, etc.), and should have the same value. However, a variety of factors, such as instrument drift, can cause mismatch.
• If one is lucky, there is only a static shift between the grids in the overlapping region and can be handled by using a common datum between the datasets.
• However, if the difference in the data in the overlapping regions are more complex, stitching the grids together using the overlapping points becomes more nebulous and can easily lead to errors in the resulting final grid.

Are you willing to help implement and maintain this feature?

Potentially. I've been having trouble finding mathematical methods and processes to handle such situations.

[leouieda]

@ckohnke thanks for opening the issue! I have a few thoughts on this:

1. For gravity and magnetics, I really don't like these stitching approaches because they can heavily bias the long wavelengths. For those, I'd rather go back to the original data and use a single equivalent source model to re-grid them all.
2. You're right about there not being a good way to do this. For grids of other data types (topography, etc), maybe we could try to blend them in some way that only introduces distortions in the boundaries and doesn't propagate them to the whole grid. Or maybe warn users if there are a lot of discrepancies.
3. GMT's grdblend does something like this already. It's not yet available in PyGMT. We could try to implement something similar to what it does using xarray. It doesn't seem like it would be too difficult to do.

[leouieda]

It would also be worth looking into what xarray has for merging datasets. Perhaps that could be enough and we don't have to implement it ourselves.

[mdtanker]

I've just had a look at this and it seems xarray.merge for overlapping coordinates of grids allows:

1. just keeping the value from the first grid: compat="override"
2. only keeping values which are NaNs in 1 grid but not the others: compat="no_conflicts
3. the remaining options seems to just omit overlapping values.

There doesn't seem to be an option for "blending" the overlapping values.

One solution might be to make any overlapping values nans with compat="no_conflicts", then interpolate the resulting NaNs, for example with the soon-to-be-implemented vd.fill_nans function

Rioxarray has a rio.merge.merge_arrays function with a method parameter:

1. you can create two grids with the counts (method="count"), and sums (method="sum") of the overlapping cells, then divide them to get a grid where the overlapping cells are the mean of the input grids.
2. you can pass a custom function to method, which could calculate the mean of the supplied grids, as shown in this stackoverflow question.

[leouieda]

Alright, seems like xarray doesn't have this. I also think it's usually very tricky to code for xarray since everything has to work for n-dimensions. I'm in favor of trying to add something to Verde for this. The main issue is how to blend the grids. GMT does a weighted average on overlapping regions and has lots of options. We could start with a simpler approach of just taking averages or trying @mdtanker's suggestion of interpolating.

[mdtanker]

I think starting with a weighted average (which defaults to equal weights) would be a sensible place to start and should be straightforward to implement. I'm happy to draft something.

[mdtanker]

Check out Georasters as well for inspiration.