diff --git a/pygmt/__init__.py b/pygmt/__init__.py index 512865264b0..4a38fe40ad6 100644 --- a/pygmt/__init__.py +++ b/pygmt/__init__.py @@ -19,7 +19,7 @@ from .sampling import grdtrack from .mathops import makecpt from .modules import GMTDataArrayAccessor, config, info, grdinfo, which -from .gridops import grdcut +from .gridops import grdcut, grdfilter from .x2sys import x2sys_init, x2sys_cross from . import datasets diff --git a/pygmt/gridops.py b/pygmt/gridops.py index 122aa42340e..4adca2f6e6b 100644 --- a/pygmt/gridops.py +++ b/pygmt/gridops.py @@ -20,13 +20,13 @@ @fmt_docstring @use_alias( - G="outgrid", - R="region", - J="projection", - N="extend", - S="circ_subregion", - Z="z_subregion", -) + G="outgrid", + R="region", + J="projection", + N="extend", + S="circ_subregion", + Z="z_subregion" + ) @kwargs_to_strings(R="sequence") def grdcut(grid, **kwargs): """ @@ -113,3 +113,85 @@ def grdcut(grid, **kwargs): result = None # if user sets an outgrid, return None return result + +@fmt_docstring +@use_alias( + G="outgrid", + F="filter", + D="distance" + ) +@kwargs_to_strings(R="sequence") + +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 +