Plot wind vectors from a XArray DataArray

Q&A PyGMT Q&A
1

Hi!

I am trying to plot temperature and wind data from a Xarray data array.

To show what I want to achieve, I have attached an example Matplotlib plot and a GMT plot that somewhat mimics it, here still without the wind vectors. The plot is made with GMT 6.1.1 and pygmt 0.2.0.

I am getting stuck trying to plot the wind vectors.

How should I proceed in pygmt? I tried:

u=ds.u.isel(isobaricInhPa=0)
v=ds.v.isel(isobaricInhPa=0)
fig.plot(lons.values, lats.values, direction=[u.values, v.values], style='V')

which yields:

    ---------------------------------------------------------------------------
    GMTInvalidInput                           Traceback (most recent call last)
    <ipython-input-174-4dd1ef6c1a08> in <module>
         36 fig.coast(projection=PROJECTION, region=REGION, frame="g", land="gray")
         37 fig.grdimage(windspeed.isel(isobaricInhPa=0), n='b', Q=True)
    ---> 38 fig.plot(lons.values, lats.values, direction=[u.values, v.values], style='V')
         39 fig.coast(projection=PROJECTION, region=REGION, frame="g", shorelines=True)
         40 fig.text(position='BL',text='d.', font=label_font)

    ~/opt/miniconda3/envs/pygmt/lib/python3.6/site-packages/pygmt/helpers/decorators.py in new_module(*args, **kwargs)
        235                 if alias in kwargs:
        236                     kwargs[arg] = kwargs.pop(alias)
    --> 237             return module_func(*args, **kwargs)
        238 
        239         new_module.aliases = aliases

    ~/opt/miniconda3/envs/pygmt/lib/python3.6/site-packages/pygmt/helpers/decorators.py in new_module(*args, **kwargs)
        372                         kwargs[arg] = separators[fmt].join(f"{item}" for item in value)
        373             # Execute the original function and return its output
    --> 374             return module_func(*args, **kwargs)
        375 
        376         return new_module

    ~/opt/miniconda3/envs/pygmt/lib/python3.6/site-packages/pygmt/base_plotting.py in plot(self, x, y, data, sizes, direction, **kwargs)
        561                 )
        562 
    --> 563             with file_context as fname:
        564                 arg_str = " ".join([fname, build_arg_string(kwargs)])
        565                 lib.call_module("plot", arg_str)

    ~/opt/miniconda3/envs/pygmt/lib/python3.6/contextlib.py in __enter__(self)
         79     def __enter__(self):
         80         try:
    ---> 81             return next(self.gen)
         82         except StopIteration:
         83             raise RuntimeError("generator didn't yield") from None

    ~/opt/miniconda3/envs/pygmt/lib/python3.6/site-packages/pygmt/clib/session.py in virtualfile_from_vectors(self, *vectors)
       1160         rows = len(arrays[0])
       1161         if not all(len(i) == rows for i in arrays):
    -> 1162             raise GMTInvalidInput("All arrays must have same size.")
       1163 
       1164         family = "GMT_IS_DATASET|GMT_VIA_VECTOR"

    GMTInvalidInput: All arrays must have same size.

Any help would be greatly appreciated!
Jelle

p.s. the fig.plot() is part of a larger pygmt call, for context:

label_font = '14p,Helvetica-Bold,black'
REGION = 'g'
PROJECTION = 'G-35/60/8c'
# Figuur D
fig.coast(projection=PROJECTION, region=REGION, frame="g", land="gray")
fig.grdimage(windspeed.isel(isobaricInhPa=0), n='b', Q=True)
fig.plot(lons.values, lats.values, direction=[u.values, v.values], style='V')
fig.coast(projection=PROJECTION, region=REGION, frame="g", shorelines=True)
fig.text(position='BL',text='d.', font=label_font)

temperature-wind-pl
temperature-wind-pl1200×800 308 KB

ecmwf_gmt
ecmwf_gmt738×610 499 KB

2

Hi @jdassink, welcome to the GMT forum!

I’m not too familiar with plotting velocity vectors, but the GMTInvalidInput error message you’re getting indicates that your arrays (lons, lats, u and v) are probably of different lengths. Perhaps do a print(len(lons.values), len(lats.values), len(u.values), len(v.values)) and see if they are the same.

If not, maybe upload a sample of your xarray.DataArray (e.g. as a zipped NetCDF file) so that we can take a look at what’s wrong?

As an aside, if you want to get the pointy colorbar like in the matplotlib plot, do fig.colorbar(position="JBC+e") (key is the +e which adds sidebar triangles). See also https://docs.generic-mapping-tools.org/latest/colorbar.html#d.

P.S. There’s a PR in progress for velo at https://github.com/GenericMappingTools/pygmt/pull/525 which may be of interest if you need to plot more complicated velocity vectors (e.g. with error ellipses and annotations).

3

Hi @weiji14, thanks for your answer.

I tried a bit more and decided to use grdvector by using the GMT clib directly.
I advanced a bit in my figure, and now am almost able to produce what I want. Basically the call that is use is:

fig.coast(projection=PROJECTION, region=REGION, frame="g", land="gray")
fig.grdimage(windspeed.isel(isobaricInhPa=0), n='b', Q=True)
fig.coast(projection=PROJECTION, region=REGION, frame="g", shorelines=True)
ds.u.isel(isobaricInhPa=0).to_netcdf('u.grd')
ds.v.isel(isobaricInhPa=0).to_netcdf('v.grd')
with pygmt.clib.Session() as session:
    session.call_module('grdvector', 
                        f'u.grd v.grd -I10d/10d -Gwhite -Q0.1i+e+gwhite -W0.8p,black -Si0.1d')
fig.text(position='BL',text='d.', font=label_font)

and so forth for the other subfigures.

I get the following plot:

gmt_20201112
gmt_20201112738×623 239 KB

It shows these weird effects around the poles. Any idea?

4

Hi!

I have been having the same error message (GMTInvalidInput: All arrays must have same size). I checked the length of each array and they are equal:

image

Here is what I try to do:

fig = pygmt.Figure()
fig.basemap(region=region, projection=proj, frame=True)
fig.grdimage("ETOPO1_Ice_g_geotiff.tif", projection=proj, frame=True, cmap="oleron")
fig.colorbar(frame='af+l"Topography (m)"')
fig.plot(
    x=lon,
    y=lat,
    direction=[[azi],[mag]],
    style="V",
    pen="1p",
    color="black",
)
fig.savefig("tectonics_map.png")
fig.show()

Many thanks!

5

Despite my PyIgnorance I strongly suspect of the above line

6

Hi Joaquim! Would you have any suggestion how I should code this instead?

Thanks!
Thomas

7

Sure, remove that trouble-maker " : " from file.

EDIT: hoops, we are mixing the threads.

8

Hmm, looking at https://www.pygmt.org/v0.6.1/api/generated/pygmt.Figure.plot.html, I’m guessing that you should do direction=[azi, mag] instead? What type of Python objects are azi and mag?

9

Fixed it! I needed to plot it as direction=[azi, mag] indeed. Also, my angles were in radians, not degrees, so that’s why I kept running into the problem of only north-facing vectors when it did work :sweat_smile:

1 Like