Matplotlib-Cartopy Streamplot приводит к QhullError с некоторыми проекциями

Question

Я хотел бы построить функции потока глобальных данных на ортогональной проекции, но, похоже, это нарушает векторное преобразование.Может быть, я что-то упускаю из ключевого слова transform, которое имеет к этому отношение?Я пробовал разные прогнозы: некоторые работали, многие нет.Можно ли использовать streamplot для глобальных данных с ортогональными (или аналогичными) проекциями?

Я использую python 3.6, numpy 1.14.3, xarray 0.10.3, matplotlib 2.2.2 и cartopy 0.16.0.

Вот пример:

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

В результате получается

~/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

Answer 1

В: Можно ли использовать стримплот на глобальных данных с ортогональными (или аналогичными) проекциями?

A: На орфографических -> Да.При условии, что массив точек для построения графика не выпадает из видимой области проекции.

Вот рабочий код, который нужно попробовать, и пример графика.

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

Полученный участок:

Answer 2

Я не думаю, что могу видеть, какова ваша целевая протяженность, но у Орфографической проекции иногда могут быть проблемы, поскольку она не является глобально репрезентативной.Я имею в виду;ваши данные кажутся глобальными, и вы пытаетесь преобразовать все это в проекцию, которая может представлять только часть земного шара в любой момент времени, и конкретная часть, которую она представляет в настоящее время, может изменить вычисления преобразования данных.

Возможно, вы можете попытаться ограничить целевые экстенты, чтобы они попадали в пределы проекции.