As part of a larger analysis I am interpolating data at points to a grid using scipy.interpolate.griddata. (scipy version 1.9.3) Since recently it worked without problems. Now I am getting QHullError QH6019 on one computer, however not on another system. Both are Linux systems with same Python version (3.9.13), venv setup etc. One difference is that on the system that throws the error qhull is installed from the system repositories, while the other does not have qhull as such (I believe scipy brings its own then).
The function that does the interpolation is:
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
with data being a geopandas dataframe (with geometry projection as EPSG25832), name the column name that holds the values to interpolate, and xgrid and ygrid are the gridpoints resulting from:
xgrid, ygrid = np.meshgrid(np.arange(xmin, xmax + xres, xres), np.arange(ymin, ymax + yres, yres))
I have seen in other question, that qhull has problems with some projections (e.g. here), however, I know that it works with the current projection from the second system where I don't get the error.
The full error message is:
---------------------------------------------------------------------------
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
