Wrap grdfilter

[carocamargo]

Hello, I would like to know if there are plans to incorporate the grdfilter function from GMT in pygmt? http://gmt.soest.hawaii.edu/doc/latest/grdfilter.html#examples

The function is used to apply different filters in space (or time) domain.

Thanks

[welcome[bot]]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. You might also want to take a look at our contributing guidelines and code of conduct.

[weiji14]

Hi @carocamargo, no plans yet for grdfilter but you're more than welcome to give it a go. A quick look at the documentation suggests that it could implemented similarly to pygmt.grdcut which also takes an input grid and outputs another grid, something like so:

https://github.com/GenericMappingTools/pygmt/blob/d3e131a863a008f14979e7d9764cc6167156376f/pygmt/gridops.py#L21-L115

Do you have an example on how grdfilter works?

[carocamargo]

I'm not sure of what happens behind the function. I use it a bit like a blackbox.. But you provide an input grid (say a world map, 1˚ resolution), and then you choose the type of the filter (e.g. boxcar, gaussian, median, ...), and the width of the filter (e.g. 600km filter), and then you just need to add a flag stating if the distance was in km, or in units for example.

What happens inside, I dont really know (as I said, it's a blackbox to me).. But you get a smoother field than your input.. I use to get rid of the 'ripple' effects from the harmonic solution (see the image example. On the left is the input file, and on the right the file after using grdfilter)..

\$ gmt gridfilter input.nc -Gfilterted_output.nc -Fm2000 -D4 -fg

[weiji14]

That's a cool plot, amazing how much changing ice volume can affect sea level around Antarctica!

Anyways, what you could try is to take the grdcut code linked above, and simply swap grdcut with grdfilter, and I wouldn't be surprised if it just works. Maybe change the docstring too. It's getting late on my timezone but I can take a look tomorrow if I have time.

[carocamargo]

Thanks! It worked modifying a bit the grdcut! So I just added the new args in the alias list, and changed the name of the function in call_module

@use_alias( G="outgrid", R="region", J="projection", N="extend", S="circ_subregion", Z="z_subregion", F="filter", D="distance" )

def grdfilter(grid, **kwargs): """

    filter a grid file in the time domain using one of the selected convolution
    or non-convolution isotropic or rectangular filters and compute distances
    using Cartesian or Spherical geometries. The output grid file can optionally
    be generated as a sub-region of the input (via -R) and/or with new increment
    (via -I) or registration (via -T). In this way, one may have “extra space” in
    the input data so that the edges will not be used and the output can be within
    one half-width of the input edges. If the filter is low-pass, then the output
    may be less frequently sampled than the input.

    Parameters
    ----------
    grid : str or xarray.DataArray
    The file name of the input grid or the grid loaded as a DataArray.
    outgrid : str or None
    The name of the output netCDF file with extension .nc to store the grid
    in.
    {F} : str
    Name of filter type you which to apply, followed by the width
    b: Box Car; c: Cosine Arch; g: Gaussian; o: Operator; m: Median; p: Maximum Likelihood probability; h: histogram
    {D}: str
    Distance flag, that tells how grid (x,y) rrlated to the filter width as follows:
    flag = p: grid (px,py) with width an odd number of pixels; Cartesian distances.

    flag = 0: grid (x,y) same units as width, Cartesian distances.

    flag = 1: grid (x,y) in degrees, width in kilometers, Cartesian distances.

    flag = 2: grid (x,y) in degrees, width in km, dx scaled by cos(middle y), Cartesian distances.

    The above options are fastest because they allow weight matrix to be computed only once. The next three options are slower because they recompute weights for each latitude.

    flag = 3: grid (x,y) in degrees, width in km, dx scaled by cosine(y), Cartesian distance calculation.

    flag = 4: grid (x,y) in degrees, width in km, Spherical distance calculation.

    flag = 5: grid (x,y) in Mercator -Jm1 img units, width in km, Spherical distance calculation.

    Returns
    -------
    ret: xarray.DataArray or None
    Return type depends on whether the *outgrid* parameter is set:
    - xarray.DataArray if *outgrid* is not set
    - None if *outgrid* is set (grid output will be stored in *outgrid*)
    """
kind = data_kind(grid)

with GMTTempFile(suffix=".nc") as tmpfile:
    with Session() as lib:
        if kind == "file":
            file_context = dummy_context(grid)
        elif kind == "grid":
            file_context = lib.virtualfile_from_grid(grid)
        else:
            raise GMTInvalidInput("Unrecognized data type: {}".format(type(grid)))

        with file_context as infile:
            if "G" not in kwargs.keys():  # if outgrid is unset, output to tempfile
                kwargs.update({"G": tmpfile.name})
            outgrid = kwargs["G"]
            arg_str = " ".join([infile, build_arg_string(kwargs)])
            lib.call_module("grdfilter", arg_str)

    if outgrid == tmpfile.name:  # if user did not set outgrid, return DataArray
        with xr.open_dataarray(outgrid) as dataarray:
            result = dataarray.load()
            _ = result.gmt  # load GMTDataArray accessor information
    else:
        result = None  # if user sets an outgrid, return None

    return result

[carocamargo]

out=pygmt.grdfilter('/Users/Desktop/ANT_trend.nc',F='m1600',D='4', G='/Users/Desktop/ANT_trend_pygmt.nc')

[weiji14]

Great! Do you want to try to start a Pull Request for this, since I saw you forked the repo already? I can help polish it up if you don't have time, but making a start shouldn't be too hard. Just modify the same gridops.py file with this new grdfilter function.

[carocamargo]

Just did (I think...It's pending the request). I just modified the and also the <__init__.py> to call also from .