Matplotlib-Cartopy Streamplot results in QhullError with some projections

Question

I'd like to plot streamfunctions of global data on an Orthographic projection, but this appears to break in the vector transform. Maybe I'm missing something about the transform keyword that deals with this? I tried with various projections: some worked, many didn't. Is it possible to use streamplot on global data with Orthographic (or similar) projections?

I'm using python 3.6, numpy 1.14.3, xarray 0.10.3, matplotlib 2.2.2, and cartopy 0.16.0.

Here's an example:

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fakelon = np.linspace(-180, 180, 288)
fakelat = np.linspace(-90, 90, 192)
u = xr.DataArray(np.random.rand(len(fakelat), len(fakelon)), coords=[fakelat, fakelon], dims=['lat', 'lon'])
v = xr.DataArray(np.random.rand(len(fakelat), len(fakelon)), coords=[fakelat, fakelon], dims=['lat', 'lon'])
x,y = np.meshgrid(u['lon'], u['lat'])
fig, ax = plt.subplots(subplot_kw={'projection':ccrs.Orthographic()})
ax.set_global()
ax.coastlines()
ax.streamplot(x, y, u.values, v.values, transform=ccrs.PlateCarree())
plt.show()

This results in

~/anaconda/envs/py3_forge/lib/python3.6/site-packages/cartopy/vector_transform.py:138: UserWarning: Some vectors at source domain corners may not have been transformed correctly
  u, v = target_proj.transform_vectors(src_crs, x, y, u, v)
~/anaconda/envs/py3_forge/lib/python3.6/site-packages/cartopy/vector_transform.py:138: RuntimeWarning: invalid value encountered in subtract
  u, v = target_proj.transform_vectors(src_crs, x, y, u, v)
---------------------------------------------------------------------------
QhullError                                Traceback (most recent call last)
<ipython-input-238-9ea7cd02e64e> in <module>()
      8 ax.coastlines()
      9 magnitude = (u ** 2 + v ** 2) ** 0.5
---> 10 ax.streamplot(x, y, u.values, v.values, transform=ccrs.PlateCarree())
     11 plt.show()

~/anaconda/envs/py3_forge/lib/python3.6/site-packages/cartopy/mpl/geoaxes.py in streamplot(self, x, y, u, v, **kwargs)
   1887         gridded = vector_scalar_to_grid(t, self.projection, regrid_shape,
   1888                                         x, y, u, v, *scalars,
-> 1889                                         target_extent=target_extent)
   1890         x, y, u, v = gridded[:4]
   1891         # If scalar fields were regridded then replace the appropriate keyword

~/anaconda/envs/py3_forge/lib/python3.6/site-packages/cartopy/vector_transform.py in vector_scalar_to_grid(src_crs, target_proj, regrid_shape, x, y, u, v, *scalars, **kwargs)
    142     # Now interpolate to a regular grid in projection space, treating each
    143     # component as a scalar field.
--> 144     return _interpolate_to_grid(nx, ny, x, y, u, v, *scalars, **kwargs)

~/anaconda/envs/py3_forge/lib/python3.6/site-packages/cartopy/vector_transform.py in _interpolate_to_grid(nx, ny, x, y, *scalars, **kwargs)
     64     for s in scalars:
     65         s_grid_tuple += (griddata(points, s.ravel(), (x_grid, y_grid),
---> 66                                   method='linear'),)
     67     return (x_grid, y_grid) + s_grid_tuple
     68 

~/anaconda/envs/py3_forge/lib/python3.6/site-packages/scipy/interpolate/ndgriddata.py in griddata(points, values, xi, method, fill_value, rescale)
    220     elif method == 'linear':
    221         ip = LinearNDInterpolator(points, values, fill_value=fill_value,
--> 222                                   rescale=rescale)
    223         return ip(xi)
    224     elif method == 'cubic' and ndim == 2:

interpnd.pyx in scipy.interpolate.interpnd.LinearNDInterpolator.__init__()

qhull.pyx in scipy.spatial.qhull.Delaunay.__init__()

qhull.pyx in scipy.spatial.qhull._Qhull.__init__()

QhullError: QH6019 qhull input error: can not scale last coordinate.  Input is cocircular
   or cospherical.   Use option 'Qz' to add a point at infinity.

While executing:  | qhull d Qbb Q12 Qc Qz Qt
Options selected for Qhull 2015.2.r 2016/01/18:
  run-id 584775470  delaunay  Qbbound-last  Q12-no-wide-dup  Qcoplanar-keep
  Qz-infinity-point  Qtriangulate  _pre-merge  _zero-centrum  Qinterior-keep
  Pgood

swatchai swatchai · Accepted Answer · 2018-05-28T13:31:42

Q: Is it possible to use streamplot on global data with Orthographic (or similar) projections?

A: On Orthographic -- > Yes. Provided that the array of points to generate the plot do not fall out of the visible area of the projection.

Here is the working code to try, and a sample plot.

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# prepare points array that are located on the ...
#  visible part of the orthographic projection
slon = np.linspace(-90., 90., 8)   # no further than +- 90 deg
slat = np.linspace(-45., 45., 6)   # take points within +-45 deg
sx,sy = np.meshgrid(slon, slat)    # create meshgrid

# prep vector data (us,vs) for stream plot
us, vs = np.ones(sx.shape), np.ones(sy.shape)
us = us * (0.5+np.random.normal(0,1))*10
vs = vs * (0.5+np.random.normal(0,1))*10

fig, ax = plt.subplots(subplot_kw={'projection': ccrs.Orthographic()})
fig.set_size_inches([8,8])

ax.set_global()
ax.stock_img()
ax.coastlines(linewidth=0.5)

# do the streamplot
ax.streamplot(sx, sy, us, vs, transform=ccrs.PlateCarree())

plt.show()  # result will be different for each run

The resulting plot:

Matplotlib-Cartopy Streamplot results in QhullError with some projections

2 Answers