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

我想我看不出你的目标范围是什么，但正交投影有时会出现问题，因为它不具有全球代表性。我的意思是;您的数据似乎是全球性的，并且您正试图将所有数据转换为在任何给定时间只能代表地球的一部分的投影，并且它当前代表的特定部分可以更改数据的转换计算。

或许您可以尝试将目标范围限制在投影范围内。

Answer 2

问：是否可以在具有正交（或类似）投影的全球数据上使用 streamplot？

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

结果图：

Matplotlib-Cartopy Streamplot 导致带有一些投影的 QhullError

Matplotlib-Cartopy Streamplot results in QhullError with some projections

python

matplotlib

cartopy