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'))