plot 2d ndarrays with grdimage

[yohaimagen]

Description of the desired feature

I find myself working with 2D numpy ndarrays often. Then, for visualization purposes, I need to input the data into an xarray DataArray and plot it using pyGMT. I think it would be very convenient if the grdimage function could accept 2D ndarrays along with x and y coordinates. Another option, if PyGMT's grdimage should not change significantly from GMT's grdimage, would be to provide a higher-level function on top of grdimage that does the same. Something that resembles matplotlib's pcolormesh.

Are you willing to help implement and maintain this feature?

Yes

[seisman]

I find myself working with 2D numpy ndarrays often. Then, for visualization purposes, I need to input the data into an xarray DataArray and plot it using pyGMT.

Just want to make sure we're on the same page. Are you writing codes like this?

import numpy as np
import xarray as xr
import pygmt

array = np.arange(200).reshape((10, 20))
data = xr.DataArray(array)
fig = pygmt.Figure()
fig.grdimage(data, frame=True)
fig.show()

[yohaimagen]

Not exactly

A simple example is a calculation for parameters x and y, which can represent longitudes and latitudes, but they can also represent anything else.

import numpy as np
import xarray as xr
import pygmt

x = np.linspace(-118, -117,  100) # longitudes
y = np.linspace(35, 36, 100) # latitudes
X, Y = np.meshgrid(x, y)
array = np.zeros_like(X)
for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        array[i, j] = f(X[i, j], Y[i, j])# f calculate something based on x and y
data_array = xr.DataArray(array, coords={'x': x, 'y': y}, dims=('y', 'x'))

fig = pygmt.Figure()
fig.grdimage(data_array, frame=True, projection='M10')
fig.show()

[yohaimagen]

I think the main point is that it will be nice to be able to plot 2D data without the need to go through raster-like data, whether it be writing it to a raster file or using xarray.

[seisman]

Not exactly

A simple example is a calculation for parameters x and y, which can represent longitudes and latitudes, but they can also represent anything else.

import numpy as np
import xarray as xr
import pygmt

x = np.linspace(-118, -117,  100) # longitudes
y = np.linspace(35, 36, 100) # latitudes
X, Y = np.meshgrid(x, y)
array = np.zeros_like(X)
for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        array[i, j] = f(X[i, j], Y[i, j])# f calculate something based on x and y
data_array = xr.DataArray(array, coords={'x': x, 'y': y}, dims=('y', 'x'))

fig = pygmt.Figure()
fig.grdimage(data_array, frame=True, projection='M10')
fig.show()

We have a similar but simpler example (https://www.pygmt.org/dev/gallery/3d_plots/grdview_surface.html), so your script can be simplified to:

x = np.linspace(-118, -117,  100) # longitudes
y = np.linspace(35, 36, 100) # latitudes
data = xr.DataArray(f(*np.meshgrid(x, y)), coords={'x': x, 'y': y}, dims=('y', 'x'))

[seisman]

I think we have three options here:

1. Do nothing and close the issue
2. Make all grid-related functions (e.g., grdimage/`grdcut``) to accept a 2-D numpy array as input. Internally, we need to (1) check if the input data is a 2-D numpy array; (2) Create a xarray.DataArray object from the 2-D numpy array, with the x/y coords determined from the array shape
3. Provide a function like matrix2grid to convert a 2-D matrix to a xarray.DataArray object. The function's definition can be something like:

   matrix2grid(data, x=None, y=None, region=None, spacing=None)

   For reference, GMT.jl also provides a function named mat2grid at https://www.generic-mapping-tools.org/GMTjl_doc/documentation/utilities/mat2grid/index.html.

I'm inclined to option 1 or 3.