Add neighborhood filter to remap

[ahijevyc]

Apply a neighborhood filter within a circular radius r (in degrees) to a UxDataset or UxDataArray.

Issue #930

Overview

This neighborhood filter is written under remap so it has the ability to remap, but perhaps the neighborhood filter could be moved up to its own function, like UxDataset.neighborhold_filter. One doesn't always need to remap. As a workaround, if you don't need to remap, just supply the same destination UxGrid as the source.

This is kind of like uxarray.UxDataArray.inverse_distance_weighted_remap , but the neighborhood is defined by distance, not a number of nearest neighbors. This is ideally suited for a variable resolution mesh, in which a constant of neighbors doesn't have a constant sized neighborhood. Another difference is that this neighborhood filter does not weight data by inverse distance.

This function shares the same logic as uxarray.UxDataArray.subset.bounding_circle to select grid elements in a circular neighborhood, but this function applies it to all elements in grid, not just one center_coordinate.

The filter function func may be a user-defined function, but uses np.mean by default. It could be min, max, np.median. It can even use np.percentile if you supply a single quantile argument with functools.partial (see below)

func = np.percentile
from functools import partial
func = partial(func, q=90)

The new function remap.apply_func.py used remap.inverse_distance_weighted.py and remap.nearest_neighbor.py as a starting point. inverse_distance_weighted.py and nearest_neighbor.py have a lot of the same code but have different documentation style. Their formatting styles should probably be brought back into agreement before this pull request is completed.

Expected Usage

import numpy as np

import uxarray

grid_path = "/glade/campaign/mmm/wmr/weiwang/cps/irma3/2020/tk707_conus/init.nc"
data_path = "/glade/campaign/mmm/wmr/weiwang/cps/irma3/mp6/tk707/diag.2017-09-07_09.00.00.nc"
uxds = uxarray.open_mfdataset(
    grid_path,
    data_path
)

# Trim domain
lon_bounds = (-74, -64)
lat_bounds = (18, 24)
uxda = uxds["refl10cm_max"].isel(Time=0).subset.bounding_box(lon_bounds, lat_bounds)
uxda


# this is how you use this function to smooth with 0.25-deg filter.
uxda_smoothed_mean = uxda.remap.apply_func(uxda.uxgrid, func=np.mean, r=0.25)


# this is another way to use this function with a maximum filter
uxda_smoothed_max = uxda.remap.apply_func(uxda.uxgrid, func=max, r=0.25)

PR Checklist

General

• An issue is linked created and linked
• Add appropriate labels
• Filled out Overview and Expected Usage (if applicable) sections

Testing

• Adequate tests are created if there is new functionality
• Tests cover all possible logical paths in your function
• Tests are not too basic (such as simply calling a function and nothing else)

Documentation

• Docstrings have been added to all new functions
• Docstrings have updated with any function changes
• Internal functions have a preceding underscore (_) and have been added to docs/internal_api/index.rst
• User functions have been added to docs/user_api/index.rst

Examples

• Any new notebook examples added to docs/examples/ folder
• Clear the output of all cells before committing
• New notebook files added to docs/examples.rst toctree
• New notebook files added to new entry in docs/gallery.yml with appropriate thumbnail photo in docs/_static/thumbnails/

[philipc2]

Hi @ahijevyc

This is a great contribution. Can you create an issue & link it to this PR (can restate parts from the description). I'll give this a detailed review after this weekend, however a couple initial thoughts:

We could have two forms of this function:

• UxDataArray.remap.neighborhood_filter() would act the same as UxDataArray.remap.apply_func() that you've implemented.
• UxDataArray.apply_neighborhood_filter() could apply the filter to the same Grid

This reminds me of our Topological Aggregations, except instead of using the connectivity to determine the "neighborhood", it uses the number of neighbors returned from a nearest neighbor query.

Could another name for this type of functionality be referred to as a "Spatial Aggregation"? A sample taking from another discussion:

"spatial" average should refer in my opinion to some sort of distance from the node to the center of the faces that fall within a radius; again, some sort of weighting based on the distance should be involved ( RBF is radial basis function , etc)

#724 (comment)

The new function remap.apply_func.py used remap.inverse_distance_weighted.py and remap.nearest_neighbor.py as a starting point. inverse_distance_weighted.py and nearest_neighbor.py have a lot of the same code but have different documentation style. Their formatting styles should probably be brought back into agreement before this pull request is completed.

Noted! I'll take a look into it.

[ahijevyc]

Thanks. And yes, I will create an issue and link it to this PR. I like the idea of - UxDataArray.apply_neighborhood_filter() which applies the filter to the same Grid Would we lose any functionality by completely separating the remapping capability from the neighborhood filter? Instead of using a function UxDataArray.remap.neighborhood_filter(), one would chain remap after UxDataArray.apply_neighborhood_filter(), like UxDataArray.apply_neighborhood_filter().remap.nearest_neighbor() ? This would require an extensive revision of my pull request, but I don't mind if it simplifies the code. By the way, should we use this call structure with the remap accessor UxDataArray.remap.nearest_neighbor() or this call structure UxDataArray.nearest_neighbor_remap()? UxDataset.nearest_neighbor_remap <https://uxarray.readthedocs.io/en/latest/user_api/generated/uxarray.UxDataset.nearest_neighbor_remap.html#uxarray.UxDataset.nearest_neighbor_remap> is in the User API, but UxDataArray.remap.nearest_neighbor() is in the User Guide on Remapping. To your comment about "spatial Aggregation", apply_neighborhood_filter()is definitely a *spatial* method, with its reliance on radius=r, but it doesn't *aggregate* in the sense of reducing the number of elements or dimensions. In other words, the output has the same number of dimensions and elements as the input. And finally, you are right; we should probably add the capability to weight neighborhood elements by inverse distance. I originally just wanted a "maximum" filter, so it didn't require weighting. But if you want a "mean" filter, weighting the elements on inverse distance, or better yet, area, makes sense.

[aaronzedwick]

"By the way, should we use this call structure with the remap accessor
UxDataArray.remap.nearest_neighbor()
or this call structure
UxDataArray.nearest_neighbor_remap()?"

To comment on this, you should use UxDataArray.remap.nearest_neighbor(). You can directly call UxDataArray.nearest_neighbor_remap(), but it will give a deprecation warning, and actually will not work anymore since we added some updates to remap, which weren't added to UxDataArray.nearest_neighbor_remap() as it was supposed to be removed awhile ago. Thanks for bringing this up! The deprecation warning has been there for several months so I think it is safe to delete this code altogether. I will put together a pr for this.

[ahijevyc]

I am preparing a new pull request #941 with the neighborhood filter on its own. It is not a part of remap.