scipy.griddata出现QhullError(QH6019)跨系统异常求助

在地理空间分析任务中，使用scipy 1.9.3的interpolate.griddata将投影为EPSG25832的geopandas数据插值到规则网格。两台配置完全相同（Python3.9.13、venv）的Linux系统中，一台因依赖系统安装的qhull抛出QhullError(QH6019)，另一台使用scipy自带的qhull则运行正常。

插值函数代码 #

def grid_interp(data, name, xgrid, ygrid):
    '''
    Interpolates point data from geopandas geoseries to a numpy 2D-array of regularly spaced grid points.
    '''
    points = np.vstack((data.geometry.x, data.geometry.y)).T
    values = griddata(
        points, data[name],
        (xgrid, ygrid),
        method=Config.interpolation_method,  # 'linear' and 'cubic' will result in nan outside of the convex hull of data points
    )
    nan_mask = np.isnan(values)  # if there are any nan points re-interpolate them using method 'nearest'

    if np.any(nan_mask):
        values2 = griddata(
            points, data[name],
            (xgrid, ygrid), method='nearest',
        )
        values[nan_mask] = values2[nan_mask]
    return values

网格生成代码 #

xgrid, ygrid = np.meshgrid(np.arange(xmin, xmax + xres, xres), np.arange(ymin, ymax + yres, yres))

完整错误信息 #

---------------------------------------------------------------------------
QhullError                                Traceback (most recent call last)
/home/ts/analysis/analysis_geospatial.ipynb Cell 11 in 8
      3 xgrid, ygrid = np.meshgrid(np.arange(xmin, xmax + xres, xres), 
      4                            np.arange(ymin, ymax + yres, yres),
      5                           )
      7 # target_values = grid_interp(station_data, target, xgrid, ygrid)
----> 8 target_values = grid_interp(station_data, 'MassConcentration', xgrid, ygrid)
      9 sedDBD_values = grid_interp(station_data, 'SedDryBulkDensity', xgrid, ygrid)
     11 target_clipped = grid_clip(target_values, poly, xgrid, ygrid)

/home/ts/analysis/analysis_geospatial.ipynb Cell 11 in 7
      2 '''
      3 Interpolates point data from geopandas geoseries to a numpy 2D-array of regularly spaced grid points.
      4 '''
      6 points = np.vstack((data.geometry.x, data.geometry.y)).T
----> 7 values = griddata(
      8     points, data[name],
      9     (xgrid, ygrid),
     10     method=Config.interpolation_method,  # 'linear' and 'cubic' will result in nan outside of the convex hull of data points
     11     # qhull_options='Qz',
     12 )
     13 nan_mask = np.isnan(values)  # if there are any nan points re-interpolate them using method 'nearest'
     15 if np.any(nan_mask):

File ~/.local/share/virtualenvs/ts/lib/python3.9/site-packages/scipy/interpolate/_ndgriddata.py:260, in griddata(points, values, xi, method, fill_value, rescale)
    258     return ip(xi)
    259 elif method == 'linear':
---> 260     ip = LinearNDInterpolator(points, values, fill_value=fill_value,
    261                               rescale=rescale)
    262     return ip(xi)
    263 elif method == 'cubic' and ndim == 2:

File interpnd.pyx:280, in scipy.interpolate.interpnd.LinearNDInterpolator.__init__()

File _qhull.pyx:1846, in scipy.spatial._qhull.Delaunay.__init__()

File _qhull.pyx:358, 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 Q12 Qt Qbb Qc Qz
Options selected for Qhull 2019.1.r 2019/06/21:
  run-id 559977496  delaunay  Q12-allow-wide  Qtriangulate  Qbbound-last
  Qcoplanar-keep  Qz-infinity-point  _pre-merge  _zero-centrum  Qinterior-keep
  Pgood  _maxoutside  0

解决方案 #

• 强制使用scipy内置qhull：差异源于系统优先调用了系统级qhull，而非scipy自带版本。可通过源码编译安装scipy确保使用内置库：

  1. 卸载虚拟环境中的scipy：pip uninstall scipy
  2. 源码安装指定版本，禁用系统qhull检测：pip install scipy==1.9.3 --no-binary :all:
     （需提前安装编译依赖：gcc、gfortran、liblapack-dev等）

• 显式添加Qhull参数Qz：根据错误提示，在griddata调用中传入qhull_options='Qz'，强制添加无穷远点解决共圆/共球面问题：

  values = griddata(
      points, data[name],
      (xgrid, ygrid),
      method=Config.interpolation_method,
      qhull_options='Qz'
  )
  

• 数据预处理添加微小扰动：打破点的严格共线/共圆性，避免数值问题：

  # 生成points后添加扰动
  points += np.random.normal(0, 1e-6, points.shape)
  

内容的提问来源于stack exchange，提问作者roble