Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

pygmt.grdfill: Interpolate across holes in a grid #3768

Closed
8 of 11 tasks
seisman opened this issue Jan 12, 2025 · 2 comments
Closed
8 of 11 tasks

pygmt.grdfill: Interpolate across holes in a grid #3768

seisman opened this issue Jan 12, 2025 · 2 comments
Labels
feature request New feature wanted
Milestone

Comments

@seisman
Copy link
Member

seisman commented Jan 12, 2025

This is the central issue to track the implementation of the pygmt.grdfill method.

GMT options

Checked: Implemented; Unchecked: To be implemented/discussed; Strikethrough: Won't implement.

@seisman seisman converted this from a draft issue Jan 12, 2025
@seisman seisman added the feature request New feature wanted label Feb 21, 2025
@seisman
Copy link
Member Author

seisman commented Mar 11, 2025

Here are my thoughts on finalizing the wrapper.

-N

This option has different long names across the GMT ecosystem.

  • PyGMT: aliased to no_data. [Need to note that common option -d is aliased to nodata]
  • GMT.jl: aliased to nodata or hole
  • GMT: aliased to hole or hole_value

This option "sets the node value used to identify a point as a member of a hole [Default is NaN].". Personally, I feel hole is a better name than no_data.

-A

This option specifie the hole-filling algorithm to use. Currently it's aliased to mode and can be used like (please refer to the docs for the actual meanings of these arguments):

  • mode="c100"
  • mode="ganother_grid.nc"
  • mode="n0.5"
  • mode="s0.3"

The syntax is not Pythonic, since we have to specify numbers in a string, and for the mode="ganother_grid.nc" case, the "another grid" must be a string, and can't accept a xarray.DataArray object.

In my mind, there are two possible ways to make it more Pythonic:

  1. Two parameters mode and value.

    • mode specifies the filling method. Valid values are "constant", "grid", "neighbor", "spline"
    • value has different meanings for different filling methods, e.g., value=100, value="another_grid.nc", value=0.5, value=0.3
  2. Four exclusive parameters

    • constantfill=100
    • gridfill="anonter_grid.nc" or gridfill=xrgrid
    • neighborfill=True or neighborfill=0.5: 0.5 is the search radius, and True means the default radius
    • splinefill=True/splinefill=0.3: 0.3 is the spline tension, True means no tension.

I prefer the 2nd option.

What do others think?

@seisman
Copy link
Member Author

seisman commented Mar 27, 2025

The -L option is still not implemented yet. Now let's focus on this option.

Behaviors in GMT

-L
Just list the rectangular subregions west east south north of each hole. No grid fill takes place and -G is ignored. Optionally, append p to instead write closed polygons for all subregions.

So, there are two possible usage for this option:

  • -L: output the bounds (in the form of (w, e, s, n)) of the holes
  • -Lp: output the closed polygons (rectangles) based on the bounds (not the exact polygons of holes)

Here are examples to understand the outputs:

$ gmt grdfill @earth_relief_20m_holes.grd -L
1.83333333333	6.16666666667	3.83333333333	8.16666666667
6.16666666667	7.83333333333	0.5	2.5
$ gmt grdfill @earth_relief_20m_holes.grd -Lp
>
1.83333333333	3.83333333333
6.16666666667	3.83333333333
6.16666666667	8.16666666667
1.83333333333	8.16666666667
1.83333333333	3.83333333333
>
6.16666666667	0.5
7.83333333333	0.5
7.83333333333	2.5
6.16666666667	2.5
6.16666666667	0.5

Need to note that -Lp write the coordinates of based on the bounds, not the exact outlines of holes:

gmt begin map
gmt grdimage @earth_relief_20m_holes.grd -Baf
gmt grdfill @earth_relief_20m_holes.grd -Lp | gmt plot -W1p
gmt end show

Image

The PyGMT syntax

  • For -L, the outputs are bounds in the form of (w, e, s, n), which fit well with a 2-D numpy array.
  • For -Lp, the outputs are multi-segment files, which PyGMT doesn't support well yet.

Since the closed polygons can be created from the bounds, I think -Lp has a low priority and we don't even have to implement it.

Then, we need to decide the parameter name for -L. GMT and GMT.jl alias -L to list, but list is a bad name in the Python world. So, what about dump or inquire (xref: #896)?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
feature request New feature wanted
Projects
Status: 4) Gallery|tutorial example
Development

No branches or pull requests

1 participant