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]

👋 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:

pygmt/pygmt/gridops.py

Lines 21 to 115 in d3e131a

	@fmt_docstring
	@use_alias(
	G="outgrid",
	R="region",
	J="projection",
	N="extend",
	S="circ_subregion",
	Z="z_subregion",
	)
	@kwargs_to_strings(R="sequence")
	def grdcut(grid, **kwargs):
	"""
	Extract subregion from a grid.
	Produce a new *outgrid* file which is a subregion of *grid*. The
	subregion is specified with *region*; the specified range must not exceed
	the range of *grid* (but see *extend*). If in doubt, run
	:meth:`pygmt.grdinfo` to check range. Alternatively, define the subregion
	indirectly via a range check on the node values or via distances from a
	given point. Finally, you can give *projection* for oblique projections to
	determine the corresponding rectangular *region* setting that will give a
	grid that fully covers the oblique domain.
	Full option list at :gmt-docs:`grdcut.html`
	{aliases}
	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.
	{J}
	{R}
	extend : bool or int or float
	Allow grid to be extended if new *region* exceeds existing boundaries.
	Give a value to initialize nodes outside current region.
	circ_subregion : str
	``'lon/lat/radius[unit][+n]'``.
	Specify an origin (*lon* and *lat*) and *radius*; append a distance
	*unit* and we determine the corresponding rectangular region so that
	all grid nodes on or inside the circle are contained in the subset.
	If **+n** is appended we set all nodes outside the circle to NaN.
	z_subregion : str
	``'[min/max][+n|N|r]'``.
	Determine a new rectangular region so that all nodes outside this
	region are also outside the given z-range [-inf/+inf]. To indicate no
	limit on *min* or *max* only, specify a hyphen (-). Normally, any NaNs
	encountered are simply skipped and not considered in the
	range-decision. Append **+n** to consider a NaN to be outside the given
	z-range. This means the new subset will be NaN-free. Alternatively,
	append **+r** to consider NaNs to be within the data range. In this
	case we stop shrinking the boundaries once a NaN is found [Default
	simply skips NaNs when making the range decision]. Finally, if your
	core subset grid is surrounded by rows and/or columns that are all
	NaNs, append **+N** to strip off such columns before (optionally)
	considering the range of the core subset for further reduction of the
	area.
	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("grdcut", 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

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 .