Plot wind vectors from a XArray DataArray

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)

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

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:

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