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.xyz2grd: Allow more flexible ways to pass x/y/z arrays #2852

Open
seisman opened this issue Dec 5, 2023 · 0 comments
Open

pygmt.xyz2grd: Allow more flexible ways to pass x/y/z arrays #2852

seisman opened this issue Dec 5, 2023 · 0 comments
Labels
enhancement Improving an existing feature

Comments

@seisman
Copy link
Member

seisman commented Dec 5, 2023

The pygmt.xyz2grd function accepts 1-D x/y/z arrays as input and produces a grid. Here, x/y are the x- and y-coordinates of grid nodes and z is the data values on grid nodes. Thus, for a 3x4 grid, x/y/z arrays all have 12 values (assuming it's gridline-registrated).

The pygmt.xyz2grd API documentation provides a simple example showing how we should prepare the x/y/z arrays.

import numpy as np
import pygmt

def func(x, y): 
    return x**2 + y**2

x, y = np.meshgrid([0, 1, 2, 3], [10.5, 11.0, 11.5, 12.0, 12.5])
z = func(x, y)
xx, yy, zz = x.flatten(), y.flatten(), z.flatten()
grid = pygmt.xyz2grd(x=xx, y=yy, z=zz, spacing=(1.0, 0.5), region=[0, 3, 10, 13])

As you can see, we have to use np.meshgrid to generate the 2-D x/y arrays, call a function (here z=x**2+y**2) to produce the 2-D z array, and use the flatten() method to convert the 2-D arrays into 1-D arrays. It's not straightforward and sometimes difficult to understand.

Here I propose to provide two more flexible ways to pass x/y/z arrays. The expected syntaxes are:

x = [0, 1, 2, 3]
y = [10.5, 11.0, 11.5, 12.0, 12.5]
z = func(*np.meshgrid(x, y))

# Case 1: 2-D z array + 1-D x/y arrays
grid = pygmt.grid(x=x, y=y, z=z)  # spacing and region are no longer needed

# Case 2: 2-D z array + region + spacing
grid = pygmt.grid(z=z, region=[0, 3, 10.5, 12.5], spacing=(1, 0.5))

Some notes:

  1. Internally we should do something like before passing data to the xyz2grd module:
    x = np.arange(region[0], region[1], spacing[0])  # required for case 2
    y = np.arange(region[2], region[3], spacing[1])  # required for case 2
    xx, yy = np.meshgrid(x, y)
    xx, yy, zz = xx.flatten(), yy.flatten(), z.flatten()
    # now pass xx/yy/zz to xyz2grd
    
  2. The -I option (spacing parameter) is required for xyz2grd. The spacings can be obtained by spacing = (x[1] - x[0], y[1] - y[0]) assuming that x and y are equal-spaced.
  3. For case 2, a special case is, only the 2-D z array is given and region/spacing are not given, then region is default to [0, nx-1, 0, ny-1], and spacing is default to [1, 1]
  4. In all cases, we need to be careful with the registration parameter.

Comments?

@seisman seisman added the enhancement Improving an existing feature label Dec 5, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement Improving an existing feature
Projects
None yet
Development

No branches or pull requests

1 participant