plot xarray Data array in pyGMT

Question

I have a grid (xyz) that covers the entire globe. The dataset consists of multiple columns containing, longitude, latitude, and several z-values. Here are the locations of each node:

Now I want to show any of the z-values using fig.grdimage. However, I am unable to construct a grid from my xyz-file.

I first read in the data using pandas. If I try to convert the dataframe to an xarray data array and plot it, I get the following error:

array = data.to_xarray()
fig = pygmt.Figure()
fig.basemap(region="d", projection="N12c", frame=True)
fig.grdimage(grid=array)
fig.coast(shorelines="0.5p,black")
fig.show()

"array" gives the following output:

Dimensions:    (index: 16200)
Coordinates:
  * index      (index) int64 0 1 2 3 4 5 ... 16194 16195 16196 16197 16198 16199
Data variables: (12/14)
    LONG       (index) float64 -179.0 -179.0 -179.0 -179.0 ... 179.0 179.0 179.0
    LAT        (index) float64 -89.0 -87.0 -85.0 -83.0 ... 83.0 85.0 87.0 89.0
    ELEVATION  (index) float64 2.965e+03 2.915e+03 ... -3.822e+03 -3.392e+03
    MOHO       (index) float64 -3.888e+04 -3.888e+04 ... -1.457e+04 -1.457e+04
    LAB        (index) float64 -1.2e+05 -1.2e+05 ... -1.245e+05 -1.245e+05
    RHO_C      (index) float64 2.851e+03 2.851e+03 ... 2.806e+03 2.806e+03
    ...         ...
    BOTTOM     (index) float64 -4.1e+05 -4.1e+05 -4.1e+05 ... -4.1e+05 -4.1e+05
    GEOID      (index) float64 -0.81 5.01 4.78 0.15 ... 3.28 3.42 1.52 -0.73
    FA         (index) float64 0.82 38.25 36.06 -17.35 ... 6.77 13.84 3.85 -2.68
    G_zz       (index) float64 0.01 0.33 0.22 -0.08 ... 0.07 0.12 0.04 -0.0
    G_xx       (index) float64 0.04 -0.1 -0.1 -0.02 ... -0.12 -0.1 -0.06 -0.04
    G_yy       (index) float64 -0.12 -0.23 -0.12 0.1 ... 0.05 -0.02 0.02 0.04
16200 16200 16200

However, I don't know how to specify which z-value of array I would like to show. The result is that GMT does not recognise the data type for grid.

I have extract the lon, lat and z-value columns and I tried to use pygmt.xyz2grd and pygmt.sphinterpolate to construct a grid, but that doesn't seem to work either. I'm convinced this should be very simple, but the pyGMT documentation (nor ChatGPT) provide a straightforward solution.

I have also tried using matplotlib's basemap (different projection):

map = Basemap(projection='ortho',lat_0=15,lon_0=160,resolution='l')     # Orthographic map centred on Western Pacific
map.drawcoastlines(linewidth=0.25)
clevs = np.linspace(0,70)
map.contourf(data.LONG, data.LAT, data.LAB, clevs, tri=True, latlon=False)
plt.title("")
plt.show()

However, this yields the following error:

<ipython-input-39-0a19769f6fd0>:5: UserWarning: The following kwargs were not used by contour: 'tri'
  map.contourf(data.LONG, data.LAT, data.LAB, clevs, tri=True, latlon=False)

Is there anyone that could help me? Many thanks in advance, Thomas

Can you show the output of `array`? Matplotlib and Cartopy is not an option? — Benjamin T. Schwertfeger, Mar 29 '23 at 16:25

yvonnefroehlich · Answer 1 · 2023-05-05T06:32:21.070

An input file with different columns (lon, lat, z_value1, z_value2, ...) can be passed directly to pygmt.xyz2grd to create a PyGMT-ready grid. The columns, which should be used for lon, lat, and z-value can be selected via the incols parameter. It is also possible to read the file first into a pandas DataFrame. A conversion to a xarray DataArray is not necessary.

This grid can then be plotted using pygmt.Figure.grdimage. Select the colormap via the cmap parameter.

Maybe this general code example can serve as an orientation.

Input data: Please copy/paste it into a txt file named "test_data.txt", as it seems not possible to upload a txt file.

lon lat z_value1 z_value2
1 2 12 458
5 4 25 265
8 3 24 123
7 4 45 234
2 5 21 401
3 1 78 201
4 7 63 310
2 8 16 520
7 6 23 275

Code example:

import pandas as pd
import pygmt


# Read data in pandas DataFrame
data_df = pd.read_csv("test_data.txt", sep=" ", header=0)

# Make PyGMT-ready grid from data
grid_test = pygmt.xyz2grd(
    #data="test_data.txt",  # pass file (add # at beginning of header line in file) OR
    data=data_df,  # pass pandas DataFrame
    region="d",
    # lon | lat | z-value
    incols="0,1,3",  # input column, zero-based indexing
    spacing="150m",  # arc-minutes, adjust for your needs and data
)


# Create figure object
fig = pygmt.Figure()

fig.grdimage(
    region="d",
    projection="N12c",
    grid=grid_test,
    cmap="batlow",  # select colormap
)

fig.coast(
    shorelines="0.2p,gray",
    frame=True,
)

fig.colorbar()

fig.show()
#fig.savefig(fname="xyz2grd_grdimage.png")

Output figure:

plot xarray Data array in pyGMT

1 Answers1