SciTools / cartopy

Cartopy - a cartographic python library with matplotlib support
https://scitools.org.uk/cartopy/docs/latest
BSD 3-Clause "New" or "Revised" License
1.43k stars 365 forks source link

Vector transform fails for global data #1608

Open karlwx opened 4 years ago

karlwx commented 4 years ago

Description

When plotting global data, vector transforms fail with many projections with a QhullError. This is similar to an issue in this stackoverflow question: https://stackoverflow.com/questions/50454322/matplotlib-cartopy-streamplot-results-in-qhullerror-with-some-projections Some map projections throw an error (Lambert Conformal, Orthographic (as in example), while others may result in no data being plotted (Plate Carree with regrid_shape specififed).

Code to reproduce

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

extent = [-123,-70,20,55]
fig = plt.figure(figsize=(12,9))
ax = fig.add_subplot(111, projection=ccrs.LambertConformal())
ax.set_extent(extent)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.STATES.with_scale('50m'))

lon = np.arange(0,360,1).astype('float')
lat = np.arange(-90,91,1).astype('float')
lons, lats = np.meshgrid(lon, lat)
u = np.ones(lats.shape) * 10.
v = np.ones(lons.shape) * 10.

ax.barbs(lons,lats,u,v,regrid_shape=12,transform=ccrs.PlateCarree())

Traceback

QhullError                                Traceback (most recent call last)
<ipython-input-51-0c1c6f67d988> in <module>
     13 v = np.ones(lons.shape) * 10.
     14 
---> 15 ax.barbs(lons,lats,u,v,regrid_shape=12,transform=ccrs.PlateCarree())

~/.conda/envs/ewall2/lib/python3.8/site-packages/cartopy/mpl/geoaxes.py in wrapper(self, *args, **kwargs)
    308 
    309         kwargs['transform'] = transform
--> 310         return func(self, *args, **kwargs)
    311     return wrapper
    312 

~/.conda/envs/ewall2/lib/python3.8/site-packages/cartopy/mpl/geoaxes.py in barbs(self, x, y, u, v, *args, **kwargs)
   1902                 args = (c,) + args[1:]
   1903             else:
-> 1904                 x, y, u, v = vector_scalar_to_grid(
   1905                     t, self.projection, regrid_shape, x, y, u, v,
   1906                     target_extent=target_extent)

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

~/.conda/envs/ewall2/lib/python3.8/site-packages/cartopy/vector_transform.py in _interpolate_to_grid(nx, ny, x, y, *scalars, **kwargs)
     65     s_grid_tuple = tuple()
     66     for s in scalars:
---> 67         s_grid_tuple += (griddata(points, s.ravel(), (x_grid, y_grid),
     68                                   method='linear'),)
     69     return (x_grid * xr + x0, y_grid * yr + y0) + s_grid_tuple

~/.conda/envs/ewall2/lib/python3.8/site-packages/scipy/interpolate/ndgriddata.py in griddata(points, values, xi, method, fill_value, rescale)
    216         return ip(xi)
    217     elif method == 'linear':
--> 218         ip = LinearNDInterpolator(points, values, fill_value=fill_value,
    219                                   rescale=rescale)
    220         return ip(xi)

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 (qh_scalelast): can not scale last coordinate to [   0,  inf].  Input is cocircular or cospherical.   Use option 'Qz' to add a point at infinity.

While executing:  | qhull d Qc Qbb Qt Qz Q12
Options selected for Qhull 2019.1.r 2019/06/21:
  run-id 1786751596  delaunay  Qcoplanar-keep  Qbbound-last  Qtriangulate
  Qz-infinity-point  Q12-allow-wide  _pre-merge  _zero-centrum  Qinterior-keep
  Pgood  _maxoutside  0

The above code can work if the latitudes are changed to: np.arange(-89,91,1) I've tried using masked arrays (data outside of visible map extent is masked), but that does not change anything. Also, setting the masked values to NaN does not work - that results in another error.

Is there a way to make this work without taking a subset of the global data?

greglucas commented 4 years ago

The regridding takes place in the target projection (axes), so some of the spatial locations are getting mapped to the infs and messing the interpolation routines up. I was running into similar issues when I opened PR https://github.com/SciTools/cartopy/pull/1636. That PR doesn't currently work on your example with the 0->360 longitudes, but if you change to -180->180 it does. One suggestion in the meantime would be to calculate and transform the regridded points yourself, which I outlined in a few steps in that PR.